text:
joint_types = [
    :Fixed,
    :Prismatic,
    :Planar,
    :FixedOrientation,
    :Revolute,
    :Cylindrical,
    :PlanarAxis,
    :FreeRevolute,
    :Orbital,
    :PrismaticOrbital,
    :PlanarOrbital,
    :FreeOrbital,
    :Spherical,
    :CylindricalFree,
    :PlanarFree
    ]

function test_get_set_data(mechanism::Mechanism)
    Nd = Dojo.data_dim(mechanism,
		attjac=false)
    data0 = rand(Nd)
    Dojo.set_data!(mechanism, data0)
    data1 = Dojo.get_data(mechanism)
    @test norm(data0 - data1) < 1.0e-8
end

@testset "Get and set data" begin
    mech = Dojo.get_snake(num_bodies=3,
		damper=1.0,
			spring=1.0,
			contact_type=:nonlinear);
    test_get_set_data(mech)

    mech = Dojo.get_snake(num_bodies=3,
		damper=1.0,
		spring=1.0,
		contact_type=:linear);
    test_get_set_data(mech)

    mech = Dojo.get_snake(num_bodies=3,
		damper=1.0,
		spring=1.0,
		contact_type=:impact);
    test_get_set_data(mech)

    mech = Dojo.get_pendulum(
		damper=1.0,
		spring=10.0);
    test_get_set_data(mech)

    mech = Dojo.get_humanoid(
		damper=1.0,
		spring=10.0,
		contact_feet=true);
    test_get_set_data(mech)

    mech = Dojo.get_humanoid(
		damper=1.0,
		spring=10.0,
		contact_feet=false);
    test_get_set_data(mech)

    mech = Dojo.get_atlas(
		damper=1.0,
		spring=10.0);
    test_get_set_data(mech)

    mech = Dojo.get_quadruped(
		damper=1.0,
		spring=10.0);
    test_get_set_data(mech)
end

################################################################################
# Analytical Jacobian
################################################################################
# Controller
function ctrl!(mechanism, k)
	nu = Dojo.input_dimension(mechanism)
	if Dojo.input_dimension(mechanism.joints[1]) == 6
		u = 0.2 * [szeros(6); mechanism.timestep * sones(nu-6)]
	else
		u = 0.2 * mechanism.timestep * sones(nu)
	end
	Dojo.set_input!(mechanism, u)
	return
end

function test_data_system(model::Symbol;
		ϵ=1.0e-6,
		tsim=0.1,
		ctrl=(m, k)->nothing,
        timestep=0.01,
		gravity=[0.0; 0.0; -9.81],
		verbose=false,
		T=Float64,
		kwargs...)

    # mechanism
    mechanism = Dojo.get_mechanism(model,
		timestep=timestep,
		gravity=gravity;
		kwargs...)
    Dojo.initialize!(mechanism, model)
    # simulate
    Dojo.simulate!(mechanism, tsim, ctrl!,
        record=false,
		verbose=false,
		opts=Dojo.SolverOptions(rtol=ϵ, btol=ϵ))

	# Finite Difference
	Nd = Dojo.data_dim(mechanism,
		attjac=false)
	data0 = Dojo.get_data(mechanism)
	sol0 = Dojo.get_solution(mechanism)
	datajac0 = Dojo.finite_difference_data_jacobian(mechanism, data0, sol0)
	attjac0 = Dojo.data_attitude_jacobian(mechanism)
	datajac0 *= attjac0

	# Analytical
	D = Dojo.create_data_matrix(mechanism.joints, mechanism.bodies, mechanism.contacts)
	Dojo.jacobian_data!(D, mechanism)
	nodes = [mechanism.joints; mechanism.bodies; mechanism.contacts]
	dimrow = length.(nodes)
	dimcol = Dojo.data_dim.(nodes)
	datajac1 = Dojo.full_matrix(D, dimrow, dimcol)

	# Test
	@testset "$(String(model))" begin
		@test norm(datajac0 - datajac1, Inf) < 1.0e-6
	end
    return nothing
end

################################################################################
# Without contact and joint limits
################################################################################
@testset "Data Jacobian without contact and limtis" begin
	for (spring, damper) in [(0.0, 0.0), (2.0, 0.3)]
		test_data_system(:sphere,
			contact=false)
		test_data_system(:box,
			contact=false)
		test_data_system(:box2D,
			contact=false)
		test_data_system(:slider,
			spring=spring,
			damper=damper)
		test_data_system(:nslider,
			spring=spring,
			damper=damper)
		test_data_system(:pendulum,
			spring=spring,
			damper=damper)
		test_data_system(:cartpole,
			spring=spring,
			damper=damper)
		test_data_system(:pendulum,
			spring=spring,
			damper=damper)
		test_data_system(:hopper,
			spring=spring,
			damper=damper,
			contact_foot=false,
			contact_body=false)
		test_data_system(:humanoid,
			spring=spring,
			damper=damper,
			contact_feet=false,
			contact_body=false)
		test_data_system(:atlas,
			spring=spring,
			damper=damper,
			contact_feet=false,
			contact_body=false)
		test_data_system(:halfcheetah,
			contact_feet=false,
			contact_body=false,
			limits=false)
		test_data_system(:walker, spring=spring,
			damper=damper,
			contact_feet=false,
			contact_body=false,
			limits=false)
		test_data_system(:quadruped, spring=spring,
			damper=damper, contact_feet=false,
			contact_body=false,
			limits=false)
		for joint_type in joint_types
			test_data_system(:snake,
				num_bodies=5,
				spring=spring,
				damper=damper,
				contact=false,
				joint_type=joint_type)
			test_data_system(:twister,
				num_bodies=5,
				spring=spring,
				damper=damper,
				contact=false,
				joint_type=joint_type)
		end
	end
end

################################################################################
# With contact and joint limits
################################################################################
@testset "Data Jacobian with contact and limtis" begin
	for (spring, damper) in [(0.0, 0.0), (2.0, 0.3)]
		test_data_system(:sphere,
			contact=true)
		test_data_system(:box,
			contact=true)
		test_data_system(:box2D,
			contact=true)
		test_data_system(:slider,
			spring=spring,
			damper=damper)
		test_data_system(:nslider,
			spring=spring,
			damper=damper)
		test_data_system(:pendulum,
			spring=spring,
			damper=damper)
		test_data_system(:cartpole,
			spring=spring,
			damper=damper)
		test_data_system(:pendulum,
			spring=spring,
			damper=damper)
		test_data_system(:hopper,
			spring=spring,
			damper=damper,
			contact_foot=true,
			contact_body=true)
		test_data_system(:humanoid,
			spring=spring,
			damper=damper,
			contact_feet=true,
			contact_body=true)
		test_data_system(:atlas,
			spring=spring,
			damper=damper,
			contact_feet=true,
			contact_body=true)
		test_data_system(:halfcheetah,
			contact_feet=true,
			contact_body=true,
			limits=true)
		test_data_system(:walker,
			spring=spring,
			damper=damper,
			contact_feet=true,
			contact_body=true,
			limits=true)
		test_data_system(:quadruped,
			spring=spring,
			damper=damper,
			contact_feet=true,
			contact_body=true,
			limits=true)
		for joint_type in joint_types
			test_data_system(:snake,
				num_bodies=5,
				spring=spring,
				damper=damper,
				contact=true,
				joint_type=joint_type)
			test_data_system(:twister,
				num_bodies=5,
				spring=spring,
				damper=damper,
				contact=true,
				joint_type=joint_type)
		end
	end
end


text:
# An implementation of CGLS for the solution of the
# over-determined linear least-squares problem
#
#  minimize ‖Ax - b‖
#
# equivalently, of the normal equations
#
#  A'Ax = A'b.
#
# CGLS is formally equivalent to applying the conjugate gradient method
# to the normal equations but should be more stable. It is also formally
# equivalent to LSQR though LSQR should be expected to be more stable on
# ill-conditioned or poorly scaled problems.
#
# This implementation is the standard formulation, as recommended by
# A. Björck, T. Elfving and Z. Strakos, Stability of Conjugate Gradient
# and Lanczos Methods for Linear Least Squares Problems.
#
# Dominique Orban, <dominique.orban@gerad.ca>
# Princeton, NJ, March 2015.

export cgls

"""Solve the regularized linear least-squares problem

  minimize ‖b - Ax‖₂² + λ ‖x‖₂²

using the Conjugate Gradient (CG) method, where λ ≥ 0 is a regularization
parameter. This method is equivalent to applying CG to the normal equations

  (A'A + λI) x = A'b

but is more stable.

CGLS produces monotonic residuals ‖r‖₂ but not optimality residuals ‖A'r‖₂.
It is formally equivalent to LSQR, though can be slightly less accurate,
but simpler to implement.
"""
function cgls(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(b,1)), λ :: Float64=0.0,
              atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6, radius :: Float64=0.0,
              itmax :: Int=0, verbose :: Bool=false) where T <: Number

  m, n = size(A);
  size(b, 1) == m || error("Inconsistent problem size");
  verbose && @printf("CGLS: system of %d equations in %d variables\n", m, n);

  x = zeros(T, n);
  r = copy(b)
  bNorm = @knrm2(m, r)   # Marginally faster than norm(b);
  bNorm == 0 && return x, SimpleStats(true, false, [0.0], [0.0], "x = 0 is a zero-residual solution");
  s = A' * M * r;
  p = copy(s);
  γ = @kdot(n, s, s)  # Faster than γ = dot(s, s);
  iter = 0;
  itmax == 0 && (itmax = m + n);

  rNorm  = bNorm;
  ArNorm = sqrt(γ);
  rNorms = [rNorm;];
  ArNorms = [ArNorm;];
  ε = atol + rtol * ArNorm;
  verbose && @printf("%5s  %8s  %8s\n", "Aprod", "‖A'r‖", "‖r‖")
  verbose && @printf("%5d  %8.2e  %8.2e\n", 1, ArNorm, rNorm);

  status = "unknown";
  on_boundary = false
  solved = ArNorm <= ε;
  tired = iter >= itmax;

  while ! (solved || tired)
    q = A * p;
    δ = @kdot(m, q, M * q)   # Faster than α = γ / dot(q, q);
    λ > 0 && (δ += λ * @kdot(n, p, p))
    α = γ / δ;

    # if a trust-region constraint is give, compute step to the boundary
    σ = radius > 0.0 ? maximum(to_boundary(x, p, radius)) : α
    if (radius > 0.0) & (α > σ)
      α = σ
      on_boundary = true
    end

    @kaxpy!(n,  α, p, x)     # Faster than x = x + α * p;
    @kaxpy!(m, -α, q, r)     # Faster than r = r - α * q;
    s = A' * M * r;
    λ > 0 && @kaxpy!(n, -λ, x, s)   # s = A' * r - λ * x;
    γ_next = @kdot(n, s, s)  # Faster than γ_next = dot(s, s);
    β = γ_next / γ;
    @kscal!(n, β, p)
    @kaxpy!(n, 1.0, s, p)    # Faster than p = s + β * p;
    # The combined BLAS calls tend to trigger some gc.
    #  BLAS.axpy!(n, 1.0, s, 1, BLAS.scal!(n, β, p, 1), 1);
    γ = γ_next;
    rNorm = @knrm2(m, r)  # Marginally faster than norm(r);
    ArNorm = sqrt(γ);
    push!(rNorms, rNorm);
    push!(ArNorms, ArNorm);
    iter = iter + 1;
    verbose && @printf("%5d  %8.2e  %8.2e\n", 1 + 2 * iter, ArNorm, rNorm);
    solved = (ArNorm <= ε) | on_boundary
    tired = iter >= itmax;
  end

  status = on_boundary ? "on trust-region boundary" : (tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol")
  stats = SimpleStats(solved, false, rNorms, ArNorms, status);
  return (x, stats);
end

ext: jl


text:
m = test_page_model()
load_parameters(m)


text:
"""
	csd2by1(Q1, Q2) computes the the Cosine-sine decomposition of
	an (m+l) by l orthogonal matrix Q = [Q1; Q2] where m <= l, such that:

```math
	Q1 = U * C * Z',    Q2 = V * S * Z'
```

- U is m-by-m orthogonal matrix,
- V is l-by-l orthogonal matrix,
- Z is l-by-l orthogonal matrix,
- C and S are  m-by-l and l-by-l ``diagonal" matrices such that
C' * C + S' * S = I. C and S have the following structures
(both alpha and beta are of length l):

(1) If m == l:
        C =  diag(alpha(1), ... , alpha(l))

        S =  diag(beta(1), ... , beta(l))

(2) If m < l:
        C =                       m                 l-m
              m ( diag(alpha(1), ... , alpha(m))     0 )

        S =                       m                 l-m
              m ( diag(beta(1), ... , beta(m)))      0 )
            l-m (                 0                  I )

# arguments:
- Q1: m by l matrix,
- Q2: l by l matrix.

# returns:
- U is m-by-m matrix,
- V is l-by-l matrix,
- Z is l-by-l matrix,
- cosine is an array of length l,
- sine is an array of length l.
"""
# csd2by1() calls csd2by1!() before making a copy of input matrices A and B,
# thus A, B will not be overwritten
function csd2by1(Q1::StridedMatrix{T}, Q2::StridedMatrix{T}) where T
	return csd2by1!(copy(Q1),copy(Q2))
end

function csd2by1!(Q1::StridedMatrix{T}, Q2::StridedMatrix{T}) where T
    m, l = size(Q1)
	if size(Q2)[1] != size(Q2)[2]
		throw(ArgumentError("Q2's row number must equal its column number"))
	elseif m > size(Q2)[1]
		throw(ArgumentError("Q1's row number must be less than that of Q2"))
	elseif l != size(Q2)[2]
		throw(ArgumentError("Q1's column number must equal that of Q2"))
	end

    # Step 1: compute the SVD of Q2
	# May cause LAPACK Exception: S/DBDSDC did not converge
	D1 = svd!(Q2)

	# Step 2.1: rearrange the values of S in non-decreasing order
	sine = reverse(D1.S)

    # Step 2.2: reverse V's col
	#     Assign array slicing to D1.U will throw an error:
	#     D1.U = @views D1.U[:, l:-1:1]
	#     ERROR: setfield! immutable struct of type SVD cannot be changed
	D1.U[:,1:l] = @views D1.U[:, l:-1:1]

    # Step 2.3: reverse Z's col
    D1.V[:,1:l] = @views D1.V[:,l:-1:1]

    # Step 3: find r such that sine(r) <= 1/sqrt(2) < sine(r+1),
	# binary search is used
	r = binSearch(sine)

    # Step 4: compute X = Q1Z
    X = Q1 * D1.V

    # Step 5: compute QR decomposition of X
    D2 = qr!(X)

    # Step 5a: sanitize R so that it only has non-negative diagonal entries
	R = D2.R
    dia1 = sign.(diag(R))
	# update dia1 if R has zeros in the diagonal
	dia1[dia1 .== T(0.0)] .= T(1.0)
	cosine = zeros(T, l)
	@views cosine[1:r] = abs.(diag(R)[1:r])
    R = Diagonal(dia1) * R
	U = Matrix{T}(D2.Q) * Diagonal(dia1)

	# Special case: exit
	if m <= r
		return U, D1.U, D1.V, cosine, sine
	end

	# Step 6: refine R
	R22 = @views R[r+1:m,r+1:l]

	# Step 6a: compute the full SVD of R22
	D3 = svd!(R22, full = true)

	# Step 6b: update the (r+1)-th to m-th cols of U
	U[:,r+1:m] = @views U[:,r+1:m]*D3.U

	# Step 6c: update the last (l-r) cols of Z
	@views D1.V[:,r+1:l] = D1.V[:,r+1:l] * D3.V

	# Step 6d: update cosine and sine
	@views sine[m+1:l] = ones(T, l-m)
	@views cosine[r+1:m] = D3.S

	# Step 7: refine V
    # Step 7a: set W
	W = @views Diagonal(sine[r+1:l])* D3.V

    # Step 7b: compute QR of W
    D4 = qr!(W)
    Qw = D4.Q * Matrix{T}(I, l-r, l-r)

    # Step 7c: sanitize D4.R so that it only has non-negative diagonal entries
    dia2 = sign.(diag(D4.R))
    Qw = Qw * Diagonal(dia2)

    # Step 7d: update V
    @views D1.U[:,r+1:l] = D1.U[:,r+1:l] * Qw

    return U, D1.U, D1.V, cosine, sine
end

# binSearch() modifies standard binary search
# to find r such that beta(r) <= 1/sqrt(2) < beta(r+1),
# return 1 if no such r exist.
function binSearch(a::Vector{T}) where T
	r = 1
	low, high = 1, length(a)
	thld = sqrt(T(0.5))
	while low <= high
		mid = floor(Integer, (high - low)/2) + low
		r = mid
		if mid+1 <= length(a) && a[mid+1] <= thld
			low = mid + 1
		elseif a[mid] > thld
			high = mid - 1
		else
			return mid
		end
	end
	return r
end


text:
"""
    lucky_image(cube; dims, q, metric=:max, register=:dft, kwargs...)

Perform classic lucky imaging for a cube of data along dimension `dims`. The selection fraction is set by the quantile `q`, where `q=0.9` corresponds to a 10% selectrion fraction. The following options are available for modifying the metric used for selection as well as how the frames are registered.

# Metrics

The following metrics, set with the keyword argument `metric` are available

* `:max` - the maximum value
* `:mean` - the mean value
* `:min` - the minimum value (useful for coronagraphic images)
* other - pass a function with the signature `f(::AbstractMatrix)` which returns a scalar value

# Registration

The following options for registering the selected frames are available, set with the `register` keyword argument. Extra options can be provided to `kwargs...`

* `:dft` - use `SubpixelRegistration.jl` to register the frames using their phase offsets. The reference frame will be the one with the highest metric, and keyword arguments like `upsample_factor` can be passed directly.
* `:max` - register to maximum value
* `:com` - register to center of mass
* `:invcom` - register to center of inverse mass (useful for coronagraphic images)
* `nothing` - cube is already co-registered

# Keyword arguments

* `dims` - the dimension along which to perform lucky imaging (required)
* `q` - the selection quantile (required)
* `metric` - the metric used for selection. Default is `:max`
* `register` - the method used for registration. Default is `:dft`
* `window` - if provided, will measure the metric and frame offsets inside a centered window with `window` side length.
* `kwargs...` - additional keyword arguments will be passed to the `register` method (e.g., `upsample_factor`)

# Examples

```julia
julia> cube = # load data ...

julia> res = lucky_image(cube; dims=3, q=0.9, upsample_factor=10)
```

# See Also
[`lucky_image!`](@ref)
"""
function lucky_image(cube::AbstractArray{T}; dims, kwargs...) where T
    pixshape = nottuple(size(cube), dims)
    storage = similar(cube, float(T), pixshape)
    return lucky_image!(storage, cube; dims, kwargs...)
end

"""
    lucky_image!(out::AbstractMatrix, cube; dims, q, kwargs...)

Perform classic lucky imaging and store the combined frame in `out`. See [`lucky_image`](@ref) for a full description.
"""
function lucky_image!(
    storage::AbstractMatrix, cube::AbstractArray{T,3};
    dims, q, minimize=false, center=center(cube), window=nothing,
    metric=:l2norm, register=:dft, shift_reference=true, kwargs...) where T
    # first, get window view of cube if preferred
    if isnothing(window)
        _cube = cube
        offset = (0.0, 0.0)
    else
        _cube = window_view(cube, dims, window, center)
        # get shift frame (possible) window to original frame center
        offset = LuckyImaging.center(cube)[1:2] .- center[1:2]
    end
    # get metric measured in each frame
    frame_gen = eachslice(_cube; dims=dims)
    if metric === :max
        _metric = map(maximum, frame_gen)
    elseif metric === :l2norm
        _metric = map(X -> mean(v -> v^2, X), frame_gen)
    elseif metric === :l1norm
        _metric = map(mean, frame_gen)
    else
        _metric = map(metric, frame_gen)
    end
    # find cutoff based off quantile
    if minimize
        cutoff = quantile(_metric, 1 - q)
        pred = ≤(cutoff)
    else
        cutoff = quantile(_metric, q)
        pred = ≥(cutoff)
    end

    if register === :dft
        # center reference with max/min
        if minimize
            reference = selectdim(_cube, dims, argmin(_metric))
            index = argmin(reference).I
        else
            reference = selectdim(_cube, dims, argmax(_metric))
            index = argmax(reference).I
        end
        delta = LuckyImaging.center(reference) .- index .+ offset
        refshift = shift_reference ? delta : zero.(delta)
    end

    # make sure array is zerod out
    fill!(storage, zero(T))
    norm_value = count(pred, _metric)
    # now, time to go through cube with frames above the cutoff
    # and combine them into the final frame.
    @inbounds for didx in axes(_cube, dims)
        # skip frames below cutoff
        pred(_metric[didx]) || continue
        # get frame
        frame = selectdim(_cube, dims, didx)
        # get index using registration method
        if !isnothing(register)
            if register === :dft
                shift = phase_offset(reference, frame; kwargs...).shift
                shift = shift .+ refshift
            elseif register === :max
                # measure shift relative to subframe
                index = argmax(frame; kwargs...).I
                shift = LuckyImaging.center(frame) .- index .+ offset
            elseif register === :com
                index = center_of_mass(frame; kwargs...)
                shift = LuckyImaging.center(frame) .- index .+ offset
            end
            # apply shift to full frame
            full_frame = selectdim(cube, dims, didx)
            registered_frame = fourier_shift(full_frame, shift)
        else
            registered_frame = selectdim(cube, dims, didx)
        end
        @. storage += registered_frame / norm_value
    end
    return storage
end


text: settings_default = @with_kw (yvar = :log_dailywages, first_idvar = :first_id, second_idvar = :second_id)


text:
type Foam_Data_Path <: Models_Data_Path
  dir_input                ::    ASCIIString
  filename                 ::    ASCIIString
  start_on_month           ::    Int64
  last_on_month            ::    Int64
end

function AFoamProblem (   ;dir_input :: ASCIIString = "/home/3Tstorage/foam/run_sd1/storage/",
                           filename :: ASCIIString = "data/control_sd1.nc",
                           first_day :: Int64 = 0, last_day :: Int64 = 0)
  return Foam_Data_Path(dir_input,filename, first_day , last_day );
end

export Foam_Data_Path
export AFoamProblem


text:
#=
Case folding for Unicode Chr types

Copyright 2017-2020 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
=#

@inline _can_upper_lat(c) = ifelse(c > 0xde, c != 0xf7, c == 0xb5)

_wide_lower_latin(ch) = (ch == 0xb5) | (ch == 0xff) | (ch == 0xdf)

_wide_upper(ch) =
    ifelse(ch == 0xb5, 0x39c,
           ifelse(ch == 0xff, 0x178, ifelse(ch == 0xdf, 0x1e9e, ch%UInt16)))

_lowercase_l(ch) = _is_upper_al(ch)  ? (ch + 0x20) : ch
_uppercase_l(ch) = _can_upper_al(ch) ? (ch - 0x20) : _wide_upper(ch)

lowercase(ch::T) where {T<:Chr} = T(_lowercase(codepoint(ch)))
uppercase(ch::T) where {T<:Chr} = T(_uppercase(codepoint(ch)))
titlecase(ch::T) where {T<:Chr} = T(_titlecase(codepoint(ch)))

lowercase(ch::ASCIIChr) = _is_upper_a(ch) ? ASCIIChr(ch + 0x20) : ch
uppercase(ch::ASCIIChr) = _is_lower_a(ch) ? ASCIIChr(ch - 0x20) : ch
titlecase(ch::ASCIIChr) = uppercase(ch)

lowercase(ch::T) where {T<:LatinChars} = T(_lowercase_l(codepoint(ch)))

_uppercase_latin(ch) = _can_upper_al(ch) ? (ch - 0x20) : ch
uppercase(ch::LatinChr) = LatinChr(_uppercase_latin(codepoint(ch)))

# Special handling for case where this is just an optimization of the first 256 bytes of Unicode,
# and not the 8-bit ISO 8859-1 character set
function uppercase(ch::_LatinChr)
    cb = codepoint(ch)
    _can_upper_al(cb) && return _LatinChr(cb - 0x20)
    # We didn't used to uppercase 0xdf, the ß character, now we do
    cb == 0xdf ? UCS2Chr(0x1e9e) : cb == 0xb5 ? UCS2Chr(0x39c) : cb == 0xff ? UCS2Chr(0x178) : ch
end
titlecase(ch::LatinChars) = uppercase(ch)

@inline function _check_tab(mask, tab, ch)
    t = (ch >>> 9)
    ((mask >>> (t & 0x7f)) & 1) != 0 && (off = tab[t+1]) != 0 &&
        (CaseTables.bitvec[off][((ch >>> 5) & 0xf) + 1] & (UInt32(1) << (ch & 0x1f))) != 0
end

@inline _get_tab(off, ch, base) =
    off == 0 ? ch : (off = CaseTables.offvec[off][((ch >>> 5) & 0x1f) + 1]) == 0 ? ch :
    (base + CaseTables.tupvec[off][(ch & 0x1f) + 1])

@inline _get_tab_bmp(mask, tab, ch) =
    (t = (ch >>> 9); ((mask >>> t) & 1) == 0 ? ch : _get_tab(tab[(t>>1)+1], ch, 0x0000))
@inline _get_tab_slp(mask, tab, ch) =
    (t = (ch >>> 9); ((mask >>> (t & 0x7f)) & 1) == 0 ? ch : _get_tab(tab[(t>>1)+1], ch, 0x10000))

@inline _upper_lat(ch) = _get_tab(CaseTables.ct.u_tab[1], ch, 0x0000)

@inline _upper_bmp(ch) = _get_tab_bmp(CaseTables.ct.can_u_flg,  CaseTables.ct.u_tab, ch)
@inline _lower_bmp(ch) = _get_tab_bmp(CaseTables.ct.can_l_flg,  CaseTables.ct.l_tab, ch)
@inline _title_bmp(ch) = _get_tab_bmp(CaseTables.ct.can_u_flg,  CaseTables.ct.t_tab, ch)
@inline _upper_slp(ch) = _get_tab_slp(CaseTables.ct.can_su_flg, CaseTables.ct.u_tab, ch)
@inline _lower_slp(ch) = _get_tab_slp(CaseTables.ct.can_sl_flg, CaseTables.ct.l_tab, ch)

@inline _can_lower_bmp(ch) = _check_tab(CaseTables.ct.can_l_flg,  CaseTables.ct.can_l_tab, ch)
@inline _can_upper_bmp(ch) = _check_tab(CaseTables.ct.can_u_flg,  CaseTables.ct.can_u_tab, ch)
@inline _can_lower_slp(ch) = _check_tab(CaseTables.ct.can_sl_flg, CaseTables.ct.can_l_tab, ch)
@inline _can_upper_slp(ch) = _check_tab(CaseTables.ct.can_su_flg, CaseTables.ct.can_u_tab, ch)
@inline _is_lower_bmp(ch)  = _check_tab(CaseTables.ct.is_l_flg,   CaseTables.ct.is_l_tab,  ch)
@inline _is_upper_bmp(ch)  = _check_tab(CaseTables.ct.is_u_flg,   CaseTables.ct.is_u_tab,  ch)
@inline _is_lower_slp(ch)  = _check_tab(CaseTables.ct.is_sl_flg,  CaseTables.ct.is_l_tab,  ch)
@inline _is_upper_slp(ch)  = _check_tab(CaseTables.ct.is_su_flg,  CaseTables.ct.is_u_tab,  ch)

@inline _is_lower_ch(ch) =
    ch <= 0x7f ? _is_lower_a(ch) :
    ch <= 0xff ? _is_lower_l(ch) :
    ch <= 0xffff ? _is_lower_bmp(ch) :
    ch <= 0x1ffff ? _is_lower_slp(ch) : false

@inline _is_upper_ch(ch) =
    ch <= 0x7f ? _is_upper_a(ch) :
    ch <= 0xff ? _is_upper_l(ch) :
    ch <= 0xffff ? _is_upper_bmp(ch) :
    ch <= 0x1ffff ? _is_upper_slp(ch) : false

@inline _can_lower_ch(ch) =
    ch <= 0x7f ? _is_upper_a(ch) :
    ch <= 0xff ? _is_upper_l(ch) :
    ch <= 0xffff ? _can_lower_bmp(ch) :
    ch <= 0x1ffff ? _can_lower_slp(ch) : false

@inline _can_upper_ch(ch) =
    ch <= 0x7f ? _is_lower_a(ch) :
    ch <= 0xff ? _can_upper_lat(ch) :
    ch <= 0xffff ? _can_upper_bmp(ch) :
    ch <= 0x1ffff ? _can_upper_slp(ch) : false

@inline _lowercase(ch) =
    ch <= 0x7f ? (ch + (_is_upper_a(ch)<<5)) :
    ch <= 0xff ? (ch + (_is_upper_l(ch)<<5)) :
    ch <= 0xffff ? _lower_bmp(ch) :
    ch <= 0x1ffff ? _lower_slp(ch) : ch

@inline _uppercase(ch) =
    ch <= 0x7f ? (_is_lower_a(ch) ? (ch - 0x20) : ch) :
    ch <= 0xff ? _upper_lat(ch) :
    ch <= 0xffff ? _upper_bmp(ch) :
    ch <= 0x1ffff ? _upper_slp(ch) : ch

@inline _titlecase(ch) =
    ch <= 0x7f ? (_is_lower_a(ch) ? (ch - 0x20) : ch) :
    ch <= 0xff ? _upper_lat(ch) :
    ch <= 0xffff ? _title_bmp(ch) :
    ch <= 0x1ffff ? _upper_slp(ch) : ch

# 0xb5, 0xdf, and 0xff cannot be uppercased in LatinCSE, although they are lowercase
@inline _can_upper_l(c) = (0xe0 <= c <= 0xfe) & (c != 0xf7)
@inline _can_upper_al(c) = _is_lower_a(c) | _can_upper_l(c)

@api develop _is_lower_ch, _is_upper_ch, _can_lower_ch, _can_upper_ch,
             _upper_bmp, _lower_bmp, _title_bmp, _upper_slp, _lower_slp,
             _can_upper_al, _can_upper_lat, _can_lower_bmp, _can_upper_bmp,
             _can_lower_slp, _can_upper_slp


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

# This file extends JuMP to complementarity constraints. It is a good example of
# how JuMP can be extended. In particular, it allows constraints of the form:
#
# @constraint(model, complements(2x - 1, x))
# @constraint(model, 2x - 1 ⟂ x)

function _build_complements_constraint(
    errorf::Function,
    F::AbstractArray{<:AbstractJuMPScalar},
    x::AbstractArray{<:AbstractVariableRef},
)
    if size(F) != size(x)
        errorf(
            "size of mapping does not match size of variables: " *
            "$(size(F)) != $(size(x)).",
        )
    end
    return VectorConstraint(vec([F x]), MOI.Complements(length(F)))
end

function _build_complements_constraint(
    errorf::Function,
    F::Containers.SparseAxisArray{<:AbstractJuMPScalar},
    x::Containers.SparseAxisArray{<:AbstractVariableRef},
)
    elements = [F[i] for i in eachindex(F)]
    for i in eachindex(F)
        if haskey(x, i)
            push!(elements, x[i])
        else
            errorf("keys of the SparseAxisArray's do not match.")
        end
    end
    return VectorConstraint(elements, MOI.Complements(length(F)))
end

function _build_complements_constraint(
    errorf::Function,
    ::AbstractArray{<:AbstractJuMPScalar},
    ::AbstractArray{<:AbstractJuMPScalar},
)
    return errorf("second term must be an array of variables.")
end

function _build_complements_constraint(
    ::Function,
    F::AbstractJuMPScalar,
    x::AbstractVariableRef,
)
    return VectorConstraint([F, x], MOI.Complements(1))
end

function _build_complements_constraint(
    errorf::Function,
    ::AbstractJuMPScalar,
    ::AbstractJuMPScalar,
)
    return errorf("second term must be a variable.")
end

function parse_one_operator_constraint(
    errorf::Function,
    ::Bool,
    ::Union{Val{:complements},Val{:⟂}},
    F,
    x,
)
    f, parse_code = _MA.rewrite(F)
    return parse_code, :(_build_complements_constraint($errorf, $f, $(esc(x))))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Assertion

# macro assert_that(expression, messages...)
#     message = isempty(messages) ? expression : messages[1]
#     if isa(message, AbstractString)
#         message = message
#     elseif !isempty(messages) && (isa(message, Expr) || isa(message, Symbol))
#         # message is an expression needing evaluating
#         message = :(Main.Base.string($(esc(message))))
#     elseif isdefined(Main, :Base) && isdefined(Main.Base, :string) && applicable(Main.Base.string, message)
#         message = Main.Base.string(message)
#     else
#         # string() might not be defined during bootstrap
#         message = :(Main.Base.string($(Expr(:quote, message))))
#     end
#
#     parts = length(messages) < 2 ? Vector{Expr}() : collect(messages[2:end])
#     part_messages = Vector{Any}()
#     if isa(parts, Vector)
#         if !isempty(parts)
#             @info "Parts: $parts"
#             for part ∈ parts
#                 if isa(part, Expr) || isa(part, Symbol)
#                     push!(part_messages, :($part))
#                 elseif isa(part, QuoteNode)
#                     @info "Part quote: $(part.value)"
#                     push!(part_messages, :( Main.Base.string( $(Expr(:quote,part.value)) ) ))
#                 else
#                     push!(part_messages, :($part))
#                 end
#             end
#         end
#     end# Main.Base.string($(esc(part))))
#
#     return :($(esc(expression)) ? $(nothing) : throw( AssertionError($(part_messages[1])) ))
# end

# macro assert_that(expression, explanation::AbstractString)
#     return quote if $expression
#             nothing
#         else
#             formula = $(string(expression))
#             explanation_text = $(explanation)
#             throw(AssertionError("Expression: '$formula'; Explanation: $explanation_text"))
#         end
#     end
# end
#
# macro assert_that(expression, explanation::AbstractString, parts)
#     return quote
#         if $expression
#             nothing
#         else
#             formula = $(string(expression))
#             explanation_text = $(explanation)
#             parts_text = $(parts)
#             throw(AssertionError("Expression: '$formula'; Explanation: $explanation_text; Parts: $parts_text"))
#         end
#     end
# end

# function assertion_parts(parts::Vector{Pair{Expr, Any}})
#     join(map(part -> "$(part.first) = $(part.second)", parts), ", ")
# end
#
# function assertion_parts(parts::Vector{Pair{Symbol, Any}})
#     join(map(part -> "$(part.first) = $(part.second)", parts), ", ")
# end
#
# function assertion_parts(parts::Vector{Pair{Any, Any}})
#     join(map(part -> "$(part.first) = $(part.second)", parts), ", ")
# end
#
# function assertion_parts(parts::Vector)
#     join(map(
#         part -> begin
#             if("$(typeof(part).name)" == "Pair")
#                 "$(part.first) = $(part.second)"
#             else
#                 "$part"
#             end
#         end, parts), ", ")
# end
#
# export @assert_that
# export assertion_parts

end # module


text:
# Defs

export linterp
function linterp(x1, y1, x2, y2, x)
    R = typeof(zero(y1) + zero(y2))
    return R(x2 - x) ./ R(x2 - x1) .* y1 + R(x - x1) ./ R(x2 - x1) .* y2
end


text:
using OpenCL: cl
using Transpiler: CLMethod, EmptyStruct
using Sugar
import GPUArrays: _gpu_call, linear_index
using Transpiler: CLMethod
using Sugar: method_nargs, getslots!, isintrinsic, getcodeinfo!, sugared
using Sugar: returntype, type_ast, rewrite_ast, newslot!, to_tuple
using Sugar: isfunction

using Base: tail

function _gpu_call(f, A::CLArray, args::Tuple, blocks_threads::Tuple{T, T}) where T <: NTuple{N, Integer} where N
    ctx = context(A)
    _args = (KernelState(), args...) # CLArrays "state"
    clfunc = CLFunction(f, _args, ctx)
    blocks, threads = blocks_threads
    global_size = blocks .* threads
    clfunc(_args, global_size, threads)
end

device_type{T}(::T) = T
device_type{T}(::Type{T}) = Type{T}
device_type{T}(x::LocalMemory{T}) = cli.LocalPointer{T}
device_type(arg::Mem.OwnedPtr{T}) where T = cli.GlobalPointer{T}

kernel_convert{T}(x::LocalMemory{T}) = cl.LocalMem(T, x.size)
kernel_convert{T}(x::Mem.OwnedPtr{T}) = x

string_kernel_convert(x::CLArray) = pointer(x)
string_kernel_convert(x) = kernel_convert(x)

function kernel_convert(x::T) where T
    contains, fields = contains_pointer(T)
    if contains
        error("Contains pointers. Please define the correct kernel_convert/reconstruct methods")
    else
        if (isbits(x) && sizeof(x) == 0 && nfields(x) == 0) || T <: Type
            return EmptyStruct()
        end
        isbits(T) || error("Only isbits types allowed for gpu kernel")
        return x
    end
end

function match_field(field_list, field)
    for elem in field_list
        if first(elem) == field
            return [Base.tail(elem)]
        end
    end
    ()
end

makeslot(m::Sugar.LazyMethod, T, name) = newslot!(m, T, name)
makeslot(m::Tuple{}, T, name) = name

function getfield_expr(m, T, name, fname)
    # julia doesn't like symbols in getfield in its ast for rount trip
    expr = Expr(:call, getfield, name, QuoteNode(fname))
    expr.typ = T
    expr
end
function getfield_expr(m::LazyMethod, T, name, fname)
    # when we have a lazy method, we're generating code for the GPU
    expr = Expr(:call, getfield, name, fname)
    expr.typ = T
    expr
end

isa_pointer_type(T) = isa(T, Type) && (T <: OwnedPtr || T <: HostPtr || T <: cli.GlobalPointer)

function contains_pointer(
        ::Type{T},
        parent_field = (),
        pointer_fields = []
    ) where T
    @assert isleaftype(T) "Arguments must be fully concrete (leaftypes)"
    # pointers must be in parameters
    hasptr = false
    for fname in fieldnames(T)
        FT = fieldtype(T, fname)
        if isa_pointer_type(FT)
            push!(pointer_fields, (parent_field..., fname))
            hasptr = true
        else
            _hasptr, pointer_fields = contains_pointer(FT, (parent_field..., fname), pointer_fields)
            hasptr |= _hasptr
        end
    end
    hasptr, pointer_fields
end

get_fields_type(T, fields::Tuple{X}) where X = fieldtype(T, first(fields))
get_fields_type(T, fields::NTuple{N, Any}) where N = get_fields_type(fieldtype(T, first(fields)), Base.tail(fields))

get_fields(x, fields::NTuple{1}) = getfield(x, first(fields))
get_fields(x, fields::NTuple{N, Any}) where N = get_fields(getfield(x, first(fields)), Base.tail(fields))

function _getast(x::CLMethod)
    if isfunction(x)
        if isintrinsic(x)
            Expr(:block)
        else
            getcodeinfo!(x) # make sure codeinfo is present
            expr = sugared(x.signature..., code_typed)
            expr.typ = returntype(x)
            expr
        end
    else
        type_ast(x.signature)
    end
end

function funcheader_string(method, args)
    sprint() do io
        cio = Transpiler.CIO(io, method)
        Transpiler.show_returntype(cio, method)
        print(cio, ' ')
        Base.show_unquoted(cio, method)
        print(cio, '(')
        for (i, elem) in enumerate(args)
            if i != 1
                print(cio, ", ")
            end
            name, T = elem.args
            Transpiler.show_type(cio, T)
            print(cio, ' ')
            Transpiler.show_name(cio, name)
        end
        print(cio, ')')
    end
end

"""
Assembles a kernel considering that you can't pass structs containing pointers.
Will generate a kernel from function f:
```
function f(arg1, contains_ptr, arg3)
    ...
end
```
in the form of
```
// ## indicating gensymed/anonymous name
__kernel void f(arg1, ##contains_ptr, ##ptr, arg3){
    contains_ptr = reconstruct(##contains_ptr, ##ptr);
    ...
}
```
The function returns the body of the function string and a list of expressions
to set up the kernel call for clEnqueueNDRangeKernel.
"""
function assemble_kernel(m::CLMethod)
    kernel_args = []
    kernel_ptrs = []
    body = Expr(:block)
    nargs = method_nargs(m)
    # declare rest of slots
    for (i, (T, name)) in enumerate(getslots!(m)[nargs+1:end])
        slot = TypedSlot(i + nargs, T)
        push!(m.decls, slot)
        push!(m, T)
        tmp = :($name::$T)
        tmp.typ = T
        push!(body.args, tmp)
    end
    st = getslots!(m)[2:nargs] # don't include self
    arg_idx = 1
    ptr_extract = []
    for (i, (T, name)) in enumerate(st)
        argslot = TypedSlot(i + 1, T)
        contains, field_list = contains_pointer(T)
        if contains # type contains pointer to global memory
            # replace variable in arguments to free up slot
            fake_arg = gensym(name)
            T̂ = predevice_type(T) # get the pointer free equivalent type
            fake_arg_slot = newslot!(m, T̂, fake_arg)
            push!(m, T̂) # include in dependencies
            push!(kernel_args, :($fake_arg_slot::$T̂)) # put pointer free variant in kernel arguments
            push!(m.decls, fake_arg_slot)
            # we extract all pointers, pass them seperately and then assemble them in the function body
            ptr_fields = []
            ptr_slots = []
            for (ptridx, fields) in enumerate(field_list)
                ptr_arg = gensym(string("ptr_", ptridx, "_", name, arg_idx))
                PT = get_fields_type(T, fields)
                PT = cli.GlobalPointer{eltype(PT)}
                ptr_slot = newslot!(m, PT, ptr_arg)
                push!(ptr_slots, ptr_slot)
                push!(kernel_ptrs, :($ptr_slot::$PT)) # add pointer to kernel arguments
                push!(ptr_fields, (fields..., ptr_slot)) # add name of ptr to fields
                push!(ptr_extract, (i, fields...))
            end
            constr = Expr(:call, reconstruct, fake_arg_slot, ptr_slots...)
            constr.typ = T
            push!(body.args, :($argslot = $constr))
        else
            push!(m.decls, argslot)
            push!(kernel_args, :($argslot::$T))
        end
    end

    append!(kernel_args, kernel_ptrs)
    real_body = _getast(m)
    body.typ = real_body.typ # use real type
    append!(body.args, real_body.args)
    body = rewrite_ast(m, body)
    io = Transpiler.CLIO(IOBuffer(), m)
    println(io, "// Inbuilds")
    println(io, "typedef char JLBool;")
    println(io, "// dependencies")
    visited = Set()
    for dep in Sugar.dependencies!(m)
        Transpiler.print_dependencies(io, dep, visited) # this is recursive but would print method itself
    end
    println(io, "// ########################")
    println(io, "// Main inner function")
    println(io, "// ", m.signature)
    print(io, "__kernel ") # mark as kernel function
    print(io, funcheader_string(m, kernel_args))
    Base.show_unquoted(io, body, 0, 0)
    fname = string(Sugar.functionname(io, m))
    String(take!(io.io)), fname, ptr_extract
end
Transpiler.functionname(x, ::typeof(getfield)) = :getfield

const compiled_functions = Dict()

function empty_compile_cache!()
    empty!(compiled_functions)
    return
end

immutable CLFunction{F, Args, Ptrs}
    kernel::cl.Kernel
end

function CLFunction(f::F, args::T, ctx = global_context()) where {T, F}
    device = getcontext!(ctx).device
    if !supports_double(device) && any(x-> isa(x, Float64), args)
        error("Float64 is not supported by your device: $device. Make sure to convert all types for the GPU to Float32")
    end
    version = cl_version(device)
    cltypes = device_type.(args)
    get!(compiled_functions, (ctx.id, f, cltypes)) do # TODO make this faster
        method = CLMethod((f, cltypes))
        source, fname, ptr_extract = assemble_kernel(method)
        options = "-cl-denorms-are-zero -cl-mad-enable -cl-unsafe-math-optimizations"
        if version > v"1.2"
            options *= " -cl-std=CL1.2"
        end
        p = try
            cl.build!(
                cl.Program(ctx, source = source),
                options = options
            )
        catch e
            println(source)
            rethrow(e)
        end
        kernel = cl.Kernel(p, fname)
        CLFunction{F, T, Tuple{ptr_extract...}}(kernel)
    end
end

"""
Compiles a kernel from f = :kernel_name => source, with arguments args.
"""
function CLFunction(f::Pair{Symbol, String}, args::T, ctx = global_context(), options = "-cl-denorms-are-zero -cl-mad-enable -cl-unsafe-math-optimizations") where T
    device = getcontext!(ctx).device
    if !supports_double(device) && any(x-> isa(x, Float64), args)
        error("Float64 is not supported by your device: $device. Make sure to convert all types for the GPU to Float32")
    end
    version = cl_version(device)
    fname, source = f
    get!(compiled_functions, (ctx.id, fname, T)) do # TODO make this faster
        if version > v"1.2"
            options *= " -cl-std=CL1.2"
        end
        p = try
            cl.build!(
                cl.Program(ctx, source = source),
                options = options
            )
        catch e
            println(source)
            rethrow(e)
        end
        kernel = cl.Kernel(p, string(fname))
        # drop kernelstate from args
        CLFunction{String, T, Tuple{}}(kernel)
    end
end

function (clf::CLFunction{F, Args, Ptrs})(
        args::Args, blocks::NTuple{N, Integer}, threads, queue = global_queue()
    ) where {F, Args, Ptrs, N}
    @assert N in (1, 2, 3) "invalid block size"
    gsize = Csize_t[b for b in blocks]
    lsize = isa(threads, Tuple) ? Csize_t[threads[i] for i=1:length(blocks)] : threads
    kernel = clf.kernel
    idx = 1
    for elem in args
        tmp = kernel_convert(elem)
        cl.set_arg!(kernel, idx, tmp)
        idx += 1
    end
    for fields in to_tuple(Ptrs)
        arg_idx, fields = first(fields), tail(fields)
        cl.set_arg!(kernel, idx, kernel_convert(get_fields(args[arg_idx], fields)))
        idx += 1
    end
    ret_event = Ref{cl.CL_event}()
    cl.@check cl.api.clEnqueueNDRangeKernel(
        queue.id, kernel.id,
        length(gsize), C_NULL, gsize, lsize,
        0, C_NULL, ret_event
    )
    return cl.Event(ret_event[], retain = false)
end

# Kernels compiled from a string have a different conversion semantic.
function (clf::CLFunction{String, Args, Tuple{}})(
        args::Args, blocks::NTuple{N, Integer}, threads, queue = global_queue()
    ) where {Args, N}
    @assert N in (1, 2, 3) "invalid block size"
    gsize = Csize_t[b for b in blocks]
    lsize = isa(threads, Tuple) ? Csize_t[threads[i] for i=1:length(blocks)] : threads
    kernel = clf.kernel
    idx = 1
    for elem in Iterators.drop(args, 1) # drop kernel state
        tmp = string_kernel_convert(elem)
        cl.set_arg!(kernel, idx, tmp)
        idx += 1
    end
    ret_event = Ref{cl.CL_event}()
    cl.@check cl.api.clEnqueueNDRangeKernel(
        queue.id, kernel.id,
        length(gsize), C_NULL, gsize, lsize,
        0, C_NULL, ret_event
    )
    return cl.Event(ret_event[], retain = false)
end


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using InteractiveDynamics, Random
using DynamicalBilliards
import GLMakie
using GLMakie: to_color, RGBf0, RGBAf0
using DynamicalSystems

# Set style to book colors
InteractiveDynamics.obcolor(::Antidot) = to_color(COLORS[3])
InteractiveDynamics.obcolor(::Obstacle) = to_color(COLORS[3])
InteractiveDynamics.obfill(::Antidot) = RGBAf0(0,0,0,0)
InteractiveDynamics.obls(::Antidot) = nothing

w = 0.2f0
bd = billiard_mushroom(1f0, w, 1f0, 0f0)

pc = MushroomTools.randomchaotic(1f0, w, 1f0)
pc = Particle(-0.01f0, 0.2f0, sqrt(3f0)/2)

pr = Particle(0.0f0, 1.2f0, 0.0f0)
pr2 = Particle(0.0f0, 1.9f0, 0.0f0)

ps = [pc, pr, pr2]
colors = [to_color(COLORS[i]) for i in (1,2,4)]

fig, bmapax = billiard_bmap_plot(bd, ps; res = (1600, 500),
colors = colors, tail = 100000, steps = 1000003, backgroundcolor = RGBf0(1,1,1),
arrowsize = 15, particle_size = 1.5,
ms = 5)

display(fig)

# interactive_billiard_bmap(bd)

# bmapax.xticklabelsize = 32
# bmapax.yticklabelsize = 32
bmapax.xlabelsize = 40
bmapax.ylabelsize = 40
bmapax.ylabelpadding = 15
# GLMakie.save(plotsdir("8", "mushroom.png"), fig)
GLMakie.save(plotsdir("8", "mushroom2.png"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using QuadGK
using LambertW
using Base.Threads

"""
Compute range y1, y2 of the Converse Jensen Inequality:
from 0 to y1: ρ1 sum to x1;
from 0 to y2: ρ2 sum to x2.
"""
function converse_jensen_range(x1::Float64, x2::Float64, xmax::Float64=1-1e-10)::Tuple{Float64,Float64}
    # Cap x1, x2 by xmax because LambertW package may have large error near 1/ℯ
    x1 = min(x1, xmax)
    x2 = min(x2, xmax)
    if 1-(1-x1)*(1-log(1-x1)) >= x2
        # case 1: λ2 <= λ1.
        # Solving 1-(1+λ)exp(-λ)=x2 gives λ2=-1-W(-(1-x2)/ℯ) (-1 branch)
        λ2 = -1-lambertw(-(1-x2)/ℯ, -1)
    else
        # case 1: λ2 > λ1
        # Solving 2-(2+λ)exp(-λ)=x1+x2 gives λ2=-2-W(-(2-x1-x2)/ℯ^2) (-1 branch)
        λ2 = -2-lambertw(-(2-x1-x2)/ℯ^2, -1)
    end
    y1 = x1
    y2 = min(1-exp(-λ2), xmax)
    return y1, y2
end

"""Compute the contribution at any y = 1-exp(-λ) to (d/dt) log(d_(x1, x2)(t))."""
function dlog_double_dy(y::Float64, y1::Float64, y2::Float64, x1::Float64, x2::Float64, t::Float64, logd::Float64)::Float64
    ρ1 = (y<=y1) ? 1.            : 0.
    ρ2 = (y<=y2) ? 1-log(1-y)-ρ1 : 0.
    w1 = exp((2*x1+x2)*t+logd)
    w2 = exp((x1+2*x2)*t+logd)
    return (w1*ρ1+w2*ρ2)/((max(w1-1,0.)*ρ1+max(w2-1,0.)*ρ2)*(1-y)+1)
end

"""Compute the value of log(q_x(t)) for all 0 <= t <=1 that is a multiple of dt, for a given x."""
function logq(x::Float64, dx::Float64=1e-3, dt::Float64=1e-3, xmax::Float64=1-1e-10, rtol::Float64=1e-5)::Array{Float64,1}

    x_grid = Array(0.:dx:1-dx)
    nx = length(x_grid)

    t_grid = Array(0.:dt:1-dt)
    nt = length(t_grid)

    logq_grid = zeros(nx,nt)

    @threads for i in 1:nx
        x1, x2 = (x>=x_grid[i]) ? [x, x_grid[i]] : [x_grid[i], x]
        y1, y2 = converse_jensen_range(x1, x2, xmax)
        logd = 0.
        for j in 1:(nt-1)
            t = t_grid[j]
            dlogd, err = quadgk(y -> dlog_double_dy(y, y1, y2, x1, x2, t, logd), 0., xmax, rtol=rtol)
            logd2      = logd-(dlogd-err)*dt
            dlogd, err = quadgk(y -> dlog_double_dy(y, y1, y2, x1, x2, t, logd2), 0., xmax, rtol=rtol)
            logd      -= (dlogd-err)*dt
            logd       = min(logd, -(x1+x2)*(t+dt))   # use trivial bound if numerical bound is worse (for small t)
            logq_grid[i,j+1] = min(logd+(x_grid[i]+dx)*(t+dt), -x*(t+dt))
        end
    end

    return vec(maximum(logq_grid, dims=1))
end

"""Compute match probability bound of an vertex with LP matched level x"""
function level2(x::Float64=1., dx::Float64=1e-3, dt::Float64=1e-3, xmax::Float64=1-1e-10, rtol::Float64=1e-5)::Float64
    log_single = 0.
    t_grid = Array(0.:dt:1-dt)
    nt = length(t_grid)
    logq_grid = logq(x, dx, dt, xmax, rtol)
    for i in 1:nt
        t = t_grid[i]
        y = 1-exp(-x*t+logq_grid[i]-log_single)
        dlog_single = (y==0.) ? x : -log(1-x*y)/y
        log_single2 = log_single - dlog_single*dt
        y = 1-exp(-x*t+logq_grid[i]-log_single2)
        dlog_single = (y==0.) ? x : -log(1-x*y)/y
        log_single -= dlog_single*dt
    end
    return 1-exp(log_single)
end

max_stars_repo_name: h-zhiyi/poisson-ocs

max_issues_repo_name: h-zhiyi/poisson-ocs

max_forks_repo_name: h-zhiyi/poisson-ocs


text:
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("colorscheme.jl"))

using InteractiveDynamics
using DynamicalSystems, OrdinaryDiffEq
import GLMakie

a = 0.2
b = 0.2
c = 5.7
amu = 0.02
ω1 = 1.1 + amu
ω2 = 1.1 - amu
k1 = 0.02
k2 = k1

y_init = [0.1, 0.2, 0., 0.11, 0.19, 0.1]

diffeq = (alg = Tsit5(), adaptive = false, dt = 0.02, reltol = 1e-9, abstol = 1e-9)

ds = Systems.coupled_roessler(y_init; ω1, ω2, a, b, c, k1, k2)

ps = Dict(
    1 => 1.0:0.01:1.5,
    2 => 1.0:0.01:1.5,
    5 => 0:0.01:6,
    6 => 0:0.001:0.2,
    7 => 0:0.001:0.2,
)
pnames = Dict(
    1 => "ωx",
    2 => "ωy",
    5 => "c",
    6 => "kx",
    7 => "ky",
)

u0s =  [ds.u0]

fig, obs = interactive_evolution_timeseries(
    ds, u0s, ps; tail = 1000, diffeq, idxs = [1,4], pnames,
    lims = ((-12,12), (-12,12))
)

ax = GLMakie.content(fig[1,1])
ax.xlabel = "x1"
ax.ylabel = "y1"

# %%
record_interaction(fig, string(@__FILE__)[1:end-2]*"mp4"; total_time = 30, sleep_time = 2)


text:
module ADSeismic

    using ADCME
    using Optim
    using PyPlot
    using PyCall
    using SparseArrays
    using LinearAlgebra
    using Statistics
    using Parameters
    using DelimitedFiles
    using MAT
    using Printf

    np = PyNULL()
    HAS_GPU = has_gpu()
    HAS_MPI = has_mpi(false)
    function __init__()
        copy!(np, pyimport("numpy"))
    end

    include("Struct.jl")
    include("Core.jl")
    include("IO.jl")
    include("NN.jl")
    include("Optim.jl")

    include("ThirdParty.jl")
    include("MPIAcoustic.jl")
    include("MPIElastic.jl")
    include("Utils.jl")

end


text:
##
cd(@__DIR__())
rx = r"([a-z][a-z][a-z]) ((\+|\-)([0-9]+))"
code = map(eachline("input.txt")) do line
    m = match(rx, line)
    String(m[1]), parse(Int,m[2])
end

mutable struct Handheld
    count::Int
    ip::Int
    running::Bool
end

Handheld() = Handheld(0, 1, false)
function start!(hh::Handheld)
    hh.count = 0
    hh.ip = 1
    hh.running = true
end
isrunning(hh::Handheld) = hh.running

function run_code!(hh::Handheld, code, jp::Int)
    start!(hh)
    visited = Set{Int}()
    inc = 1
    exited = false
    c, i = code[jp]
    if c == "jmp"
        code[jp] = ("nop", i)
    elseif c == "nop"
        code[jp] = ("jmp", i)
    end
    while isrunning(hh)
        if hh.ip < 1 || hh.ip > length(code)
            @error "ip out of bounds"
        end
        c, i = code[hh.ip]
        if hh.ip ∈ visited
            @info "Second run: $(hh.count)" hh
            break
        end
        push!(visited, hh.ip)
        if c == "acc"
            hh.count += i
            @debug "acc $i; ip = $(hh.ip); count = $(hh.count)"
        elseif c == "jmp"
            inc = i
            @debug "jmp $i; ip = $(hh.ip); count = $(hh.count)"
        elseif c == "nop"
            @debug "NOP $i; ip = $(hh.ip); count = $(hh.count)"
        else
            @error "Unknown instruction: $c; ip = $(hh.ip)"
            break
        end
        hh.ip += inc
        if hh.ip > length(code)
            hh.running = false
            @debug "should stop" hh
            if inc == 1
                exited = true
            end
        end
        inc = 1
    end
    @info "count is: $(hh.count)"
    return exited
end

##
ENV["JULIA_DEBUG"] = nothing

##
function find_error(hh::Handheld, code)
    jmps = findall(x->x[1]=="jmp", code)
    nops = findall(x->x[1]=="nop", code)
    e = false
    for j in jmps
        e = run_code!(hh, copy(code), j)
        if e
            @info "Found error; count = $(hh.count)" hh
            break
        end
    end
    e && return
    for n in nops
        e = run_code!(hh, copy(code), n)
        if e
            @info "Found error; count = $(hh.count)" hh
            break
        end
    end
end

##
hh = Handheld()
find_error(hh, code)

## PARTE 2


text:
immutable PingFrame
    is_ack::Bool
    data::Array{UInt8, 1}
end

==(a::PingFrame, b::PingFrame) =
    a.is_ack == b.is_ack &&
    a.data == b.data

function decode_ping(header, payload)
    @assert length(payload) == 8
    @assert header.stream_identifier == 0x0

    is_ack = header.flags & 0x1 == 0x1
    return PingFrame(is_ack, payload)
end

function encode_ping(frame)
    flags = 0x0 | (frame.is_ack ? 0x1 : 0x0)

    return wrap_payload(frame.data, PING, flags, 0x0)
end


text:
"""
Convert a polynomial to a string for printing.
"""
function poly_to_string(poly::Polynomial)
    io = IOBuffer();
    printpoly(io, poly, descending_powers=true)
    return String(take!(io))
end

function show_time_delay(io::IO, time_delay::Union{Int, Float64})
    if time_delay != 0.0
        print(io, " e^($(-time_delay)*s)")
    end
end

function Base.show(io::IO, tf::TransferFunction)
    top = poly_to_string(tf.numerator)
    bottom = poly_to_string(tf.denominator)

    # is this a rational polynomial? then just print the numerator / constant
    if degree(tf.denominator) == 0
        # normalized numerator
        p = tf.numerator / tf.denominator[0]
        print(io, poly_to_string(p))
        show_time_delay(io, tf.time_delay)
        return nothing
    end

    nb_dashes = maximum([length(top), length(bottom)])

    # print numerator
    nb_slack = floor(Int, (nb_dashes - length(top)) / 2) # for centering
    for i = 1:nb_slack
        print(io, " ")
    end
    println(io, top)

    # print dashes to separate numerator and denominator
    for _ = 1:nb_dashes
        print(io, "-")
    end
    show_time_delay(io, tf.time_delay)
    print(io, "\n")

    # print denominator
    nb_slack = floor(Int, (nb_dashes - length(bottom)) / 2) # for centering
    for i = 1:nb_slack
        print(io, " ")
    end
    print(io, bottom)
    return nothing
end

function Base.show(io::IO, g::ClosedLoopTransferFunction)
    println(io, "closed-loop transfer function.
      top
    -------
    1 + g_ol")
    println(io, "\n  top = ")
    show(io, g.top)
    println(io, "\n\n  g_ol = ")
    show(io, g.g_ol)
end

function Base.show(io::IO, pc::PController)
    println(io, "P controller")
    println(io, "\tcontroller gain Kc = ", pc.Kc)
end

function Base.show(io::IO, pic::PIController)
    println(io, "PI controller")
    println(io, "\tcontroller gain Kc = ", pic.Kc)
    println(io, "\tintegral time constant τI = ", pic.τI)
end

function Base.show(io::IO, pidc::PIDController)
    println(io, "PID controller")
    println(io, "\tcontroller gain Kc = ", pidc.Kc)
    println(io, "\tintegral time constant τI = ", pidc.τI)
    println(io, "\tderivative time constant τD = ", pidc.τD)
    println(io, "\tderivative filter α = ", pidc.α)
end

function Base.show(io::IO, m::Margins)
    println(io, "-- gain/phase margin info--")
    @printf(io, "\tcritical frequency ω_c [rad/time]:       %.5f\n", m.ω_c)
    @printf(io, "\tgain crossover frequency ω_g [rad/time]: %.5f\n", m.ω_g)
    @printf(io, "\tgain margin:                             %.5f\n", m.gain_margin)
    @printf(io, "\tphase margin:                            %.5f\n", m.phase_margin)
end

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
using OpenQuantumBase, LinearAlgebra, Test

num_qubits = 3
H₃ = random_ising(num_qubits)
@test size(H₃) == (2^num_qubits, 2^num_qubits)
@test isdiag(H₃)
@test !issparse(H₃)
@test issparse(random_ising(num_qubits, sp=true))

@test alt_sec_chain(1,0.5,1,3) == σz⊗σz⊗σi + 0.5*σi⊗σz⊗σz


text:
mutable struct VIPData
    F::Array{JuMP.NonlinearExpression,1}
    var::Array{VariableRef,1}
    relation::Dict{VariableRef, JuMP.NonlinearExpression}
end

function VIPModel(optimizer=Ipopt.Optimizer)
    m = JuMP.Model(with_optimizer(optimizer, max_cpu_time=60.0, print_level=1))     ### print_level=0 ?

    m.ext[:VIP] = VIPData(Array{JuMP.NonlinearExpression}[], Array{VariableRef}[], Dict{VariableRef, JuMP.NonlinearExpression}() )
    return m
end

function getVIPData(m::Model)
    if haskey(m.ext, :VIP)
        return m.ext[:VIP]::VIPData
    else
        error("The 'getVIPData' function is only for VIP models as in VariationalInequality.jl")
    end
end

@eval const $(Symbol("@mapping")) = $(Symbol("@NLexpression"))

macro innerproduct(args...)
  if length(args) != 3
    error("3 arguments are required in @innerproduct(...)")
  end

  m = esc(args[1])
  F_name = esc(args[2])
  var = esc(args[3])

  code = quote
    innerproduct($m, $F_name, $var)
  end

  return code
end

################################################################################
# correspond interface
# The most basic/important one. All other interfaces call this method.
function innerproduct(m::Model, expression::JuMP.NonlinearExpression, variable::VariableRef)
    data = getVIPData(m)
    data.relation[variable] = expression
end

# Alternative
function innerproduct(m::Model, variable::VariableRef, expression::JuMP.NonlinearExpression)
    innerproduct(m, expression, variable)
end

#=
# Do we need the below?
function correspond(m::Model, variables::Array{VariableRef,1}, expressions::Array{JuMP.NonlinearExpression,1})
    correspond(m, expressions, variables)
end
function correspond(m::Model, expressions::Array{JuMP.NonlinearExpression,1}, variables::Array{VariableRef,1})
    @assert length(expressions) == length(variables)
    for i in 1:length(variables)
        correspond(m, expressions[i], variables[i])
    end
end
=#

function innerproduct(m::Model, expressions::Array{JuMP.NonlinearExpression}, variables::Array{VariableRef})
    expressions = collect(expressions)
    variables = collect(variables)
    @assert length(expressions) == length(variables)
    for i in 1:length(variables)
        innerproduct(m, expressions[i], variables[i])
    end
end

function innerproduct(m::Model, variables::Array{VariableRef}, expressions::Array{JuMP.NonlinearExpression})
    innerproduct(m, expressions, variables)
end

#=
# Do we need the below?
function correspond(m::Model, expressions::JuMP.JuMPArray, variables::JuMP.JuMPArray)
    variables = collect(variables.innerArray)
    expressions = collect(expressions.innerArray)
    @assert length(expressions) == length(variables)
    for i in 1:length(variables)
        correspond(m, expressions[i], variables[i])
    end
end
=#
################################################################################

#=
function setVIP(m::Model, expressions::Array{JuMP.NonlinearExpression,1}, variables::Array{VariableRef,1})
    @assert length(expressions) == length(variables)

    data = getVIPData(m)
    data.F = expressions
    data.var = variables

    return m
end
=#

#=
function solveVIP(m::Model    ; step_size=0.01,
                                algorithm=:fixed_point,
                                tolerance=1e-10,
                                max_iter=1000            )

    F = getVIPData(m).F
    var = getVIPData(m).var
    @assert length(F) == length(var)
    dim = length(F)

    x0 = zeros(dim)

    return solveVIP(m, x0, step_size=step_size, algorithm=algorithm, tolerance=tolerance, max_iter=max_iter)

end
=#

# Copied from JuMP.jl
# internal method that doesn't print a warning if the value is NaN
_get_value(v::VariableRef) = (has_values(v.model) ? value(v) : NaN)

function clear_values(m)
    relation = getVIPData(m).relation
    for variable in keys(relation)
        setvalue(variable, NaN)
    end
end

function initial_projection(m::JuMP.Model)
    relation = getVIPData(m).relation

    initial_values = Dict{VariableRef, Float64}()
    for variable in keys(relation)
        val = _get_value(variable)
        if isnan(val)
            val = 0.0
        end
        initial_values[variable] = val
    end

    @objective(m, Min, sum( ( variable - initial_values[variable] )^2 for variable in keys(relation) ) )

    optimize!(m)
    status = termination_status(m)

    @show status
    @show m

    @assert status == MOI.LOCALLY_SOLVED || status == MOI.OPTIMAL
end

function gap_function(m)
    relation = getVIPData(m).relation
    var = get_variables(relation)
    y = get_current_x(relation)
    Fy = get_current_F(relation)

    @objective(m, Max, sum( Fy[j] * ( y[j] - var[j] ) for j in 1:length(var)) )
    optimize!(m)
    return getobjectivevalue(m)
end

function save_solution(m)
    relation = getVIPData(m).relation
    var = get_variables(relation)
    x = get_current_x(relation)
    F = get_current_F(relation)

    solution = Dict(zip(var, x))
    F_value = Dict(zip(var, F))
    gap = gap_function(m)

    return solution, F_value, gap
end

function solveVIP(m::Model; step_size=0.01, algorithm=:fixed_point, tolerance=1e-6, max_iter=1000)
    initial_projection(m)

    if algorithm==:fixed_point
        _fixed_point(m, step_size, tolerance, max_iter)
    elseif algorithm==:extra_gradient
        _extra_gradient(m, step_size, tolerance, max_iter)
    # elseif algorithm==:hyperplane
    #     _hyperplane(m, step_size, tolerance, max_iter)
    end

    # sol, Fval, gap = save_solution(m)
    return save_solution(m)
end


text:
# methods/inplace.jl
#
# Method-based access to tensor operations using inplace definitions.

tensorcopy!(A, IA, C, IC) = tensoradd!(1, A, IA, 0, C, IC)

function tensoradd!(α, A, IA, β, C, IC)
    checkindices(A, IA)
    checkindices(C, IC)
    indCinA = add_indices(IA, IC)
    add!(α, A, Val{:N}, β, C, indCinA)
end

function tensortrace!(α, A, IA, β, C, IC)
    checkindices(A, IA)
    checkindices(C, IC)
    indCinA, cindA1, cindA2 = trace_indices(IA, IC)
    trace!(α, A, Val{:N}, β, C, indCinA, cindA1, cindA2)
    return C
end

function tensorcontract!(α, A, IA, conjA, B, IB, conjB, β, C, IC; method::Symbol = :BLAS)
    # Updates C as β*C+α*contract(A, B), whereby the contraction pattern
    # is specified by IA, IB and IC. The iterables IA(B, C)
    # should contain a unique label for every index of array A(B, C), such that
    # common I of A and B correspond to indices that will be contracted.
    # Common I between A and C or B and C indicate the position of the
    # uncontracted indices of A and B with respect to the indices of C, such
    # that the output array of the contraction can be added to C. Every label
    # should thus appear exactly twice in the union of IA, IB and
    # IC and the associated indices of the tensors should have identical
    # size.
    # Array A and/or B can be also conjugated by setting conjA and/or conjB
    # equal  to 'C' instead of 'N'.
    # The parametric argument method can be specified to choose between two
    # different contraction strategies:
    # -> method = :BLAS : permutes tensors (requires extra memory) and then
    #                   calls built-in (typically BLAS) multiplication
    # -> method = :native : memory-free native julia tensor contraction

    checkindices(A, IA)
    checkindices(B, IB)
    checkindices(C, IC)

    conjA == 'N' || conjA == 'C' || throw(ArgumentError("Value of conjA should be 'N' or 'C' instead of $conjA"))
    conjB == 'N' || conjB == 'C' || throw(ArgumentError("Value of conjB should be 'N' or 'C' instead of $conjB"))
    CA = conjA == 'N' ? :N : :C
    CB = conjB == 'N' ? :N : :C

    oindA, cindA, oindB, cindB, indCinoAB = contract_indices(IA, IB, IC)

    if method == :BLAS
        contract!(α, A, Val{CA}, B, Val{CB}, β, C, oindA, cindA, oindB, cindB, indCinoAB,Val{:BLAS})
    elseif method == :native
        contract!(α, A, Val{CA}, B, Val{CB}, β, C, oindA, cindA, oindB, cindB, indCinoAB,Val{:native})
    else
        throw(ArgumentError("unknown contraction method"))
    end
    return C
end

function tensorproduct!(α, A, IA, B, IB, β, C, IC)
    isempty(intersect(IA, IB)) || throw(LabelError("not a valid tensor product"))
    tensorcontract!(α, A, IA, 'N', B, IB, 'N', β, C, IC; method = :native)
end


text:
using MCObservables

function test_VectorObservable(MCS, SObs, VObs)
  inform("Start testing $VObs")

  n = 2

  s = zeros(SObs, n)
  v = VObs()

  for i in 1:MCS
    x = zeros(n)
    for j in 1:n
      x[j] = randn()
      s[j] << x[j]
    end
    v << x
  end

  vjk = jackknife(v)

  for i in 1:n
    @test mean(s[i]) ==  mean(v)[i]
    @test stderror(s[i]) == stderror(v)[i]
    sjk = jackknife(s[i])
    @test_approx_eq_eps mean(sjk) mean(vjk)[i] 1e-10
    @test_approx_eq_eps stderror(sjk) stderror(vjk)[i] 1e-10
  end

  inform("Success testing $VObs")
end

result_simpleObservable = test_VectorObservable(8192, SimpleObservable, SimpleVectorObservable)
result_simpleObservable = test_VectorObservable(8192, BinningObservable, BinningVectorObservable)


text:
# This file is a part of GeDetROOT.jl, licensed under the MIT License (MIT).

import Compat.Test
Test.@testset "Package GeDetROOT" begin
end


text:
function CategoricalValue{T, R}(level::Integer, pool::CategoricalPool{T, R})
    return CategoricalValue(convert(R, level), pool)
end

Base.convert{T, R}(::Type{CategoricalValue{T, R}}, x::CategoricalValue{T, R}) = x
Base.convert{T}(::Type{CategoricalValue{T}}, x::CategoricalValue{T}) = x
Base.convert(::Type{CategoricalValue}, x::CategoricalValue) = x

Base.promote_rule{S, T, R}(::Type{CategoricalValue{S, R}}, ::Type{T}) = promote_type(S, T)
Base.promote_rule{S, T}(::Type{CategoricalValue{S}}, ::Type{T}) = promote_type(S, T)
Base.promote_rule{T}(::Type{CategoricalValue}, ::Type{T}) = T

# To fix ambiguities with definitions from Base
Base.convert{S}(::Type{Nullable{S}}, x::CategoricalValue{Nullable}) =
    convert(Nullable{S}, index(x.pool)[x.level])
Base.convert{S}(::Type{Nullable}, x::CategoricalValue{S}) = convert(Nullable{S}, x)
Base.convert{T}(::Type{Nullable{CategoricalValue{Nullable{T}}}},
                x::CategoricalValue{Nullable{T}}) =
    Nullable(x)
Base.convert{T}(::Type{Ref}, x::CategoricalValue{T}) = RefValue{T}(x)

Base.convert{S, T}(::Type{S}, x::CategoricalValue{T}) = convert(S, index(x.pool)[x.level])

function Base.show{T}(io::IO, x::CategoricalValue{T})
    if @compat(get(io, :compact, false))
        print(io, repr(index(x.pool)[x.level]))
    elseif ordered(x.pool)
        @printf(io, "%s %s (%i/%i)",
                typeof(x),
                repr(index(x.pool)[x.level]),
                order(x.pool)[x.level], length(x.pool))
    else
        @printf(io, "%s %s",
                typeof(x),
                repr(index(x.pool)[x.level]))
    end
end

if VERSION < v"0.5.0-dev+1936"
    Base.showcompact{T}(io::IO, x::CategoricalValue{T}) =
        print(io, repr(index(x.pool)[x.level]))
end

@compat Base.:(==)(x::CategoricalValue, y::CategoricalValue) =
    index(x.pool)[x.level] == index(y.pool)[y.level]

# To fix ambiguities with Base
@compat Base.:(==)(x::CategoricalValue, y::WeakRef) = index(x.pool)[x.level] == y
@compat Base.:(==)(x::WeakRef, y::CategoricalValue) = y == x

@compat Base.:(==)(x::CategoricalValue, y::Any) = index(x.pool)[x.level] == y
@compat Base.:(==)(x::Any, y::CategoricalValue) = y == x

Base.isequal(x::CategoricalValue, y::CategoricalValue) =
    isequal(index(x.pool)[x.level], index(y.pool)[y.level])
Base.isequal(x::CategoricalValue, y::Any) = isequal(index(x.pool)[x.level], y)
Base.isequal(x::Any, y::CategoricalValue) = isequal(y, x)

Base.hash(x::CategoricalValue, h::UInt) = hash(index(x.pool)[x.level], h)

function Base.isless{S, T}(x::CategoricalValue{S}, y::CategoricalValue{T})
    error("CategoricalValue objects with different pools cannot be tested for order")
end

function Base.isless{T}(x::CategoricalValue{T}, y::CategoricalValue{T})
    if x.pool !== y.pool
        error("CategoricalValue objects with different pools cannot be tested for order")
    elseif !ordered(x.pool) # !ordered(y.pool) is implied by x.pool === y.pool
        error("Unordered CategoricalValue objects cannot be tested for order; use the ordered! function on the parent array to change this")
    else
        return isless(order(x.pool)[x.level], order(y.pool)[y.level])
    end
end

Base.get{T,R}(x::CategoricalValue{T,R}) = convert(T,x)


text:
"""Reference for parameters update when present"""
struct UpdateRef{T}
    access_ref::Symbol
    accessor_func::Union{Nothing, String}
end

function UpdateRef{T}(name::Symbol) where {T <: Union{JuMP.VariableRef, PJ.ParameterRef}}
    return UpdateRef{T}(name, nothing)
end

function UpdateRef{T}(
    name::AbstractString,
    accessor_func::Union{Nothing, String} = nothing,
) where {T <: Union{JuMP.VariableRef, PJ.ParameterRef}}
    return UpdateRef{T}(Symbol(name), accessor_func)
end

function UpdateRef{T}(
    ::Type{U},
    name::AbstractString,
) where {T <: Union{JuMP.VariableRef, PJ.ParameterRef}, U <: PSY.Component}
    return UpdateRef{T}(encode_symbol(U, name), nothing)
end

function UpdateRef{T}(
    name::AbstractString,
    accessor_func::AbstractString,
) where {T <: PSY.Component}
    # Combine these three fields together in order to guarantee uniqueness.
    return UpdateRef{T}(encode_symbol(T, name, accessor_func), accessor_func)
end

function get_accessor_func(ref::UpdateRef{T}) where {T <: PSY.Component}
    if isnothing(ref.accessor_func)
        throw(IS.InvalidValue("accessor_func is not defined for $ref"))
    end

    return ref.accessor_func
end


text:
module IonosphereWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Class = @jimport java.lang.Class
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const DateComponents = @jimport org.orekit.time.DateComponents
const EstimatedIonosphericModel = @jimport org.orekit.models.earth.ionosphere.EstimatedIonosphericModel
const Field = @jimport org.hipparchus.Field
const FieldAbsoluteDate = @jimport org.orekit.time.FieldAbsoluteDate
const FieldGeodeticPoint = @jimport org.orekit.bodies.FieldGeodeticPoint
const FieldNeQuickParameters = @jimport org.orekit.models.earth.ionosphere.FieldNeQuickParameters
const FieldSpacecraftState = @jimport org.orekit.propagation.FieldSpacecraftState
const GeodeticPoint = @jimport org.orekit.bodies.GeodeticPoint
const GlobalIonosphereMapModel = @jimport org.orekit.models.earth.ionosphere.GlobalIonosphereMapModel
const InputStream = @jimport java.io.InputStream
const IonosphericMappingFunction = @jimport org.orekit.models.earth.ionosphere.IonosphericMappingFunction
const IonosphericModel = @jimport org.orekit.models.earth.ionosphere.IonosphericModel
const JString = @jimport java.lang.String
const KlobucharIonoCoefficientsLoader = @jimport org.orekit.models.earth.ionosphere.KlobucharIonoCoefficientsLoader
const KlobucharIonoModel = @jimport org.orekit.models.earth.ionosphere.KlobucharIonoModel
const List = @jimport java.util.List
const NeQuickModel = @jimport org.orekit.models.earth.ionosphere.NeQuickModel
const NeQuickParameters = @jimport org.orekit.models.earth.ionosphere.NeQuickParameters
const Object = @jimport java.lang.Object
const RealFieldElement = @jimport org.hipparchus.RealFieldElement
const SingleLayerModelMappingFunction = @jimport org.orekit.models.earth.ionosphere.SingleLayerModelMappingFunction
const SpacecraftState = @jimport org.orekit.propagation.SpacecraftState
const TimeScale = @jimport org.orekit.time.TimeScale
const TopocentricFrame = @jimport org.orekit.frames.TopocentricFrame

export EstimatedIonosphericModel
export FieldNeQuickParameters
export GlobalIonosphereMapModel
export IonosphericMappingFunction
export IonosphericModel
export KlobucharIonoCoefficientsLoader
export KlobucharIonoModel
export NeQuickModel
export NeQuickParameters
export SingleLayerModelMappingFunction
export equals
export get_alpha
export get_b1_bot
export get_b1_top
export get_b2_bot
export get_be_bot
export get_be_top
export get_beta
export get_class
export get_h0
export get_hm_e
export get_hm_f1
export get_hm_f2
export get_layer_amplitudes
export get_nm_f2
export get_parameters
export get_parameters_drivers
export get_supported_names
export get_tec
export hash_code
export load_data
export load_klobuchar_ionospheric_coefficients
export mapping_factor
export notify
export notify_all
export path_delay
export stec
export still_accepts_data
export to_string
export wait

include("estimated_ionospheric_model.jl")
include("field_ne_quick_parameters.jl")
include("global_ionosphere_map_model.jl")
include("ionospheric_mapping_function.jl")
include("ionospheric_model.jl")
include("klobuchar_iono_coefficients_loader.jl")
include("klobuchar_iono_model.jl")
include("ne_quick_model.jl")
include("ne_quick_parameters.jl")
include("single_layer_model_mapping_function.jl")

end


text:
######################################################
#               ALGORITMO FORWARD PROPAGATION
# Algoritmo 3 descrito en la memoria. Capítulo 5.
######################################################
export forward_propagation
"""
forward_propagation (h::AbstractOneLayerNeuralNetwork, activation_function, x::Vector{Real})
Only use an activation function
"""
function forward_propagation(h::AbstractOneLayerNeuralNetwork,activation_function, x)
    s = h.W1 * push!(copy(x),1)
    ∑= map(activation_function,s)
    return h.W2 * ∑
end


text:
"""
FunctionBlocks
"""

export StepBlock

mutable struct StepBlock <: AbstractFunctionBlock
    inport::Vector{InPort}
    outport::Vector{OutPort}
    steptime::AbstractSignalLine
    initialvalue::AbstractSignalLine
    finalvalue::AbstractSignalLine

    function StepBlock(::Type{Tv} = Float64; steptime = 0, initialvalue = 0, finalvalue = 0) where Tv
        @createblock new(Vector{InPort}(), Vector{OutPort}(),
            signal(steptime, Tv),
            signal(initialvalue, Tv),
            signal(finalvalue, Tv)
            ) 0 1
    end
end

"""
IO
"""

function Base.show(io::IO, blk::StepBlock)
    println(io, "Step()")
end

"""
to expr
"""

function _toexpr(blk::StepBlock)
    steptime = _toexpr(blk.steptime)
    initialvalue = _toexpr(blk.initialvalue)
    finalvalue = _toexpr(blk.finalvalue)
    :((t < $steptime) ? $initialvalue : $finalvalue)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using BinaryBuilderBase

name = "finufft"
version = v"2.0.4"
julia_compat = "1.6"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/flatironinstitute/finufft/archive/v$(version).zip", "2434f694b4fbdbeb65c77f65d784a1712852130b9c61e15999555a2e2cf1a9fa")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/finufft*/
# Overwrite LIBSFFT such that we do not require fftw3_threads or fftw3_omp for OMP support. Since the libraries in FFTW_jll already provide for threading, we do not loose anything.
# Make use of the -DFFTW_PLAN_SAFE flag to allow for multiple threads using finufft at the same time.
make lib CFLAGS="-fopenmp -fPIC -O3 -funroll-loops -fcx-limited-range -Iinclude" CXXFLAGS="-fopenmp -fPIC -O3 -funroll-loops -fcx-limited-range -Iinclude -std=c++14 -DFFTW_PLAN_SAFE" LIBSFFT="-lfftw3 -lfftw3f -lm"
mv lib/libfinufft.so "${libdir}/libfinufft.${dlext}"
"""

# Build for all supported platforms
platforms_x86_64 = filter(p -> p.tags["arch"]=="x86_64", supported_platforms())
platforms_other  = filter(p -> p.tags["arch"]!="x86_64", supported_platforms())

# Expand for microarchitectures on x86_64 (library doesn't have CPU dispatching)
# Tests on Linux/x86_64 yielded a slow binary with avx512 for some reason, so disable that
platforms_x86_64 = expand_microarchitectures(platforms_x86_64)
platforms_x86_64 = filter(p -> p.tags["march"] != "avx512", platforms_x86_64)

platforms = vcat(platforms_x86_64, platforms_other)
# The products that we will ensure are always built
products = [
    LibraryProduct("libfinufft", :libfinufft)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FFTW_jll", uuid="f5851436-0d7a-5f13-b9de-f02708fd171a"))
    Dependency(PackageSpec(name="CompilerSupportLibraries_jll", uuid="e66e0078-7015-5450-92f7-15fbd957f2ae"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8", julia_compat=julia_compat)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# other optimized string routines copied from julia/base/strings/util.jl

function Base.startswith(a::Union{DenseStringViewAndSub,StringAndSub},
                         b::Union{DenseStringViewAndSub,StringAndSub})
    cub = ncodeunits(b)
    if ncodeunits(a) < cub
        false
    elseif _memcmp(a, b, sizeof(b)) == 0
        nextind(a, cub) == cub + 1
    else
        false
    end
end

function Base.endswith(a::Union{DenseStringViewAndSub,StringAndSub},
                       b::Union{DenseStringViewAndSub,StringAndSub})
    cub = ncodeunits(b)
    astart = ncodeunits(a) - ncodeunits(b) + 1
    if astart < 1
        false
    elseif GC.@preserve(a, _memcmp(pointer(a, astart), b, sizeof(b))) == 0
        thisind(a, astart) == astart
    else
        false
    end
end

function Base.chomp(s::StringViewAndSub)
    i = lastindex(s)
    if i < 1 || codeunit(s,i) != 0x0a
        return @inbounds SubString(s, 1, i)
    elseif i < 2 || codeunit(s,i-1) != 0x0d
        return @inbounds SubString(s, 1, prevind(s, i))
    else
        return @inbounds SubString(s, 1, prevind(s, i-1))
    end
end

Base.replace(str::DenseStringViewAndSub, pat_repl::Pair{<:AbstractChar}; count::Integer=typemax(Int)) =
    replace(str, isequal(first(pat_repl)) => last(pat_repl); count=count)

Base.replace(str::DenseStringViewAndSub, pat_repl::Pair{<:Union{Tuple{Vararg{<:AbstractChar}},
                                            AbstractVector{<:AbstractChar},Set{<:AbstractChar}}};
        count::Integer=typemax(Int)) =
    replace(str, in(first(pat_repl)) => last(pat_repl), count=count)

import Base: _pat_replacer, _free_pat_replacer

function Base.replace(str::DenseStringViewAndSub, pat_repl::Pair; count::Integer=typemax(Int))
    pattern, repl = pat_repl
    count == 0 && return str
    count < 0 && throw(DomainError(count, "`count` must be non-negative."))
    n = 1
    e = lastindex(str)
    i = a = firstindex(str)
    pattern = _pat_replacer(pattern)
    r = something(findnext(pattern,str,i), 0)
    j, k = first(r), last(r)
    if j == 0
        _free_pat_replacer(pattern)
        return str
    end
    out = IOBuffer(sizehint=floor(Int, 1.2sizeof(str)))
    while j != 0
        if i == a || i <= k
            GC.@preserve str unsafe_write(out, pointer(str, i), UInt(j-i))
            Base._replace(out, repl, str, r, pattern)
        end
        if k < j
            i = j
            j > e && break
            k = nextind(str, j)
        else
            i = k = nextind(str, k)
        end
        r = something(findnext(pattern,str,k), 0)
        r === 0:-1 || n == count && break
        j, k = first(r), last(r)
        n += 1
    end
    _free_pat_replacer(pattern)
    write(out, SubString(str,i))
    String(take!(out))
end


text:
using Test

using ErlPort.ErlTerms


text:
################################################
"""
$(TYPEDEF)

Data structure containing DenseSystem used to calculate
test functions for boundary flux calculations.

$(TYPEDFIELDS)
"""
mutable struct TestFunctionFactory{Tv}

    """
    Original system
    """
    system::AbstractSystem{Tv}

    """
    Test function system
    """
    tfsystem::DenseSystem{Tv}
end

################################################
"""
$(TYPEDSIGNATURES)

Constructor for TestFunctionFactory from System
"""
function TestFunctionFactory(system::AbstractSystem{Tv}) where Tv
    physics=Physics(
        flux=function(f,u,edge)
        f[1]=u[1]-u[2]
        end,
        storage=function(f,u,node)
        f[1]=u[1]
        end
    )
    tfsystem=System(system.grid,physics,unknown_storage=:dense)
    enable_species!(tfsystem,1,[i for i=1:num_cellregions(system.grid)])
    return TestFunctionFactory{Tv}(system,tfsystem)
end

############################################################################
"""
$(TYPEDSIGNATURES)

Create testfunction which has Dirichlet zero boundary conditions  for boundary
regions in bc0 and Dirichlet one boundary conditions  for boundary
regions in bc1.
"""
function testfunction(factory::TestFunctionFactory{Tv}, bc0, bc1) where Tv

    u=unknowns(factory.tfsystem)
    f=unknowns(factory.tfsystem)
    u.=0
    f.=0

    factory.tfsystem.boundary_factors.=0
    factory.tfsystem.boundary_values.=0

    for i=1:length(bc1)
        factory.tfsystem.boundary_factors[1,bc1[i]]=Dirichlet
        factory.tfsystem.boundary_values[1,bc1[i]]=-1
    end

    for i=1:length(bc0)
        factory.tfsystem.boundary_factors[1,bc0[i]]=Dirichlet
        factory.tfsystem.boundary_values[1,bc0[i]]=0
    end

    eval_and_assemble(factory.tfsystem,u,u,f,Inf)

    _initialize!(u,factory.tfsystem)
    lufact=LinearAlgebra.lu(factory.tfsystem.matrix)
    ldiv!(vec(u),lufact,vec(f))
    return vec(u)
end

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

"""
$(SIGNATURES)

Calculate test function integral for transient solution.
"""
function integrate(system::AbstractSystem{Tv,Ti},tf::Vector{Tv},U::AbstractMatrix{Tv},
                   Uold::AbstractMatrix{Tv}, tstep::Real) where {Tv,Ti}
    grid=system.grid
    nspecies=num_species(system)
    integral=zeros(Tv,nspecies)
    res=zeros(Tv,nspecies)
    src=zeros(Tv,nspecies)
    stor=zeros(Tv,nspecies)
    storold=zeros(Tv,nspecies)
    tstepinv=1.0/tstep

    physics=system.physics
    node=Node{Tv,Ti}(system)
    bnode=BNode{Tv,Ti}(system)
    edge=Edge{Tv,Ti}(system)
    @create_physics_wrappers(physics,node,bnode,edge)

    UKL=Array{Tv,1}(undef,2*nspecies)
    UK=Array{Tv,1}(undef,nspecies)
    UKold=Array{Tv,1}(undef,nspecies)

    geom=grid[CellGeometries][1]
    csys=grid[CoordinateSystem]
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]

    if haskey(grid,CellEdges)
        cellx=grid[CellEdges]
        edgenodes=grid[EdgeNodes]
        has_celledges=true
    else
        cellx=grid[CellNodes]
        edgenodes=local_celledgenodes(geom)
        has_celledges=false
    end

    for icell=1:num_cells(grid)
        for iedge=1:num_edges(geom)
            _fill!(edge,cellx,edgenodes,cellregions,iedge,icell, has_celledges)

            @views UKL[1:nspecies].=U[:,edge.node[1]]
            @views UKL[nspecies+1:2*nspecies].=U[:,edge.node[2]]

            res.=zero(Tv)
            fluxwrap(res,UKL)
            for ispec=1:nspecies
                if system.node_dof[ispec,edge.node[1]]==ispec && system.node_dof[ispec,edge.node[2]]==ispec
                    integral[ispec]+=system.celledgefactors[iedge,icell]*res[ispec]*(tf[edge.node[1]]-tf[edge.node[2]])
                end
            end
        end

        for inode=1:num_nodes(geom)
            _fill!(node,cellnodes,cellregions,inode,icell)
            begin
                res.=zero(Tv)
                stor.=zero(Tv)
                storold.=zero(Tv)
                for ispec=1:nspecies
                    UK[ispec]=U[ispec,node.index]
                    UKold[ispec]=Uold[ispec,node.index]
                end
                reactionwrap(res,UK)
                sourcewrap(src)
                storagewrap(stor,UK)
                storagewrap(storold,UKold)
            end
            for ispec=1:nspecies
                if system.node_dof[ispec,node.index]==ispec
                    integral[ispec]+=system.cellnodefactors[inode,icell]*(res[ispec]-src[ispec]+(stor[ispec]-storold[ispec])*tstepinv)*tf[node.index]
                end
            end
        end
    end
    return integral
end

############################################################################
"""
$(SIGNATURES)

Calculate test function integral for steady state solution.
"""
integrate(system::AbstractSystem{Tv,Ti},tf::Vector{Tv},U::AbstractMatrix{Tv}) where {Tv,Ti} =integrate(system,tf,U,U,Inf)

############################################################################
"""
$(SIGNATURES)

Steady state part of test function integral.
"""
function integrate_stdy(system::AbstractSystem{Tv,Ti},tf::Vector{Tv},U::AbstractArray{Tu,2}) where {Tu,Tv,Ti}
    grid=system.grid
    nspecies=num_species(system)
    integral=zeros(Tu,nspecies)
    res=zeros(Tu,nspecies)
    src=zeros(Tu,nspecies)
    stor=zeros(Tu,nspecies)

    physics=system.physics
    node=Node{Tv,Ti}(system)
    bnode=BNode{Tv,Ti}(system)
    edge=Edge{Tv,Ti}(system)
    @create_physics_wrappers(physics,node,bnode,edge)

    UKL=Array{Tu,1}(undef,2*nspecies)
    UK=Array{Tu,1}(undef,nspecies)
    geom=grid[CellGeometries][1]
    csys=grid[CoordinateSystem]
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]

    if haskey(grid,CellEdges)
        cellx=grid[CellEdges]
        edgenodes=grid[EdgeNodes]
        has_celledges=true
    else
        cellx=grid[CellNodes]
        edgenodes=local_celledgenodes(geom)
        has_celledges=false
    end

    for icell=1:num_cells(grid)
        for iedge=1:num_edges(geom)
            _fill!(edge,cellx,edgenodes,cellregions,iedge,icell, has_celledges)

            @views UKL[1:nspecies].=U[:,edge.node[1]]
            @views UKL[nspecies+1:2*nspecies].=U[:,edge.node[2]]
            res.=zero(Tv)

            fluxwrap(res,UKL)
            for ispec=1:nspecies
                if system.node_dof[ispec,edge.node[1]]==ispec && system.node_dof[ispec,edge.node[2]]==ispec
                    integral[ispec]+=system.celledgefactors[iedge,icell]*res[ispec]*(tf[edge.node[1]]-tf[edge.node[2]])
                end
            end
        end

        for inode=1:num_nodes(geom)
            _fill!(node,cellnodes,cellregions,inode,icell)

            res.=zeros(Tv)
            src.=zeros(Tv)
            @views UK.=U[:,node.index]
            reactionwrap(res,UK)
            sourcewrap(src)
            for ispec=1:nspecies
                if system.node_dof[ispec,node.index]==ispec
                    integral[ispec]+=system.cellnodefactors[inode,icell]*(res[ispec]-src[ispec])*tf[node.index]
                end
            end
        end
    end
    return integral
end

############################################################################
"""
$(SIGNATURES)

Calculate transient part of test function integral.
"""
function integrate_tran(system::AbstractSystem{Tv,Ti},tf::Vector{Tv},U::AbstractArray{Tu,2}) where {Tu,Tv,Ti}
    grid=system.grid
    nspecies=num_species(system)
    integral=zeros(Tu,nspecies)
    res=zeros(Tu,nspecies)
    stor=zeros(Tu,nspecies)

    physics=system.physics
    node=Node{Tv,Ti}(system)
    bnode=BNode{Tv,Ti}(system)
    edge=Edge{Tv,Ti}(system)
    @create_physics_wrappers(physics,node,bnode,edge)

    UK=Array{Tu,1}(undef,nspecies)
    geom=grid[CellGeometries][1]
    csys=grid[CoordinateSystem]
    coord=grid[Coordinates]
    cellnodes=grid[CellNodes]
    cellregions=grid[CellRegions]

    for icell=1:num_cells(grid)
        for inode=1:num_nodes(geom)
            _fill!(node,cellnodes,cellregions,inode,icell)

            res.=zeros(Tv)
            stor.=zeros(Tv)
            @views UK.=U[:,node.index]

            storagewrap(stor,U[:,node.index])
            for ispec=1:nspecies
                if system.node_dof[ispec,node.index]==ispec
                    integral[ispec]+=system.cellnodefactors[inode,icell]*stor[ispec]*tf[node.index]
                end
            end
        end
    end
    return integral
end


text:
# module Stats
# using lectura
# Tabla::DataFrame
# numero= "23"
using DataFrames
push!(LOAD_PATH,"../src/") #eso es por si Julia no encuentra el ./src
Tabla = DataFrame()
function Stats(Path::String,Tipo::String,Formato="csv")
  # global Tabla = lectura.lecturas(Path,Tipo,Formato)
end
function consulta(queries::Vector{String})
  print(numero)
end
# end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to display printing and manipulation.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _draw_continuation_row(display::Display, io::IO, tf::TextFormat, text_crayon::Crayon, border_crayon::Crayon, cols_width::Vector{Int}, vlines::Vector{Int}, alignment::Symbol)

Draw the continuation row when the table has filled the vertical space
available. This function prints in each column the character `⋮` with the
alignment in `alignment`.

"""
function _draw_continuation_row(display::Display, io::IO, tf::TextFormat,
                                text_crayon::Crayon, border_crayon::Crayon,
                                cols_width::Vector{Int}, vlines::Vector{Int},
                                alignment::Symbol)

    # In case of a continuation row, we want the last character to indicate that
    # the table continues both vertically and horizontally in case the text is
    # cropped.
    old_cont_char = display.cont_char
    display.cont_char = '⋱'

    num_cols = length(cols_width)

    0 ∈ vlines && _p!(display, border_crayon, tf.column)

    @inbounds for j = 1:num_cols
        data_ij_str, data_ij_len = _str_aligned("⋮", alignment, cols_width[j])
        data_ij_str  = " " * data_ij_str * " "
        data_ij_len += 2

        flp = j == num_cols

        _p!(display, text_crayon, data_ij_str, false, data_ij_len)
        _pc!(j ∈ vlines, display, border_crayon, tf.column, "", flp, 1, 0)
        _eol(display) && break
    end

    _nl!(display, io)

    display.cont_char = old_cont_char

    return nothing
end

"""
    _draw_line!(display::Display, io::IO, left::Char, intersection::Char, right::Char, row::Char, border_crayon::Crayon, cols_width::Vector{Int}, vlines::Vector{Int})

Draw a vertical line in internal line buffer of `display` and then flush to the
io `io`.

"""
function _draw_line!(display::Display, io::IO, left::Char, intersection::Char,
                     right::Char, row::Char, border_crayon::Crayon,
                     cols_width::Vector{Int}, vlines::Vector{Int})

    # We does not want to add ellipsis when drawing lines.
    old_cont_char           = display.cont_char
    old_cont_space_char     = display.cont_space_char
    display.cont_char       = row
    display.cont_space_char = row

    num_cols = length(cols_width)

    0 ∈ vlines && _p!(display, border_crayon, left)

    @inbounds for i = 1:num_cols
        # Check the alignment and print.
        _p!(display, border_crayon, row^(cols_width[i]+2)) && break

        i != num_cols &&
            _pc!(i ∈ vlines, display, border_crayon, intersection, "")
    end

    _pc!(num_cols ∈ vlines, display, border_crayon, right, "", true)
    _nl!(display, io)

    display.cont_char       = old_cont_char
    display.cont_space_char = old_cont_space_char

    return nothing
end

"""
    _eol(display::Display)

Return `true` if the cursor is at the end of line or `false` otherwise.

"""
_eol(display::Display) =
    (display.size[2] > 0) && (display.col >= display.size[2])

"""
    _nl!(display::Display, io::IO)

Flush the internal line buffer of `display` into `io`.

"""
function _nl!(display::Display, io::IO)
    # Update the information about the current columns and row of the display.
    display.row += 1
    display.col  = 0

    # Flush the current line to the buffer removing any trailing space.
    str = rstrip(String(take!(display.buf_line)))
    println(io, str)
end

"""
    _p!(display::Display, crayon::Crayon, str::Char, final_line_print::Bool = false, lstr::Int = -1)
    _p!(display::Display, crayon::Crayon, str::String, final_line_print::Bool = false, lstr::Int = -1)

Print `str` into the internal line buffer of `display` using the Crayon `crayon`
with the display information in `display`. The parameter `final_line_print` must
be set to `true` if this is the last string that will be printed in the line.
This is necessary for the algorithm to select whether or not to include the
continuation character.

The size of the string can be passed to `lstr` to save computational burden. If
`lstr = -1`, then the string length will be computed inside the function.

The line buffer can be flushed to an `io` using the function `_nl!`.

"""
_p!(display::Display, crayon::Crayon, str::Char, final_line_print::Bool = false,
    lstr::Int = -1) = _p!(display, crayon, string(str), final_line_print, lstr)

function _p!(display::Display, crayon::Crayon, str::String,
             final_line_print::Bool = false, lstr::Int = -1)

    # Get the size of the string if required.
    lstr < 0 && (lstr = textwidth(str))

    # `sapp` is a string to be appended to `str`. This is used to add `⋯` if the
    # text must be wrapped. Notice that `lapp` is the length of `sapp`.
    sapp = ""
    lapp = 0

    # When printing a line, we must verify if the display has bounds. This is
    # done by looking if the horizontal size is positive.
    @inbounds if display.size[2] > 0

        # If we are at the end of the line, then just return.
        _eol(display) && return true

        Δ = display.size[2] - (lstr + display.col)

        # Check if we can print the entire string.
        if Δ <= 0
            # If we cannot, then create a wrapped string considering how many
            # columns are left.
            if lstr + Δ - 2 > 0
                # Here we crop the string considering the available space
                # reserving 2 characters for the line continuation indicator.
                str  = _crop_str(str, lstr + Δ - 2)
                lstr = lstr + Δ - 2
                sapp = display.cont_space_char * display.cont_char
                lapp = 2
            elseif display.size[2] - display.col == 2
                # If there are only 2 characters left, then we must only print
                # " ⋯".
                str  = ""
                lstr = 0
                sapp = display.cont_space_char * display.cont_char
                lapp = 2
            elseif display.size[2] - display.col == 1
                # If there are only 1 character left, then we must only print
                # "⋯".
                str  = ""
                lstr = 0
                sapp = string(display.cont_char)
                lapp = 1
            else
                # This should never be reached.
                @error("Internal error!")
                return true
            end

            # In this case, we reached the end of display. Thus, remove any
            # trailing spaces.
            sapp = rstrip(sapp)
            length(sapp) == 0 && (str = rstrip(str))

        elseif !final_line_print
            # Here we must verify if this is the final printing on this line. If
            # it is, then we should just check if the entire string fits on the
            # available size. Otherwise, we must see if, after printing the
            # current string, we will have more than 1 space left. If not, then
            # we just add the continuation character sequence.

            if Δ == 1
                str   = lstr > 1 ? _crop_str(str, lstr - 1) : ""
                lstr -= 1
                sapp  = display.cont_space_char * display.cont_char
                lapp  = 2

                # In this case, we reached the end of display. Thus, remove any
                # trailing spaces.
                sapp = rstrip(sapp)
                length(sapp) == 0 && (str = rstrip(str))
            end
        end
    end

    # Print the with correct formating.
    #
    # Notice that all text printed with crayon is reset right after the string.
    # Hence, if the crayon is empty (`_default_crayon`) or if it is a reset,
    # then we can just print as if the terminal does not support color.
    if (crayon != _default_crayon) && (crayon != _reset_crayon) && display.has_color
        print(display.buf_line, crayon, str, _reset_crayon, sapp)
    else
        print(display.buf_line, str, sapp)
    end

    # Update the current columns.
    display.col += lstr + lapp

    # Return if we reached the end of line.
    return _eol(display)
end

"""
    _pc!(cond::Bool, display::Display, io::IO, crayon::Crayon, str_true::Union{Char,String}, str_false::Union{Char,String}, final_line_print::Bool = false, lstr_true::Int = -1, lstr_false::Int = -1)

If `cond == true` then print `str_true`. Otherwise, print `str_false`. Those
strings will be printed into the internal line buffer of `display` using the
Crayon `crayon` with the display information in `display`. The parameter
`final_line_print` must be set to `true` if this is the last string that will be
printed in the line. This is necessary for the algorithm to select whether or
not to include the continuation character.

The size of the strings can be passed to `lstr_true` and `lstr_false` to save
computational burden. If they are `-1`, then the string lengths will be computed
inside the function.

"""
function _pc!(cond::Bool, display::Display, crayon::Crayon,
              str_true::Union{Char,String}, str_false::Union{Char,String},
              final_line_print::Bool = false, lstr_true::Int = -1,
              lstr_false::Int = -1)
    if cond
        return _p!(display, crayon, str_true, final_line_print, lstr_true)
    else
        return _p!(display, crayon, str_false, final_line_print, lstr_false)
    end
end


text: println("Hello,World!!")#>printing hello world

max_stars_repo_name: ashwani-rathee/Audio

max_issues_repo_name: ashwani-rathee/Audio

max_forks_repo_name: ashwani-rathee/Audio


text:
cube, angles = BetaPictoris[:cube, :pa]

@testset "Interface" begin
    @test GreeDS().kernel == PCA()
    @test GreeDS().threshold == 0
    @test GreeDS(7).kernel == PCA(7)
end

@testset "Decomposition" begin

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        A, w = @inferred ADI.fit(GreeDS(PCA(ncomps=N)), cube; angles=angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (size(cube, 1), N)
    end

    # default is to use whole cube
    S = reconstruct(GreeDS(), cube; angles=angles)
    @test size(S) == size(cube)
end

@testset "ADI Trivial" begin
    data = ones(10, 100, 100)
    angs = zeros(10)

    reduced_5 = GreeDS(PCA(5))(data, angs)
    reduced_10 = GreeDS(PCA(10))(data, angs)

    @test size(reduced_5) == size(reduced_10) == (100, 100)
    @test all(x -> isapprox(x, 0, atol=1e-9), reduced_5)
    @test all(x -> isapprox(x, 0, atol=1e-9), reduced_10)
end

@testset "RDI Trivial" begin

    S = subtract(GreeDS(1), cube; angles=angles, ref=zero(cube))
    @test S ≈ cube rtol=2e-1
end

@testset "Decomposition - TPCA" begin
    # get sizes correct for ncomps
    for N in [1, 3, 5]
        A, w = ADI.fit(GreeDS(TPCA(N)), cube; angles=angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (size(cube, 1), N)
    end

    # default is to use whole cube
    S = reconstruct(GreeDS(TPCA()), cube; angles=angles)
    @test size(S) == size(cube)
end

@testset "Decomposition - NMF" begin
    data = ones(10, 100, 100)
    angs = zeros(10)

    reduced_8 = GreeDS(NMF(8))(data, angs)
    reduced_10 = GreeDS(NMF(10))(data, angs)

    @test size(reduced_8) == size(reduced_10) == (100, 100)
    @test all(x -> isapprox(x, 0, atol=1e-8), reduced_8)
    @test all(x -> isapprox(x, 0, atol=1e-7), reduced_10)
end


text:
module API
using ..Telegram

include("telegram_api.jl")

for (f, fdoc) in TELEGRAM_API
    fs=string(f)
    @eval begin
        function $f(client::Telegram.TelegramClient = Telegram.DEFAULT_OPTS.client; kwargs...)
            params = Dict{Symbol, Any}(kwargs)
            params[:chat_id] = get(params, :chat_id, Telegram.chatid(client))
            params[:parse_mode] = get(params, :parse_mode, client.parse_mode)
            Telegram.query(client, $fs, params = params)
        end

        @doc $fdoc $f
        export $f
    end
end

end # module