text:
function Base.showarg(io::IO, A::BSplineInterpolation{T,N,TW,ST}, toplevel) where {T,N,TW,ST}
    print(io, "interpolate(")
    Base.showarg(io, A.coefs, false)
    print(io, ", ")
    show(io, itpflag(A))
    if toplevel
        print(io, ") with element type ",T)
    else
        print(io, ')')
    end
end

function Base.showarg(io::IO, A::GriddedInterpolation{T,N,TC,ST,K}, toplevel) where {T,N,TC,ST,K}
    print(io, "interpolate(")
    _showknots(io, A.knots)
    print(io, ", ")
    Base.showarg(io, A.coefs, false)
    print(io, ", ")
    show(io, itpflag(A))
    if toplevel
        print(io, ") with element type ",T)
    else
        print(io, ')')
    end
end

_showknots(io, A) = Base.showarg(io, A, false)
function _showknots(io, tup::NTuple{N,Any}) where N
    print(io, '(')
    for (i, A) in enumerate(tup)
        if A isa AbstractRange
            show(io, A)
        else
            Base.showarg(io, A, false)
        end
        i < N && print(io, ',')
    end
    N == 1 && print(io, ',')
    print(io, ')')
end

function Base.showarg(io::IO, A::ScaledInterpolation{T}, toplevel) where {T}
    print(io, "scale(")
    Base.showarg(io, A.itp, false)
    print(io, ", ", A.ranges, ')')
    if toplevel
        print(io, " with element type ",T)
    end
end

function Base.showarg(io::IO, A::MonotonicInterpolation{T, TCoeffs, Tel, Type, K}, toplevel) where {T, TCoeffs, Tel, Type, K}
    print(io, "interpolate(")
    _showknots(io, A.knots)
    print(io, ", ")
    Base.showarg(io, A.A, false)
    print(io, ", ")
    show(io, A.it)
    if toplevel
        print(io, ") with element type ",T)
    else
        print(io, ')')
    end
end

function Base.showarg(io::IO, A::Extrapolation{T,N,TI,IT,ET}, toplevel) where {T,N,TI,IT,ET}
    print(io, "extrapolate(")
    Base.showarg(io, A.itp, false)
    print(io, ", ")
    show(io, etpflag(A))
    print(io, ')')
    if toplevel
        print(io, " with element type ",T)
    end
end

function Base.showarg(io::IO, A::FilledExtrapolation{T,N,TI,IT}, toplevel) where {T,N,TI,IT}
    print(io, "extrapolate(")
    Base.showarg(io, A.itp, false)
    print(io, ", ", A.fillvalue, ')')
    if toplevel
        print(io, " with element type ",T)
    end
end

"""
    show_ranged(io, X, knots)

A replacement for the default array-`show` for types that may not have the canonical evaluation points.
`rngs` is the tuple of knots along each axis.
"""
function show_ranged(io::IO, X, knots)
    summary(io, X)
    isempty(X) && return
    print(io, ":")
    if !(haskey(io, :compact)) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if get(io, :limit, false) && eltype(X) === Method
        io = IOContext(io, :limit => false)
    end
    if get(io, :limit, false) && (displaysize(io))[1] - 4 <= 0
        return print(io, " …")
    else
        println(io)
    end
    io = IOContext(io, :typeinfo => eltype(X))
    Base.print_array(io, [X(x...) for x in Iterators.product(knots...)])
end

getknots(X::BSplineInterpolation)  = axes(X)
getknots(X::ScaledInterpolation)   = X.ranges
getknots(X::GriddedInterpolation)  = X.knots
getknots(X::AbstractExtrapolation) = getknots(parent(X))

Base.show(io::IO, ::MIME{Symbol("text/plain")}, X::ScaledInterpolation)   = show_ranged(io, X, getknots(X))
Base.show(io::IO, ::MIME{Symbol("text/plain")}, X::GriddedInterpolation)  = show_ranged(io, X, getknots(X))
Base.show(io::IO, ::MIME{Symbol("text/plain")}, X::AbstractExtrapolation) = show_ranged(io, X, getknots(X))

# these are called by show(X)
Base.show(io::IO, X::ScaledInterpolation)   = show_ranged(io, X, getknots(X))
Base.show(io::IO, X::GriddedInterpolation)  = show_ranged(io, X, getknots(X))
Base.show(io::IO, X::AbstractExtrapolation) = show_ranged(io, X, getknots(X))


text:
using ElectronDisplay
ElectronDisplay.CONFIG.showable = showable
ElectronDisplay.CONFIG.single_window = true
ElectronDisplay.CONFIG.focus = false
using ImageMagick, FileIO
using WGLMakie, AbstractPlotting, JSServe, Test
using Pkg

# ImageIO seems broken on 1.6 ... and there doesn't
# seem to be a clean way anymore to force not to use a loader library?
filter!(x-> x !== :ImageIO, FileIO.sym2saver[:PNG])
filter!(x-> x !== :ImageIO, FileIO.sym2loader[:PNG])
AbstractPlotting.set_theme!(resolution=(400, 400))

path = normpath(joinpath(dirname(pathof(AbstractPlotting)), "..", "test", "ReferenceTests"))
Pkg.develop(PackageSpec(path = path))
using ReferenceTests
using ReferenceTests: nice_title
excludes = Set([
    "Streamplot animation",
    "Transforming lines",
    "image scatter",
    "Line GIF",
    "surface + contour3d",
    # Hm weird, looks like some internal JSServe error missing an Observable:
    "Errorbars x y low high",
    "Rangebars x y low high",
    # These are a bit sad, since it's just missing interpolations
    "FEM mesh 2D",
    "FEM polygon 2D",
    # missing transparency & image
    "Wireframe of a Surface",
    "Image on Surface Sphere",
    "Surface with image",
    # Marker size seems wrong in some occasions:
    "Hbox",
    "UnicodeMarker",
    # Not sure, looks pretty similar to me! Maybe blend mode?
    "Test heatmap + image overlap",
    "Stars"
])

database = ReferenceTests.load_database()
filter!(database) do (name, entry)
    !(entry.title in excludes) &&
    nice_title(entry) !== "short_tests_83" &&
    nice_title(entry) !== "short_tests_78"
end
recorded = joinpath(@__DIR__, "recorded")
rm(recorded; force=true, recursive=true); mkdir(recorded)
ReferenceTests.record_tests(database; recording_dir=recorded)
ReferenceTests.reference_tests(recorded; difference=0.06)

ext: jl


text:
@testset "901.online-stock-span.jl" begin
    s = StockSpanner()
    @test next!(s, 100) == 1
    @test next!(s, 80) == 1
    @test next!(s, 60) == 1
    @test next!(s, 70) == 2
    @test next!(s, 60) == 1
    @test next!(s, 75) == 4
    @test next!(s, 85) == 6
end


text:
function ssdDistance(referenceImage::regImage,templateImage::regImage,
                     transformedGrid::Array{Float64,1};
                     doDerivative::Bool=false,doHessian::Bool=false,options::regOptions=regOptions(),centeredGrid::Array{Float64,1}=zeros(1))

  # get relevant options
  parametricOnly = options.parametricOnly
  matrixFree = options.matrixFree

  # check centered grid
  if(size(centeredGrid)[1]==1)
    centeredGrid = getCellCenteredGrid(referenceImage).data
  end

  # check grid
  if(checkForOddNumberOfGridPoints(referenceImage.data,transformedGrid))
    error("StaggeredGrids are not supported. Please transform to cell centered first.")
  end

  # interpolation of the template image at transformed grid points
  transformedImage, dX_transformedImage, dY_transformedImage =
      interpolateImage(templateImage,transformedGrid,doDerivative=true)

  # measure the ssd distance
  N = prod(size(referenceImage.data)) # product of sizes

  residual = Array(transformedImage .- referenceImage.data)[:]
  prodH = prod(referenceImage.voxelsize)
  functionValue = 0.5 * prodH * residual' * residual
  # calculate ssd derivatives matrix free?
  if(matrixFree)
    dFunctionValue,d2FunctionValue = ssdDerivativesMatrixFree(dX_transformedImage,dY_transformedImage,residual,centeredGrid,prodH,N,parametricOnly,doDerivative,doHessian)
  else
    dFunctionValue,d2FunctionValue = ssdDerivativesMatrixBased(dX_transformedImage,dY_transformedImage,residual,centeredGrid,prodH,N,parametricOnly,doDerivative,doHessian)
  end

  return [functionValue,dFunctionValue,d2FunctionValue,(dX_transformedImage,dY_transformedImage)]

end

function ssdDerivativesMatrixBased(dX_transformedImage::Array{Float64,1},
                                   dY_transformedImage::Array{Float64,1},
                                   residual::Array{Float64,1},
                                   centeredGrid::Array{Float64,1},
                                   pixelSize::Number,N::Number,
                                   parametricOnly::Bool,doDerivative::Bool,doHessian::Bool)

  if(doDerivative)
    dTransformedImage = spdiagm((dX_transformedImage,dY_transformedImage),[0,N])
    if(parametricOnly)
      Q = sparse([centeredGrid[1:N] centeredGrid[N+1:end] ones(N)])
      Q = [Q spzeros(size(Q)[1],size(Q)[2]); spzeros(size(Q)[1],size(Q)[2]) Q]
      dTransformedImage = dTransformedImage * Q
    end
    dFunctionValue = pixelSize .* dTransformedImage' * residual
  else
    dFunctionValue = 0.0
    dTransformedImage = 0.0
  end

  if(doHessian)
    d2FunctionValue = pixelSize .* dTransformedImage' * dTransformedImage
  else
    d2FunctionValue = 0.0
  end

  return dFunctionValue, d2FunctionValue

end

function ssdDerivativesMatrixFree(dX_transformedImage::Array{Float64,1},
                                  dY_transformedImage::Array{Float64,1},
                                  residual::Array{Float64,1},
                                  centeredGrid::Array{Float64,1},
                                  pixelSize::Number,N::Number,
                                  parametricOnly::Bool,doDerivative::Bool,doHessian::Bool)

  if(doDerivative)
    dFunctionValue = BLAS.scal(2*N,pixelSize,dTtransposedMultiplication(N,dX_transformedImage,dY_transformedImage,residual),1)
    if(parametricOnly)
      dFunctionValue = QtransposedMultiplication(N,dFunctionValue,centeredGrid)
    end
  else
    dFunctionValue = 0.0
    dTransformedImage = 0.0
  end

  if(doHessian)
    if(parametricOnly)
      d2FunctionValue = computeParametricHessian(N,pixelSize,dX_transformedImage,dY_transformedImage,centeredGrid)
    else
      d2FunctionValue(x) = hessianFunction(N,pixelSize,dX_transformedImage,dY_transformedImage,x)
    end
  else
    d2FunctionValue = 0.0
  end

  return dFunctionValue,d2FunctionValue

end


text:
"""

                newtonroot(f[, f']; x0 = 0.8, tolerance = 1.e-13, maxiter = 1000)

This function computes the local root of a univariate smooth function `f(x)` using the Newton's Method. This method starts with a  `x_0` guess, a function `f(x)` and the first-derivative  `f'(x)`. The algorithm follows: `` x^{n+1} = x^n - f(x^n)/f'(x^n) `` until `` |x^{n+1} - x^n|`` is below the tolerance level.

...
# Arguments
- `f(x)::Function`: the function you are trying to solve.
- `f'(x)::Function`: the first derivative of the function you are trying to solve. If it is not given the program estimate `f'(x)` using ForwardDiff.
- `x0::Float64`: the initial value for the algorithm. The default value is 0.8.
- `tolerance::Float64` : the minimum tolerance for the algorithm to converge. The default value is 1.0e-13.
- `maxiter::Int64`: the maximum number of iterations allowed. The default value is 1000.
---
# Examples
```julia-repl
julia> f(x)  =  (x - 2)^2
f (generic function with 1 method)

julia> df(x) = 2(x - 2)
df (generic function with 1 method)

julia> newtonroot(f,df)
(value = 1.9999999999999318, normdiff = 6.816769371198461e-14, iter = 45)

julia> a = newtonroot(f)[1]
1.9999999999999318

julia> newtonroot(f)[1] ≈ 2
true
```
...

"""
function newtonroot(f, f_prime; x0=0.8, tolerance=1.0E-13, maxiter=1000)
    # Initial values
    v_old = x0
    normdiff = Inf
    iter = 1

    # Running the algorithm
    while (normdiff > tolerance)  && (iter < maxiter)
        if !(f_prime(v_old) == 0)
            v_new = v_old - (f(v_old) / f_prime(v_old))
            normdiff = norm(v_new - v_old)
            # replace and continue
            v_old = v_new
            iter = iter + 1
        else
            return(value = nothing, normdiff = nothing, iter = nothing)
        end
    end
    if (iter < maxiter)
        return(value = v_old, normdiff = normdiff, iter = iter)
    else
        return(value = nothing, normdiff = nothing, iter = nothing)
    end
end

function newtonroot(f; x0=0.8, tolerance=1.0E-13, maxiter=1000)
    # Initial values
    v_old = x0
    normdiff = Inf
    iter = 1

    # Taking the derivative
    D(f) = x -> ForwardDiff.derivative(f, x)
    f_prime = D(f)

    # Running the algorithm

    while (normdiff > tolerance)  && (iter < maxiter)
        if !(f_prime(v_old) == 0)
            v_new = v_old - (f(v_old) / f_prime(v_old))
            normdiff = norm(v_new - v_old)
            # replace and continue
            v_old = v_new
            iter = iter + 1
        else
            return(value = nothing, normdiff = nothing, iter = nothing)
        end
    end
    if (iter < maxiter)
        return(value = v_old, normdiff = normdiff, iter = iter)
    else
        return(value = nothing, normdiff = nothing, iter = nothing)
    end
end


text:
"""
    MovementModel{T<:Real}
        M::Matrix

Model for movement preferences.
"""
struct MovementModel{T<:Real, Ti<:Integer}
    M::Matrix{T}
    Ωsize::Tuple{Ti, Ti}
end

function MovementModel(domain::AbstractFisheryDomain,
                       habitat::AbstractMatrix,
                       pref_fn::Function,
                       dist_fn::Function)
    mvt_n = prod(domain.n)
    mvt = Matrix{eltype(habitat)}(undef, mvt_n, mvt_n)
    for jdx in 1:mvt_n, idx in 1:mvt_n
        mvt[idx, jdx] = pref_fn(habitat[idx]) *
            dist_fn(domain.distances[idx, jdx])
    end
    mvt ./= sum(mvt; dims = 1)
    MovementModel(mvt, domain.n)
end

function MovementModel(B::Bathymetry,
                       distance::Distributions.UnivariateDistribution)
    mvt = pdf.(dist, B.Ω.distances)
    mvt ./= sum(mvt, 2)
    MovementModel(Matrix(mvt'), size(B.Ω))
end

function MovementModel(B::Bathymetry,
                       distance::Distributions.UnivariateDistribution,
                       depth::Distributions.UnivariateDistribution)
    distmvt = pdf.(distance, B.Ω.distances)
    dpthmvt = pdf.(depth, vec(B.bathymetry))
    mvt = distmvt .* dpthmvt'
    mvt ./= sum(mvt; dims = 2)
    MovementModel(Matrix(mvt'), size(B.Ω))
end

function (M::MovementModel)(P::PopState)
    p = M.M * vecstate(P)
    PopState(reshape(p, M.Ωsize...))
end

"""
    eqdist(M::MovementModel)

Find the equilibrium distribution under a given movement model. Uses Arpack.jl
to find the eigenvector associated with the largest eigenvalue (which will be
one due to the row-standardization of the movement operator).
"""
function eqdist(M::MovementModel)
    n = M.Ωsize
    tol = 1e-10
    eq = zeros((0,))
    while true
        λ, ϕ = eigs(M.M; nev = 1, which = :LM, ritzvec = true,
                    tol = tol, maxiter = 1_000, v0 = eq)
        eq = vec(real(ϕ))
        if any(eq .< 0)
            tol *= 1e-1
        else
            break
        end
    end

    eq ./= sum(eq)
    PopState(reshape(eq, n...))
end

"""
    eqdist(M::MovementModel, B0::Real)

Finds the equilibrium spatial population distribution given a movement
operator `M` and total biomass `B0`.
"""
function eqdist(M::MovementModel, B0::Real)
    eq = eqdist(M)
    PopState(B0 .* eq.P)
end

struct MovementRate{F}
    dist::F

    function MovementRate(dist::F) where F<:Function
        new{F}(dist)
    end
end

(mr::MovementRate)(domain::AbstractFisheryDomain) = mr.dist.(domain.distances)

function MovementModel(habitat::Habitat{T, N} where T, habpref::HabitatPreference{F, N}, rate::MovementRate, domain::AbstractFisheryDomain) where {N, F}
    hpref = habpref(habitat; normalize = true)
    dist = rate(domain)

    movop = vec(hpref) .* dist
    movop ./= sum(movop; dims = 1)

    MovementModel(movop, size(domain))
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 30fce022-fdb3-11ea-1059-f97f38c19c0f
using Pkg, DrWatson

# ╔═╡ 30fd1a06-fdb3-11ea-2118-a1facb5194e2
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 3105a0f4-fdb3-11ea-22d8-8b6e68afa993
for i in 5:7
	include(projectdir("models", "05", "m5.$(i)s.jl"))
end

# ╔═╡ 50bbda36-fdb2-11ea-0702-a7ae407a74a7
md"## Clip-05-39s.jl"

# ╔═╡ 30fda02a-fdb3-11ea-2037-6f772dd59330
md"### snippet 5.39"

# ╔═╡ 8f2d61f8-8346-11eb-12a0-c5989da2d4db
md"
!!! note
	In below cell `include()` is used. If this cell needs to be executed again, the notebook needs to be restarted.
"

# ╔═╡ 3108f92a-fdb3-11ea-0959-9b1b6d094a4c
begin
	(s1, p1) = plot_model_coef([m5_5s, m5_6s, m5_7s], [:a, :bN, :bM];
		title="Masked relationships: bN & bM Normal estimates")
	plot(p1;dpi = 460)
end

# ╔═╡ 310f3100-fdb3-11ea-0c35-b3014ae7e73d
s1

# ╔═╡ 3115a60c-fdb3-11ea-16d2-75f8753f7f0f
begin
	(s2, p2) = plot_model_coef([m5_5s, m5_6s, m5_7s], [:a, :bN, :bM];
		title="Masked relationships: bN & bM Normal estimates")
	plot(p2; dpi=460)
 end

# ╔═╡ 31163fc2-fdb3-11ea-03a4-9dc90236e2d6
s2

# ╔═╡ 312018a0-fdb3-11ea-171d-49ed6fe7bbcb
md"## End of clip-05-39s.jl"

# ╔═╡ Cell order:
# ╟─50bbda36-fdb2-11ea-0702-a7ae407a74a7
# ╠═30fce022-fdb3-11ea-1059-f97f38c19c0f
# ╠═30fd1a06-fdb3-11ea-2118-a1facb5194e2
# ╟─30fda02a-fdb3-11ea-2037-6f772dd59330
# ╟─8f2d61f8-8346-11eb-12a0-c5989da2d4db
# ╠═3105a0f4-fdb3-11ea-22d8-8b6e68afa993
# ╠═3108f92a-fdb3-11ea-0959-9b1b6d094a4c
# ╠═310f3100-fdb3-11ea-0c35-b3014ae7e73d
# ╠═3115a60c-fdb3-11ea-16d2-75f8753f7f0f
# ╠═31163fc2-fdb3-11ea-03a4-9dc90236e2d6
# ╟─312018a0-fdb3-11ea-171d-49ed6fe7bbcb


text:
c = 'Q'

c_and = Char(38)

c_at = Int('@')

c_CR = '\n'

str1 = "Hello world!\nJulia is fun\n"

print(str1)

str2 = "I \"have\" \$50, A\\b\n"
print(str2)

str3 = """
    This is some multi-line text.
    Carriage return is inserted for new lines.
    Double quotes can be used as-is: "
    But backslash and dollar still need extra backslashes: \\, \$
    The indentation of the lines is determined by the position of the final triple-quote.
"""
print(str3)

str4a = "Hello"
str4b = "World"
str4ab = string(str4a, " ", str4b, "\n")  # Concatenation
str5ab = str4a * " " * str4b * "\n"       # Same thing

print(str4ab)
print(str5ab)

str6ab = "$str4a $str4b\n"
print(str6ab)

vec = rand(3)
println("Random vector: $vec")
println("sin of 45 degrees = $(sind(45))")

"12345 "^9

str = "abcdefghij"
str[7]

str[7:end]

str[7:7]

sub = SubString(str, 7, 10)

for i = 1:length(str)
    println("Character #$i = '$(str[i])'")
end

str[4] = 'A'    # Error - cannot change strings

function is_palindrome(str)
    return str[end:-1:1] == str
end

strings = ["racecar", "Sit on a Potato Pan, Otis", "sitonapotatopanotis"]
for str in strings
    println("\"$str\": ", is_palindrome(str))
end

function is_palindrome_recursive(str)
    if length(str) ≤ 1
        return true
    elseif str[1] == str[end]
        return is_palindrome_recursive(str[2:end-1])
    else
        return false
    end
end

for str in strings
    println("\"$str\": ", is_palindrome_recursive(str))
end


text:
using Base.Test
using MPFR
using MPC

# real constructors
x = MPCComplex{53,53}()
x = MPCComplex(12)
y = MPCComplex(x)
@test x == y
y = MPCComplex(0xc)
@test x == y
y = MPCComplex(12.)
@test x == y
y = MPCComplex(BigInt(12))
@test x == y
y = MPCComplex(BigFloat(12))
@test x == y
y = MPCComplex(MPFRFloat(12))
@test x == y
y = MPCComplex("12")
@test x == y
y = MPCComplex(float32(12.))
@test x == y
y = MPCComplex(12//1)
@test x == y

# complex constructors
x = MPCComplex(12, 42)
y = MPCComplex(x)
@test x == y
y = MPCComplex(0xc, 0x2a)
@test x == y
y = MPCComplex(12., 42.)
@test x == y
y = MPCComplex(BigInt(12), BigInt(42))
@test x == y
y = MPCComplex(BigFloat(12), BigFloat(42))
@test x == y
y = MPCComplex(MPFRFloat(12), MPFRFloat(42))
@test x == y
y = MPCComplex("(12 42)")
@test x == y
y = MPCComplex(float32(12.), float32(42))
@test x == y
y = MPCComplex(12//1, 42//1)
@test x == y
y = MPCComplex(12 + 42im)
@test x == y

# real/imag
x = MPCComplex(12, 42)
y = MPFRFloat(12)
z = MPFRFloat(42)
@test real(x) == y
@test imag(x) == z
y = MPCComplex(x)
@test real(x) == real(y)
@test imag(x) == imag(y)

# +
x = MPCComplex(12,42)
@test (x + x) == MPCComplex(24, 84)
@test (x + 2) == MPCComplex(14, 42)
@test (x + (2+ 1im)) == MPCComplex(14, 43)
@test (x + 2 + 1im) == MPCComplex(14, 43)

# integer_valued
@test !integer_valued(MPCComplex(2,3))
@test integer_valued(MPCComplex(2))


text:
module AntennaWrapper

using JavaCall

const AbsoluteDate = @jimport org.orekit.time.AbsoluteDate
const Antenna = @jimport org.orekit.gnss.antenna.Antenna
const AntexLoader = @jimport org.orekit.gnss.antenna.AntexLoader
const Class = @jimport java.lang.Class
const DataProvidersManager = @jimport org.orekit.data.DataProvidersManager
const Enum = @jimport java.lang.Enum
const ExtendedPVCoordinatesProvider = @jimport org.orekit.utils.ExtendedPVCoordinatesProvider
const Frame = @jimport org.orekit.frames.Frame
const Frequency = @jimport org.orekit.gnss.Frequency
const FrequencyPattern = @jimport org.orekit.gnss.antenna.FrequencyPattern
const GNSSAttitudeProvider = @jimport org.orekit.gnss.attitude.GNSSAttitudeProvider
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Map = @jimport java.util.Map
const Object = @jimport java.lang.Object
const OneDVariation = @jimport org.orekit.gnss.antenna.OneDVariation
const Optional = @jimport java.util.Optional
const PhaseCenterVariationFunction = @jimport org.orekit.gnss.antenna.PhaseCenterVariationFunction
const ReceiverAntenna = @jimport org.orekit.gnss.antenna.ReceiverAntenna
const SatelliteAntenna = @jimport org.orekit.gnss.antenna.SatelliteAntenna
const SatelliteSystem = @jimport org.orekit.gnss.SatelliteSystem
const SatelliteType = @jimport org.orekit.gnss.antenna.SatelliteType
const TimeScale = @jimport org.orekit.time.TimeScale
const TimeSpanMap = @jimport org.orekit.utils.TimeSpanMap
const TwoDVariation = @jimport org.orekit.gnss.antenna.TwoDVariation
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D

export Antenna
export AntexLoader
export FrequencyPattern
export OneDVariation
export PhaseCenterVariationFunction
export ReceiverAntenna
export SatelliteAntenna
export SatelliteType
export TwoDVariation
export build_attitude_provider
export compare_to
export describe_constable
export equals
export find_satellite_antenna
export get_class
export get_cospar_id
export get_declaring_class
export get_eccentricities
export get_frequencies
export get_name
export get_phase_center_variation
export get_prn_number
export get_receivers_antennas
export get_satellite_code
export get_satellite_system
export get_satellite_type
export get_satellites_antennas
export get_serial_number
export get_sinex_code
export get_type
export get_valid_from
export get_valid_until
export hash_code
export name
export notify
export notify_all
export ordinal
export parse_satellite_type
export to_string
export value
export value_of
export values
export wait

include("antenna.jl")
include("antex_loader.jl")
include("frequency_pattern.jl")
include("one_d_variation.jl")
include("phase_center_variation_function.jl")
include("receiver_antenna.jl")
include("satellite_antenna.jl")
include("satellite_type.jl")
include("two_d_variation.jl")

end


text:
module Deductive

using DataFrames, PrettyTables

export Predicate, truthtable, tableau, prove
export Ē, Ā

include("./utils.jl")
include("./common_algorithms.jl")
include("./expression.jl")
include("./manipulation.jl")
include("./search.jl")

include("./db/index.jl")
include("./db/expressions_index.jl")

include("./proof_utilities.jl")
include("./inference_rules.jl")
include("./transformation_proof.jl")
include("./assertion_proof.jl")
include("./atomic.jl")

include("./set/set.jl")
include("./set/set_operator.jl")
include("./set/set_proof.jl")
include("./set/special_definitions.jl")
include("./set/set_equality.jl")

include("./latex.jl")

# function FreeVariable(sym::Symbol, metadata::Symbol)
#     Sym{Any, Symbol}(sym, metadata)
# end
# function FreeVariable(sym::Symbol)
#     FreeVariable(sym, :free)
# end
FreeVariableType = LogicalSymbol
# LogicalSymbol = Sym{Bool}
# SB = Union{LogicalSymbol, Term}

# Predicate(sym::Symbol) = (st) -> begin
#     Term{Bool}(LogicalSymbol(sym), [st]; metadata=:predicate)
# end

# quantifiers
# todo: fix quantifier logic (see https://github.com/ctrekker/Deductive.jl/issues/2)
Ē = LogicalOperation((n, m) -> nothing, :Ē, 2, false, false)
Ā = LogicalOperation((n, m) -> nothing, :Ā, 2, false, false)
# function Ē(x::FreeVariableType, st::AbstractExpression)
#     quantified_var = Sym{Any, Symbol}(Symbolics.tosymbol(x), :quantified)
#     AbstractExpression(Ē, [quantified_var, replace(st, x => quantified_var)])
# end
# function Ā(x::FreeVariableType, st::AbstractExpression)
#     quantified_var = Sym{Any, Symbol}(Symbolics.tosymbol(x), :quantified)
#     Term(Ā, [quantified_var, substitute(st, x => quantified_var)])
# end

# function substitute_quantified(term::Term, substitution::FreeVariableType)
#     term_args = arguments(term)
#     substitute(term_args[2], term_args[1] => substitution)
# end

# free variable unary "placeholder"
# marker for substitution with a skolem variable
_f = LogicalSymbol(:_f, :definitionallyfree)

function truthtable(st::AbstractExpression)
    vars = sort([variables(st)...])
    combinations = 2 ^ length(vars)
    st_sym = Symbol(string(st))

    # prevent duplicate table keys if the expression is, for example, "a"
    if st_sym == Symbol(first(vars))
        st_sym = Symbol("_" * string(st))
    end

    table = DataFrame(Dict(Symbol(var) => Bool[] for var ∈ vars)..., Dict(st_sym => Bool[])...)

    for i ∈ 0:(combinations-1)
        variable_values = _dict_from_combination_index(vars, i)
        result = evaluate(st, variable_values)
        push!(table, merge(Dict(Symbol(k) => v for (k, v) ∈ variable_values), Dict(st_sym => result)))
    end

    table
end

p, q = LogicalSymbol.([:p, :q])
const demorgan_and = ¬(p ∧ q) => ¬p ∨ ¬q
const demorgan_or = ¬(p ∨ q) => ¬p ∧ ¬q
const double_negative = ¬¬p => p
const material_implication = p → q => ¬p ∨ q
const negated_material_implication = ¬(p → q) => p ∧ ¬q
const material_equivalence = p ⟷ q => (p ∧ q) ∨ (¬p ∧ ¬q)
const negated_material_equivalence = ¬(p ⟷ q) => (p ∧ ¬q) ∨ (¬p ∧ q)

const tableau_replacement_rules = Pair{AbstractExpression, AbstractExpression}[
    demorgan_and,
    demorgan_or,
    double_negative,
    material_implication,
    negated_material_implication,
    material_equivalence,
    negated_material_equivalence
]

# These rules require some extra things, such as substitution by skolem constants after rule application
x, y = LogicalSymbol.([:x, :y])
deny_universal =     ¬Ā(x, y) => ¬y
assert_existential =  Ē(x, y) =>  y
assert_universal =    Ā(x, y) =>  y
deny_existential =   ¬Ē(x, y) => ¬y

tableau(proposition::AbstractExpression; kw...) = tableau([proposition]; kw...)
tableau(propositions...; kw...) = tableau([propositions...]; kw...)
function tableau(propositions::Union{Set, Vector}; skolem_vars=[], proof::Vector{ProofLine}=ProofLine[])
    if length(proof) == 0
        for proposition ∈ propositions
            push!(proof, ProofLine(length(proof) + 1, proposition, "Assumption"))
        end
    end
    simplified_propositions = AbstractExpression[repeated_chain_simplify(p, tableau_replacement_rules) for p ∈ propositions]

    for proposition ∈ simplified_propositions
        if isnothing(find_proof_line_by_statement(proof, proposition))
            push!(proof, ProofLine(length(proof) + 1, proposition, "Replacement Rule <TODO>"))
        end
    end

    return _tableau_simplified(simplified_propositions; skolem_vars=skolem_vars, proof=proof)
end

function _tableau_simplified(propositions::Vector{AbstractExpression}; skolem_vars=[], proof::Vector{ProofLine}=ProofLine[])
    # compile list of all free variables in all propositions
    free_vars = collect(Iterators.flatten([variables(p) for p ∈ propositions]))
    free_vars = filter([istree(v) ? first(arguments(v)) : v for v ∈ free_vars]) do v
        v.metadata == :free
    end

    # when a quantifier is simplified we find existing skolem variables and substitute them in for universal free variables
    # function populate_skolems(term::Term, unary_operator=identity)
    #     placeholder_assertion = unary_operator(substitute_quantified(term, _f))
    #     realized_assertions = [unary_operator(substitute_quantified(term, var)) for var ∈ Iterators.flatten([skolem_vars, free_vars])]

    #     Set([placeholder_assertion, realized_assertions...])
    # end

    # when a quantifier is simplified a new skolem variable can be created in some cases
    # this function creates a new skolem variable and substitutes it into all free variable predicates
    # function create_skolem(term::Term, reduced_propositions, unary_operator=identity)
    #     new_skolem_var = FreeVariable(Symbol("c" * string(length(skolem_vars) + 1)), :skolem)
    #     new_assertion = unary_operator(substitute_quantified(term, new_skolem_var))
    #     realized_placeholders = [substitute(st, _f => new_skolem_var) for st ∈ reduced_propositions]

    #     new_skolem_var, Set([realized_placeholders..., new_assertion])
    # end

    ordered_propositions = [propositions...]
    for p ∈ propositions
        pv = p

        if !istree(pv) || metadata(pv) == :predicate
            # check for contradiction
            contradiction_index = findfirst([isequal(¬p, q) for q ∈ ordered_propositions])
            if !isnothing(contradiction_index)
                contradictory_statement = ordered_propositions[contradiction_index]
                contradiction_references = [
                    find_proof_line_by_statement(proof, pv),
                    find_proof_line_by_statement(proof, contradictory_statement)
                ]

                if any(isnothing.(contradiction_references))
                    @error "Contradiction found, but one or more references is missing"
                else
                    push!(proof, ProofLine(length(proof) + 1, p ∧ contradictory_statement, "Contradiction", contradiction_references))
                end

                return false
            end
        else
            term_op = operation(pv)
            term_args = arguments(pv)
            reduced_propositions = setdiff(propositions, Set([p]))

            if length(term_args) == 2 # binary operation
                line_ref = find_proof_line_by_statement(proof, pv)

                if term_op == ∧
                    # break up term and make one recursive call

                    # add simplification to proof
                    push!(proof, ProofLine(length(proof) + 1, term_args[1], "Simplification", line_ref))
                    push!(proof, ProofLine(length(proof) + 1, term_args[2], "Simplification", line_ref))

                    if !tableau(reduced_propositions ∪ Set(term_args); skolem_vars=skolem_vars, proof=proof)
                        return false
                    end
                    break
                elseif term_op == ∨
                    # consider both arguments of term and make two recursive calls

                    # case-by-case checking
                    case1_linenum, case2_linenum = length(proof) + 1, length(proof) + 2
                    push!(proof, ProofLine(case1_linenum, term_args[1], "Case 1", line_ref))
                    push!(proof, ProofLine(case2_linenum, term_args[2], "Case 2", line_ref))

                    if !tableau(reduced_propositions ∪ Set([term_args[1]]); skolem_vars=skolem_vars, proof=proof) && !tableau(reduced_propositions ∪ Set([term_args[2]]); skolem_vars=skolem_vars, proof=proof)
                        push!(proof, ProofLine(length(proof) + 1, ¬pv ∧ pv, "Contradiction", [proof[case1_linenum], proof[case2_linenum]]))

                        return false
                    end
                    break
                elseif term_op == Ā  # asserted universal quantifier
                    # create free variable
                    if !tableau(reduced_propositions ∪ populate_skolems(pv); skolem_vars=skolem_vars, proof=proof)
                        return false
                    end
                    break
                elseif term_op == Ē  # asserted existential quantifier
                    # substitute skolem variable
                    new_skolem_var, new_propositions = create_skolem(pv, reduced_propositions)
                    if !tableau(reduced_propositions ∪ new_propositions; skolem_vars=[skolem_vars..., new_skolem_var], proof=proof)
                        return false
                    end
                    break
                end
            else # unary operation (¬, logical statements)
                subterm = first(term_args)

                if istree(subterm)
                    subterm_op = operation(subterm)

                    # check for negated quantifiers
                    if subterm_op == Ā  # denied universal quantifier
                        new_skolem_var, new_propositions = create_skolem(subterm, reduced_propositions, ¬)
                        if !tableau(reduced_propositions ∪ new_propositions; skolem_vars=[skolem_vars..., new_skolem_var], proof=proof)
                            return false
                        end
                        break
                    elseif subterm_op == Ē  # denied existential quantifier
                        if !tableau(reduced_propositions ∪ populate_skolems(subterm, ¬); skolem_vars=skolem_vars, proof=proof)
                            return false
                        end
                        break
                    end
                end
            end
        end
    end

    return true
end

prove(proposition::AbstractExpression; kw...) = prove([proposition]; kw...)
prove(propositions...; kw...) = prove([propositions...]; kw...)
function prove(propositions::Union{Set, Vector})
    proof = ProofLine[]
    tableau(propositions; proof=proof)

    proof
end

_dict_from_combination_index(variables, index::Int) = Dict([variables[i] => Bool(index >> (i-1) & 1) for i ∈ 1:length(variables)])

end # module


text:
#include("test_lists.jl")
include("test_trees.jl")
include("test_lvltree.jl")
include("test_blocktree.jl")


text:
"""
Coordinate transform function `ctparent` convert the coordinates of `point = (x1, y1)` on a child frame `c1` to its corresponding coordinates `(x0, y0)` on a parent frame `c0`.

For example,
saying you have a circle centered at `(x1, y1)` of the coordinate system `c1`, where `c1` is located on `contextdim` of `c0`:
```julia
c0 = context();
c1 = context(contextdim...;units=UnitBox(0,0,2,2));
circle1 = circle(x1,y1,0.1);
compose(c0, (c1, circle1));
```
Using the transformed coordinates `(x0,y0)=ctparent((x1,y1),contextdim;UnitBox=(0,0,2,2))` we can draw the same circle as `circle1` by `compose(c0, circle(x0,y0,0.1))`.

Broadcasting example:
```julia
ctparent.([(x1,y1), (x2,y2)], [contextdim]; UnitBox=(0,0,2,2))
```
"""
function ctparent(point::NTuple{2,<:Real},contextdim::NTuple{4,<:Real};UnitBox::NTuple{4,<:Real}=(0,0,1,1))
    x1, y1 = point;
    xstart1a, ystart1a, width1a, height1a = contextdim;
    xstart1b, ystart1b, width1b, height1b = UnitBox;
    Δx1 = x1-xstart1b;
    Δy1 = y1-ystart1b;
    cfactorx = width1a/1; # default context()'s width/height is 1
    cfactory = height1a/1;
    ubfactorx = 1/width1b; # default UnitBox()'s width/height is 1
    ubfactory = 1/height1b;
    x0 = xstart1a + Δx1*cfactorx*ubfactorx;
    y0 = ystart1a + Δy1*cfactory*ubfactory;
    return (x0, y0)
end

"""
`ctchildwidth(width0, contextdim::NTuple{4,<:Real};UnitBox::NTuple{4,<:Real}=(0,0,1,1))` convert the width in `c0` to the length it should be in `c1` that the line:
```
c0 = context();
c1 = context(contextdim...;units=UnitBox(0,0,2,2));
width1=ctchildwidth(width0, contextdim;UnitBox=(0,0,2,2));
compose(c0, (c1, line([(0,0),(width1,0)])) )
```
has the length looks visually the same as:
```
compose(c0, line([(0,0),(width0,0)]) )
```
"""
function ctchildwidth(width0, contextdim::NTuple{4,<:Real};UnitBox::NTuple{4,<:Real}=(0,0,1,1))
	# width0/contextdim[3] = width1/UnitBox[3]
	return width0*UnitBox[3]/contextdim[3] # width1
end

"""
`height1 = ctchildheight(height0, contextdim::NTuple{4,<:Real};UnitBox::NTuple{4,<:Real}=(0,0,1,1))` works the same as `ctchildwidth()` but rescale in y direction. See `ctchildwidth`'s documentation.
"""
function ctchildheight(height0, contextdim::NTuple{4,<:Real};UnitBox::NTuple{4,<:Real}=(0,0,1,1))
	# height0/contextdim[4] = height1/UnitBox[4]
	return height0*UnitBox[4]/contextdim[4] # height1
end

dfthw = 0.015; # default head width
dfthh = 0.025; # default head height
"""
`myarrow(startpoint, arrangle, arrlength; headwidth=$dfthw, headheight=$dfthh) |> polygon`
"""
function myarrow(startpoint, arrangle, arrlength; headwidth=dfthw, headheight=dfthh)
	rot(θ,point) = (point[1]*cos(θ) - point[2]*sin(θ), point[1]*sin(θ) + point[2]*cos(θ));
	xo,yo = (0.0,0.0);
	head_x =xo + arrlength;
	headpoint = (head_x, yo);
	sidepoint_a = headpoint.+(-headheight,yo-headwidth/2);
	sidepoint_b = headpoint.+(-headheight,yo+headwidth/2);
	function lim_x(sidepoint_a)
		sidepoint_a = maximum.(([sidepoint_a[1], xo], [sidepoint_a[2]]))
	end
	sidepoint_a = lim_x(sidepoint_a);
	sidepoint_b = lim_x(sidepoint_b);

	eq_triangle = [sidepoint_a,
				   sidepoint_b,
				   headpoint]

	tailwidth = headwidth/5;
	tail_a = (xo, yo+tailwidth);
	tail_b = (xo, yo-tailwidth);
	tail_c = (sidepoint_b[1], yo-tailwidth); # sidepoint_a[1] should equal sidepoint_b[1]
	tail_d = (sidepoint_a[1], yo+tailwidth);

	pts = rot.([arrangle], [tail_a, tail_b, tail_c, sidepoint_a, headpoint, sidepoint_b, tail_d]);
	xs = [pt[1] for pt in pts] .+ startpoint[1];
	ys = [pt[2] for pt in pts] .+ startpoint[2];

	return tuple.(xs, ys)
end

"""
`myarrow(startpoint, arrangle, arrlength, contextdim::NTuple{4,<:Real}; UnitBox::NTuple{4,<:Real}=(0,0,1,1))` is the same as `myarrow(startpoint, arrangle, arrlength; headwidth=$dfthw, headheight=$dfthh)` but the arrows `headwidth` and `headheight` will be converted to make sure their size is consistent.
"""
function myarrow(startpoint, arrangle, arrlength, contextdim::NTuple{4,<:Real}; UnitBox::NTuple{4,<:Real}=(0,0,1,1), headwidth=dfthw, headheight=dfthh)
	hh = ctchildwidth(headheight, contextdim; UnitBox=UnitBox); # headheight is "along x axis" at zero degree, so we have to use `ctchildwidth` to covert.
	hw = ctchildheight(headwidth, contextdim; UnitBox=UnitBox);
	return myarrow(startpoint, arrangle, arrlength; headwidth=hw, headheight=hh)
end

"""
`composelabels(xstr, ystr, contextdim)` add x and y labels on the parent frame of `contextdim`.

```julia
composelabels(labels::Vector{<:AbstractString}, contextdim::NTuple{4,<:Real}; labelloc::Vector{Symbol}=[:bottom,:left], labeldir::Vector{Symbol}=[:h,:h], labelcolor=("black","black"))
```
- supported `labelloc`: `:bottom`,`:top`,`:left`,`:right`.
- `labeldir=:h` means the label is horizontally placed; `labeldir=:h`, vertically placed.

# Example:
```julia
contextdim = (0.1,0.4,0.8,0.55); # the position of `cplot1`
cplot1 = (context(contextdim..., units=UnitBox(UB1...)),
blablabla...); # a plot perhaps

cplot1withlabels = (context(), composelabels(["time","velocity"],contextdim;labelloc=[:left,:top], labeldir=[:v,:h])..., cplot1);
compose(context(), cplot1withlabels)
```
"""
function composelabels(labels::Vector{<:AbstractString}, contextdim::NTuple{4,<:Real}; labelloc::Vector{Symbol}=[:bottom,:left], labeldir::Vector{Symbol}=[:h,:h], labelcolor=("black","black"))
    labh = 0.03; # labelheight
    textpos = (0.5,0.5); # of created label context
    (pad1x, pad1y, width1, height1) = contextdim;

    lenl = length(labels);

    if length(labelloc) != lenl || length(labeldir) != lenl
        error("`labelloc`/`labeldir` has to be the same length as `labels` for specifying labels' location/orientation.");
    end

    contexts = [];
    for i = 1:lenl
        textrot = ();
        textalign = ();
        # txtpt = (1.0,0.5);
        # (pad1x, pad1y, width1, height1) = contextdim;
        # labh = 0.03;
        #
        # leftlabel = (a1,pad1y,pad1x-a1,height1);
        # cplot1a = (context(), (context(leftlabel...),Compose.text(txtpt...,"hello world",hright,vcenter,Rotation(-1/2*pi,0.5,0.5)) ,(context(), rectangle(), stroke("pink"),fillopacity(0.3)))
        # ,cplot1);

        if labelloc[i] == :left
            a1 = maximum([0.0, pad1x-labh]);
            contextdim_txt = (a1,pad1y,pad1x-a1,height1);
            textpos = (0.9,0.5);
            textalign= (hright, vcenter);
        end

        if labelloc[i] == :bottom
            a1 = minimum([1.0-pad1y-height1, labh]);
            contextdim_txt = (pad1x, pad1y+height1, width1, a1);
            textpos = (0.5, 0.1);
            textalign= (hcenter,vtop);
        end

        if labelloc[i] == :top
            a1 = maximum([0.0, pad1y-labh]);
            contextdim_txt = (pad1x, a1, width1, pad1y-a1);
            textpos = (0.5,0.9);
            textalign= (hcenter,vbottom);
        end

        if labelloc[i] == :right
            a1 = minimum([1.0-pad1x-width, labh]);
            contextdim_txt = (pad1x+width1, pad1y, a1, height1);
            textpos = (0.1, 0.5);
            textalign= (hleft,vcenter);
        end

        if labeldir[i] == :v
            textrot = (Rotation(-pi/2,0.5,0.5),);
            textpos = (0.5, 0.5);
            textalign= (hcenter,vcenter);
        end
        # convert dictionary to named tuple
        # TODO: write this technique to Notes
        # opts = NamedTuple{Tuple(keys(contextopt_d))}(values(contextopt_d));

        push!(contexts,
              (context(contextdim_txt...),
               Compose.text(textpos..., labels[i], textalign..., textrot...),fill(labelcolor[i])),
        );

    end
    return contexts
end

"""
`pos2points(xstart, ystart, width, height)` convert the position to a vector of points `[(xstart, ystart), (xstart, ystart+height), (xstart+width, ystart+height), (xstart+width,ystart)]`.
"""
function pos2points(xstart, ystart, width, height)
    return [(xstart, ystart), (xstart, ystart+height), (xstart+width, ystart+height), (xstart+width,ystart)]
end

function pos2points(pos::NTuple{4,<:Real})
    pos2points(pos...);
end

"""
`pointshift(point::NTuple{2,<:Real}, shift2::NTuple{2,<:Real})` returns `(point[1]+shift2[1], point[2]+shift2[2])`.
"""
function pointshift(point::NTuple{2,<:Real}, shift2::NTuple{2,<:Real})
    return (point[1]+shift2[1], point[2]+shift2[2])
end


text:
# Script

const SIGHASH_ALL = UInt32(1)
const SIGHASH_NONE = UInt32(2)
const SIGHASH_SINGLE = UInt32(3)

const SCRIPT_TYPES = Dict([
    ("P2PKH", [0x6f, 0x00]),
    ("P2SH", [0xc4, 0x05])
])

# Block

GENESIS_BLOCK_HASH = Dict([
    (false, hex2bytes("000000000019d6689c085ae165831e934ff763ae46a2a6c172b3f1b60a8ce26f")),
    (true, hex2bytes("000000000933ea01ad0ee984209779baaec3ced90fa3f408719526f8d77f4943"))
])

# Network

USER_AGENT = "/bitcoin.jl:0.1/"

const DEFAULT = Dict([
    ("version", 70015),
    ("services", zero(UInt64)),
    ("ip", ip"0.0.0.0"),
    ("port", Dict([(false, 8333),
                   (true, 18333)])),
    ("rpcport", Dict([(false, 8332),
                   (true, 18332)])),
    ("start_height", zero(UInt32)),
    ("relay", true)
])

const IPV4_PREFIX = append!(fill(0x00, 10), [0xff, 0xff])

const NODE_URL = "btc.brane.cc"

const TX_DATA_TYPE = 1
const BLOCK_DATA_TYPE = 2
const FILTERED_BLOCK_DATA_TYPE = 3
const COMPACT_BLOCK_DATA_TYPE = 4

const DATA_MESSAGE_NAME = Dict([
    (1, "MSG_TX"),
    (2, "MSG_BLOCK"),
    (3, "MSG_FILTERED_BLOCK"),
    (4, "MSG_CMPCT_BLOCK"),
    (5, "MSG_WITNESS_BLOCK"),
    (6, "MSG_WITNESS_TX"),
    (7, "MSG_FILTERED_WITNESS_BLOCK")])

const DATA_MESSAGE_TYPE = Dict([
    ("MSG_TX", 1),
    ("MSG_BLOCK", 2),
    ("MSG_FILTERED_BLOCK", 3),
    ("MSG_CMPCT_BLOCK", 4),
    ("MSG_WITNESS_BLOCK", 5),
    ("MSG_WITNESS_TX", 6),
    ("MSG_FILTERED_WITNESS_BLOCK", 7)])

"""
NETWORK_MAGIC is testnet if `true`
"""
const NETWORK_MAGIC = Dict([
    (false, 0xd9b4bef9)
    (true, 0x0709110b)
])

SERVICES_NAME = Dict([
    (0x0000000000000001, "NODE_NETWORK"),
    (0x0000000000000002, "NODE_GETUTXO"),
    (0x0000000000000004, "NODE_BLOOM"),
    (0x0000000000000008, "NODE_WITNESS"),
    (0x0000000000000400, "NODE_NETWORK_LIMITED")
])


text:
ENV["PATH"] = "/usr/local/opt/llvm/bin:$(ENV["PATH"])"

using Clang

cd(Pkg.dir("LibPQ", "src", "headers"))

context = wrap_c.init(;
    clang_includes=["/usr/local/opt/libpq/include"],
    common_file="libpq_common.jl",
)

context.headers = [
    "/usr/local/opt/libpq/include/libpq-fe.h",
]

run(context)


text:
function _withtemp(f, file)
    try f(file)
    catch err
        rethrow()
    finally
        try rm(file; force = true)
        catch
        end
    end
end

# Runs a function at a commit on a repo and afterwards goes back
# to the original commit / branch.
function _withcommit(f, repo, commit)
    original_commit = _shastring(repo, "HEAD")
    LibGit2.transact(repo) do r
        branch = try LibGit2.branch(r) catch err; nothing end
        try
            LibGit2.checkout!(r, _shastring(r, commit))
            f()
        catch err
            rethrow(err)
        finally
            if branch !== nothing
                LibGit2.branch!(r, branch)
            else
                LibGit2.checkout!(r, original_commit)
            end
        end
    end
end

_shastring(r::LibGit2.GitRepo, targetname) = string(LibGit2.revparseid(r, targetname))
_shastring(dir::AbstractString, targetname) = LibGit2.with(r -> _shastring(r, targetname), LibGit2.GitRepo(dir))

_benchinfo(str) = printstyled(stdout, "PkgBenchmark: ", str, "\n"; color = Base.info_color())
_benchwarn(str) = printstyled(stdout, "PkgBenchmark: ", str, "\n"; color = Base.info_color())

############
# Markdown #
############

_idrepr(id) = (str = repr(id); str[coalesce(findfirst(isequal('['), str), 0):end])
_intpercent(p) = string(ceil(Int, p * 100), "%")
_resultrow(ids, t::BenchmarkTools.Trial, col_widths) =
    _resultrow(ids, minimum(t), col_widths)

_update_col_widths!(col_widths, ids, t::BenchmarkTools.Trial) =
    _update_col_widths!(col_widths, ids, minimum(t))

function _resultrow(ids, t::BenchmarkTools.TrialEstimate, col_widths)
    t_tol = _intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = _intpercent(BenchmarkTools.params(t).memory_tolerance)
    timestr = BenchmarkTools.time(t) == 0 ? "" : string(BenchmarkTools.prettytime(BenchmarkTools.time(t)), " (", t_tol, ")")
    memstr = BenchmarkTools.memory(t) == 0 ? "" : string(BenchmarkTools.prettymemory(BenchmarkTools.memory(t)), " (", m_tol, ")")
    gcstr = BenchmarkTools.gctime(t) == 0 ? "" : BenchmarkTools.prettytime(BenchmarkTools.gctime(t))
    allocstr = BenchmarkTools.allocs(t) == 0 ? "" : string(BenchmarkTools.allocs(t))
    return "| $(rpad("`"*_idrepr(ids)*"`", col_widths[1])) | $(lpad(timestr, col_widths[2])) | $(lpad(gcstr, col_widths[3])) | $(lpad(memstr, col_widths[4])) | $(lpad(allocstr, col_widths[5])) |"
end

function _update_col_widths!(col_widths, ids, t::BenchmarkTools.TrialEstimate)
    t_tol = _intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = _intpercent(BenchmarkTools.params(t).memory_tolerance)
    timestr = BenchmarkTools.time(t) == 0 ? "" : string(BenchmarkTools.prettytime(BenchmarkTools.time(t)), " (", t_tol, ")")
    memstr = BenchmarkTools.memory(t) == 0 ? "" : string(BenchmarkTools.prettymemory(BenchmarkTools.memory(t)), " (", m_tol, ")")
    gcstr = BenchmarkTools.gctime(t) == 0 ? "" : BenchmarkTools.prettytime(BenchmarkTools.gctime(t))
    allocstr = BenchmarkTools.allocs(t) == 0 ? "" : string(BenchmarkTools.allocs(t))
    idrepr = "`"*_idrepr(ids)*"`"
    for (i, s) in enumerate((idrepr, timestr, gcstr, memstr, allocstr))
        w = length(s)
        if (w > col_widths[i]) col_widths[i] = w end
    end
end

function _resultrow(ids, t::BenchmarkTools.TrialJudgement, col_widths)
    t_tol = _intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = _intpercent(BenchmarkTools.params(t).memory_tolerance)
    t_ratio = @sprintf("%.2f", BenchmarkTools.time(BenchmarkTools.ratio(t)))
    m_ratio =  @sprintf("%.2f", BenchmarkTools.memory(BenchmarkTools.ratio(t)))
    t_mark = _resultmark(BenchmarkTools.time(t))
    m_mark = _resultmark(BenchmarkTools.memory(t))
    timestr = "$(t_ratio) ($(t_tol)) $(t_mark)"
    memstr = "$(m_ratio) ($(m_tol)) $(m_mark)"
    return "| $(rpad("`"*_idrepr(ids)*"`", col_widths[1])) | $(lpad(timestr, col_widths[2])) | $(lpad(memstr, col_widths[3])) |"
end

function _update_col_widths!(col_widths, ids, t::BenchmarkTools.TrialJudgement)
    t_tol = _intpercent(BenchmarkTools.params(t).time_tolerance)
    m_tol = _intpercent(BenchmarkTools.params(t).memory_tolerance)
    t_ratio = @sprintf("%.2f", BenchmarkTools.time(BenchmarkTools.ratio(t)))
    m_ratio =  @sprintf("%.2f", BenchmarkTools.memory(BenchmarkTools.ratio(t)))
    t_mark = _resultmark(BenchmarkTools.time(t))
    m_mark = _resultmark(BenchmarkTools.memory(t))
    timestr = "$(t_ratio) ($(t_tol)) $(t_mark)"
    memstr = "$(m_ratio) ($(m_tol)) $(m_mark)"
    idrepr = "`"*_idrepr(ids)*"`"
    for (i, s) in enumerate((idrepr, timestr, memstr))
        w = length(s)
        if (w > col_widths[i]) col_widths[i] = w end
    end
end

_resultmark(sym::Symbol) = sym == :regression ? _REGRESS_MARK : (sym == :improvement ? _IMPROVE_MARK : "")

const _REGRESS_MARK = ":x:"
const _IMPROVE_MARK = ":white_check_mark:"


text:
# Stuff in here should probably move to NaiveNASflux once mature enough...

struct SizePseudoTransparent{T} <: NaiveNASlib.DecoratingTrait
    base::T
end
NaiveNASlib.base(t::SizePseudoTransparent) = t.base

NaiveNASlib.all_in_Δsize_graph(mode, ::SizePseudoTransparent, args...) = NaiveNASlib.all_in_Δsize_graph(mode, SizeInvariant(), args...)

NaiveNASlib.nout(::SizePseudoTransparent, v::AbstractVertex) = calc_outsize(v)
NaiveNASlib.nin(::SizePseudoTransparent, v::AbstractVertex) = nin(SizeInvariant(), v)

calc_outsize(v::AbstractVertex, vo::AbstractVertex=v) = calc_outsize(base(v), vo)
calc_outsize(v::CompVertex, vo::AbstractVertex) = calc_outsize(v.computation, vo)

"""
    reshape_keepshape(x, shape)

Same as `reshape` except that `shape[i] = 0` results in `size(y, i+offs) == size(x, i)` where `y` is the output of this function and `offs = ndims(x) - length(shape)`.

This is basically trying to be compliant to the ONNX Reshape operator although the description there of how to interpret a shape of `0` is a bit vague.
"""
function reshape_keepshape(x, shape)
    offs = ndims(x) - length(shape)
    newshape = Flux.Zygote.@ignore map(enumerate(shape)) do (ind, new)
        new == 0 && return size(x, ind+offs)
        return new
    end
    return reshape(x, newshape...)
end

"""
    Reshape{T}
    Reshape(dims...; activation_dim=actdim(guess_layertype(length(dims))))
    Reshape(dims; activation_dim=actdim(guess_layertype(length(dims))))

Reshape operation wrapped in a struct for the sole purpose of handling size changes in neighbouring vertices.

Ensures that changes in output size of the input vertex are possible to reshape into the shape given by `dims`.

Will treat `Integer`s in dims as fixed, meaning that if `dims[activation_dim] isa Integer` the input size of the output vertices will be fixed as well.
"""
mutable struct Reshape{T}
    adim::Int
    dims::T
end
Reshape(dims...; activation_dim=actdim(guess_layertype(length(dims)))) = Reshape(activation_dim, dims)
Reshape(dims; activation_dim=actdim(guess_layertype(length(dims)))) = Reshape(activation_dim, dims)

(r::Reshape)(x) = any(==(0), r.dims) ? reshape_keepshape(x, r.dims) : reshape(x, r.dims)

function calc_outsize(r::Reshape, v)
    invertex = first(inputs(v))
    outshape = r.dims[r.adim]
    outshape == 0 && return nout(invertex)
    outshape isa Integer && return outshape
    outshape isa Colon && return 0 # Must be set later...
    error("Unknown outshape: " + outshape)
end

function NaiveNASlib.Δsize!(r::Reshape{<:Tuple}, ins::AbstractVector, outs::AbstractVector; kwargs...)
    r.dims = Tuple(map(enumerate(r.dims)) do (i, s)
        s isa Colon && return s
        i == r.adim && return length(outs)
        return s
    end)
end

NaiveNASflux.layer(r::Reshape) = r
NaiveNASflux.actdim(r::Reshape) = r.adim
NaiveNASflux.actrank(r::Reshape) = length(r.dims)

# If Reshape is wrapped in an AbstractMutableComp we will hit this method instead due to how NaiveNASflux unwraps things
NaiveNASlib.compconstraint!(case, s::NaiveNASlib.AbstractJuMPΔSizeStrategy, ::Type{<:Reshape}, data) = NaiveNASlib.compconstraint!(case, s, layer(data.vertex), data)

# Special case: Reshape to 2D with variable batch size. Maybe a more general case when this is ok is hiding here somewhere...
NaiveNASlib.compconstraint!(case, s::NaiveNASlib.AbstractJuMPΔSizeStrategy, r::Reshape{Tuple{Int, Colon}}, data) = reshape_nout_constraint(s, Colon(), r, filter(vin -> vin in keys(data.noutdict), inputs(data.vertex)), data)

function NaiveNASlib.compconstraint!(case, s::NaiveNASlib.AbstractJuMPΔSizeStrategy, r::Reshape, data)
    ins = filter(vin -> vin in keys(data.noutdict), inputs(data.vertex))

    reshape_nout_constraint(s, r.dims[r.adim], r, ins, data)
    isempty(ins) && return

    fixeddims = filter(dim -> dim isa Integer && dim != 0, collect(r.dims))
    isempty(fixeddims) && return

    if length(fixeddims) == length(r.dims)
        # No size change possible!
        @constraint(data.model, [i=1:length(ins), j=1:length(fixeddims)], data.noutdict[ins[i]] == nout(ins[i]))
        return
    end

    # Make sure that nout of inputs is an integer multiple of at least one of the output sizes
    # It is unfortunately too restrictive to just enforce fv_dims[j] * fixeddims[j] == nout[ins[i]] ∀ i, j
    # as it is enough that this is fulfilled for one j

    # Unfortunately, something like  fv_dims[j] * fixeddims[j] == nout[ins[i]] * b[i] where b is a binary variable is not possible as it is nonlinear.
    # Instead we use i slack(ish) variables with lower bound 0 and force the smallest of them to be 0 (i.e at least one dimension has zero slack).
    fv_dims = @variable(data.model, [1:length(fixeddims)], integer=true)
    slack = @variable(data.model, [1:length(fixeddims)], integer=true, lower_bound=0)
    @constraint(data.model, [i=1:length(ins), j=1:length(fixeddims)], fv_dims[j] * fixeddims[j] ==  data.noutdict[ins[i]] + slack[j])

    # Force the smallest slack to be 0 through a big-M strategy
    atleast_one = @variable(data.model, [1:length(fixeddims)], binary=true)
    min_slack = @variable(data.model, integer=true)
    M = 10000
    @constraint(data.model,[i=1:length(fixeddims)], -M * (1 - atleast_one[i]) <= min_slack - slack[i])
    @constraint(data.model,[i=1:length(fixeddims)], min_slack - slack[i] <= M * (1-atleast_one[i]))
    @constraint(data.model, sum(atleast_one) == 1)
    @constraint(data.model, min_slack == 0)
end

# Case 1: Output size is fixed so we will only put constraints on the inputs to ensure reshaping is possible
function reshape_nout_constraint(s, outsize::Integer, r, ins, data)
    if outsize == 0

        fix_size = true
        for iv in ins
            inrank = unique(actrank(iv))[] # Should be one element or else we are fked
            inadim = unique(actdim(iv))[] # Should be one element or else we are fked

            if inrank - length(r.dims) + r.adim + 1== inadim
                # Dimension to keep size of happens to be input layers activation dimension.
                # This is basically a SizeInvariant vertex
                @constraint(data.model, [i=1:length(ins)], data.noutdict[data.vertex] == data.noutdict[iv])
                fix_size = false
            end
        end

        if fix_size
            @constraint(data.model, [i=1:length(ins)], data.noutdict[data.vertex] == nout(data.vertex))
        end

    else
        @constraint(data.model, data.noutdict[data.vertex] == outsize)
    end
end

# Case 2: Outsize is not fixed so we need to change it so that the output vertices change their input size
function reshape_nout_constraint(s, outsize::Colon, r, ins, data)
    @constraint(data.model, [i=1:length(ins)], data.noutdict[data.vertex] / nout(data.vertex) == data.noutdict[ins[i]] / nout(ins[i]))
end

"""
    Flatten
    Flatten(dim)

Flattens input array along `dim` to a 2D array.

If input has shape `(d_1, d_2, ... d_n)` then the output will have shape `(d_1 X d_2 ... d_(dim-1), d_dim X d_(dim+1) ... X dn)`.

If `dim = 0`, the shape of the output is `(d_1 X d_2 ... d_n, 1)`, where the shape of the input is `(d_0, d_1, ... d_n)`.
"""
struct Flatten
    dim::Int
end

(f::Flatten)(x) = flatten(x, f.dim)
function flatten(x, dim)
    dim == 0 && return reshape(x, :, 1)
    xs = size(x)
    absdim = dim < 0 ? length(xs) + dim : dim
    return reshape(x, prod(xs[1:absdim]), prod(xs[absdim+1:end]))
end

NaiveNASflux.layer(f::Flatten) = f
NaiveNASflux.actdim(f::Flatten) = 1
NaiveNASflux.actrank(f::Flatten) = 1

calc_outsize(::Flatten, v) = 0

# If Flatten is wrapped in an AbstractMutableComp we will hit this method instead due to how NaiveNASflux unwraps things
NaiveNASlib.compconstraint!(case, s::NaiveNASlib.AbstractJuMPΔSizeStrategy, ::Type{<:Flatten}, data) = NaiveNASlib.compconstraint!(case, s, layer(data.vertex), data)

function NaiveNASlib.compconstraint!(case, ::AbstractJuMPΔSizeStrategy, f::Flatten, data)
    ins = filter(vin -> vin in keys(data.noutdict), inputs(data.vertex))
    for iv in ins
        inadim = unique(actdim(iv))[] # Should be one element or else we are fked

        absdim = f.dim > 0 ? f.dim : f.dim + unique(actrank(iv))[] + 1 # Should be one element or else we are fked

        if f.dim == 0 || inadim <= absdim
            # Size of input affect output size
            @constraint(data.model, [i=1:length(ins)], data.noutdict[data.vertex] / nout(data.vertex) ==  data.noutdict[iv] / nout(iv))
        end
        # else: Size of input affects batch size => no constraint needed
    end
end


text:
# Object-based wrapper for level-1 BLIS fused-vector routines.
#

@blis_ccall_object(axpy2v,
                   Cvoid,
                   αx, xObj,
                   αy, xObj,
                   x,  xObj,
                   y,  xObj,
                   z,  xObj)
@blis_ccall_object(axpyf,
                   Cvoid,
                   α, xObj,
                   a, xObj,
                   x, xObj,
                   y, xObj)
@blis_ccall_object(dotxf,
                   Cvoid,
                   α, xObj,
                   a, xObj,
                   x, xObj,
                   β, xObj,
                   y, xObj)
# TODO: Document of dotaxpyv and dotxaxpyf is not consistent
#       with header files. Prospone.


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

name = "mpv"
version = v"0.32.0"

# Collection of sources required to complete build #
sources = [
    ArchiveSource("https://github.com/mpv-player/mpv/archive/v0.32.0.tar.gz", "9163f64832226d22e24bbc4874ebd6ac02372cd717bef15c28a0aa858c5fe592")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd mpv-*
ln -s /usr/bin/pkg-config /usr/bin/$target-pkg-config

# libstc++ required for c++ libs when using C compiler
if [[ "${nbits}" == 32 ]]; then
    export LDFLAGS="-L${prefix}/lib -liconv -Wl,-rpath-link,/opt/${target}/${target}/lib"
elif [[ "${target}" != *-apple-* ]]; then
    export LDFLAGS="-L${prefix}/lib -liconv -Wl,-rpath-link,/opt/${target}/${target}/lib64"
else
    export LDFLAGS="-L${prefix}/lib -liconv"
fi

# pkg-config files for ffmpeg are in $bindir for windows
if [[ "${target}" == *-mingw* ]]; then
    cp $bindir/pkgconfig/lib*  /workspace/destdir/lib/pkgconfig/
fi
python3 bootstrap.py

# No opengl on MacOS
if [[  "${target}" == *-apple-* ]]; then
    TARGET=$target ./waf --prefix=${prefix} --disable-manpage-build --enable-sdl2 --disable-gl configure
else
    TARGET=$target ./waf --prefix=${prefix} --disable-manpage-build --enable-sdl2 configure
fi
./waf build -j${nproc}
./waf install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
filter!(!Sys.isfreebsd, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("mpv", :mpv)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Libiconv_jll", uuid="94ce4f54-9a6c-5748-9c1c-f9c7231a4531"))
    Dependency(PackageSpec(name="SDL2_jll", uuid="ab825dc5-c88e-5901-9575-1e5e20358fcf"))
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="FFMPEG_jll", uuid="b22a6f82-2f65-5046-a5b2-351ab43fb4e5"))
    Dependency(PackageSpec(name="Lua_jll", uuid="a4086b1d-a96a-5d6b-8e4f-2030e6f25ba6"); compat="~5.3.6")
    Dependency(PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8"))
    Dependency(PackageSpec(name="Xorg_libXrandr_jll", uuid="ec84b674-ba8e-5d96-8ba1-2a689ba10484"))
    Dependency(PackageSpec(name="Xorg_libXinerama_jll", uuid="d1454406-59df-5ea1-beac-c340f2130bc3"))
    Dependency(PackageSpec(name="Libglvnd_jll", uuid="7e76a0d4-f3c7-5321-8279-8d96eeed0f29"))
    Dependency(PackageSpec(name="Xorg_libX11_jll", uuid="4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"))
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
    Conjunction{N, VT<:AbstractVector{<:Tuple{<:Predicate, <:AbstractVector{Int}}}} <: Predicate{N}

A conjunction of predicates of arity `N`.

### Fields

- `conjuncts` -- vector of conjuncts
"""
struct Conjunction{N, VT<:AbstractVector{<:Tuple{<:Predicate, <:AbstractVector{Int}}}} <: Predicate{N}
    conjuncts::VT

    function Conjunction(conjuncts::VT; N::Int=maximum([length(v) for (p, v) in conjuncts])
                        ) where {VT<:AbstractVector{<:Tuple{<:Predicate, <:AbstractVector{Int}}}}
        return new{Val{N}, VT}(conjuncts)
    end
end

# constructor for unary conjunction with only unary conjuncts
function Conjunction(conjuncts::AbstractVector{<:Predicate{Val{1}}})
    return Conjunction([(c, [1]) for c in conjuncts])
end

function Base.:(==)(c1::Conjunction, c2::Conjunction)
    return c1.conjuncts == c2.conjuncts
end

# function-like evaluation
@inline function (c::Conjunction)(args...)
    evaluate(c, args...)
end

function evaluate(c::Conjunction, args...)
    assert_same_length(c, args...)
    @inbounds for (conjunct, n_args) in c.conjuncts
        if !evaluate(conjunct, args[n_args]...)
            return false
        end
    end
    return true
end


text:
Array{T}(::UndefInitializer, ::Tuple{Infinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{Integer, Infinity}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{Infinity, Integer}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{Integer, Infinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{Infinity, Integer}) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Tuple{Infinity, Infinity}) where {T,N} = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Tuple{Infinity, Infinity}) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Array{T}(::UndefInitializer, ::Infinity, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Infinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Array{T}(::UndefInitializer, ::Integer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Matrix{T}(::UndefInitializer, ::Infinity, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Infinity, ::Integer) where T = throw(ArgumentError("Cannot create infinite Array"))
Matrix{T}(::UndefInitializer, ::Integer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

Vector{T}(::UndefInitializer, ::Tuple{Infinity}) where T = throw(ArgumentError("Cannot create infinite Array"))
Vector{T}(::UndefInitializer, ::Infinity) where T = throw(ArgumentError("Cannot create infinite Array"))

similar(A::AbstractArray, ::Type{T}, axes::Tuple{OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(A::AbstractArray, ::Type{T}, axes::Tuple{OneToInf{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{Infinity}) where T = cache(Zeros{T}(dims))
similar(A::AbstractArray, ::Type{T}, dims::Tuple{Infinity,Infinity}) where T = cache(Zeros{T}(dims))

similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneToInf{Int},OneTo{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, axes::Tuple{OneTo{Int},OneToInf{Int}}) where T = cache(Zeros{T}(axes))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{Infinity}) where T = cache(Zeros{T}(dims))
similar(::Type{<:AbstractArray{T}}, dims::Tuple{Infinity,Infinity}) where T = cache(Zeros{T}(dims))

zeros(::Type{T}, ::Tuple{Infinity}) where T = cache(Zeros{T}(∞))
zeros(::Type{T}, nm::Tuple{Integer, Infinity}) where T = cache(Zeros{T}(nm...))
zeros(::Type{T}, nm::Tuple{Infinity, Integer}) where T = cache(Zeros{T}(nm...))
zeros(::Type{T}, nm::Tuple{Infinity, Infinity}) where T = cache(Zeros{T}(nm...))

fill(x, ::Tuple{Infinity}) = cache(Fill(x,∞))
fill(x, nm::Tuple{Integer, Infinity}) = cache(Fill(x,nm...))
fill(x, nm::Tuple{Infinity, Integer}) = cache(Fill(x,nm...))
fill(x, nm::Tuple{Infinity, Infinity}) = cache(Fill(x,nm...))

# This gets called when infinit number of columns
print_matrix_row(io::IO,
        X::AbstractVecOrMat, A::Vector,
        i::Integer, cols::AbstractVector{<:Infinity}, sep::AbstractString) = nothing

print_matrix_vdots(io::IO, vdots::AbstractString,
        A::Vector, sep::AbstractString, M::Integer, ::NotANumber) = nothing

# Avoid infinite loops on maximum
Base.mapreduce_impl(f, op, A::AbstractArray, ifirst::Integer, ::Infinity) =
    throw(ArgumentError("Cannot call mapreduce on an infinite length $(typeof(A))"))

function show_delim_array(io::IO, itr::AbstractArray, op, delim, cl,
                          delim_one, i1, ::Infinity)
    print(io, op)
    l = 20
    if !show_circular(io, itr)
        recur_io = IOContext(io, :SHOWN_SET => itr)
        if !haskey(io, :compact)
          recur_io = IOContext(recur_io, :compact => true)
        end
        first = true
        i = i1
        if 20 >= i1
          while true
              if !isassigned(itr, i)
                  print(io, undef_ref_str)
              else
                  x = itr[i]
                  show(recur_io, x)
              end
              i += 1
              if i > l
                  print(io, delim)
                  print(io, ' ')
                  print(io, '…')
                  delim_one && first && print(io, delim)
                  break
              end
              first = false
              print(io, delim)
              print(io, ' ')
          end
        end
    end
    print(io, cl)
end

#####
# FillArrays
#####

# Lazy Broadacasting
for typ in (:Ones, :Zeros, :Fill)
    @eval begin
        BroadcastStyle(::Type{$typ{T,N,NTuple{N,<:OneToInf}}}) where {T,N} = LazyArrayStyle{N}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:Any,<:OneToInf}}}) where {T} = LazyArrayStyle{2}()
        BroadcastStyle(::Type{$typ{T,2,<:Tuple{<:OneToInf,<:Any}}}) where {T} = LazyArrayStyle{2}()
    end
end

BroadcastStyle(::Type{<:Diagonal{T,<:AbstractFill{T,1,Tuple{OneToInf{I}}}}}) where {T,I} = LazyArrayStyle{2}()

#####
# Diagonal
#####

BroadcastStyle(::Type{<:Diagonal{<:Any,<:AbstractInfUnitRange}}) = LazyArrayStyle{2}()

#####
# Vcat length
#####

function getindex(f::Vcat{T,1}, k::Infinity) where T
    length(f) == ∞ || throw(BoundsError(f,k))
    ∞
end

_gettail(k, a::Number, b...) = k ≤ 1 ? tuple(a, b...) : _gettail(k - length(a), b...)
_gettail(k, a, b...) = k ≤ length(a) ? tuple(a[k:end], b...) : _gettail(k - length(a), b...)
_vcat(a) = a
_vcat(a, b, c...) = Vcat(a, b, c...)
_unsafe_getindex(::IndexLinear, A::Vcat, r::InfUnitRange) = _vcat(_gettail(first(r), A.args...)...)


text:
module Tst

using Base.Test

import LsoBase

# include ./tst/*
LsoBase.includedir(dirname(@__FILE__)*"/tst")

"""
    tst()

    Run all tests.
"""
function all()
    for exp in names(Tst, true) # iterate over exported members
        if (exp != :Tst && exp != :__META__ && exp != :eval && exp != :all && !startswith(string(exp), "_"))
            eval( :(($exp)()) )  # call exp as function
        end
    end
end

end

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
# Load DataFrame package
require("DataFrames")
using DataFrames

# Load a CSV file into a DataFrame
df = read_table("demo/toy_example.csv")

# Basic indexing
df[1, :]
df["A"]
df[1, "A"]
df[1:2, "A"]
df[1:2, ["A", "B"]]

# Use the with() function to evaluate expressions relative to a DataFrame
with(df, :(A + C))


text:
using AvailablePotentialEnergyFramework, PyCall,PyPlot, Statistics, JLD; tkr = pyimport("matplotlib.ticker")# import FormatStrFormatter

#using , Statistics, DataStructures
#using JLD, SAMtools, ProgressMeter

#const ApeHelperFunctions = SAMtools

data_dir = "/global/cscratch1/sd/aramreye/for_postprocessing/ApeBudgetOutputs/"

# apes_1daysmooth_f5e-4_2km_1000km_homoAll_homoSfc_part_1.jld
function gatherapes(data_dir,exp_name)
            int_ape = Float32[]
            int_ape_ke = Float32[]
            int_ape_dia = Float32[]
            int_ape_rad = Float32[]
            int_ape_ub2 =Float32[]
            int_ape_vb2 = Float32[]
            int_ape_wn2 = Float32[]
            int_ape_rate = Float32[]
            int_ape_Fs = Float32[]
    for iter in 1:24
    file_name = string(exp_name,iter,".jld")
    file_path = joinpath(data_dir,file_name)

        file = jldopen(file_path, "r", mmaparrays=true)
            append!(int_ape,read(file, "int_APE"))
            append!(int_ape_ke,read(file, "int_KE"))
            append!(int_ape_dia,read(file, "int_DIA"))
            append!(int_ape_rad,read(file, "int_RAD"))
            append!(int_ape_ub2,read(file, "int_Ub2"))
            append!(int_ape_vb2,read(file, "int_Vb2"))
            append!(int_ape_wn2,read(file, "int_WN2"))
            append!(int_ape_rate,read(file, "int_APE_rate"))
            append!(int_ape_Fs,read(file, "APE_Fs"))
        close(file)
    end
    residual = int_ape_rate .-  int_ape_rad .-  int_ape_dia .+  int_ape_ub2 .+ int_ape_vb2 .+  int_ape_wn2 .- int_ape_Fs
    return int_ape, int_ape_ke, int_ape_rad, int_ape_dia, int_ape_ub2.+int_ape_vb2, int_ape_wn2, int_ape_Fs, residual, int_ape_rate
end

alldata = [gatherapes(data_dir,"f5e-4_2km_1000km_control")
 gatherapes(data_dir,"f5e-4_2km_1000km_homoRad")
gatherapes(data_dir,"f5e-4_2km_1000km_homoSfc")
gatherapes(data_dir,"f5e-4_2km_1000km_homoRad_homoSfc")];

int_APE_control = gatherapes(data_dir,"f5e-4_2km_1000km_control")
int_APE_homoRad = gatherapes(data_dir,"f5e-4_2km_1000km_homoRad")
int_APE_homoSfc = gatherapes(data_dir,"f5e-4_2km_1000km_homoSfc")
int_APE_homoRad_homoSfc = gatherapes(data_dir,"f5e-4_2km_1000km_homoRad_homoSfc");

size(int_APE_homoRad_homoSfc)

smooth_window = 12
smooth_window_2 = 10
int_APE_control_smooth = filter_array_time.(int_APE_control,smooth_window,1)
int_APE_homoRad_smooth = filter_array_time.(int_APE_homoRad,smooth_window,1)
int_APE_homoSfc_smooth = filter_array_time.(int_APE_homoSfc,smooth_window,1)
int_APE_homoRad_homoSfc_smooth = filter_array_time.(int_APE_homoRad_homoSfc,smooth_window,1);

times = range(0,stop=100,length=1200)

experiment_labels = ["Control" "HomoSfc" "HomoRad" "HomoAll"]
budget_labels = ["APE" "KE" "Rad" "Conv" "Advec" "WN2" "Fs" "Residual" "APE change rate"]
plot_labels = ["a" "b" "c" "d" "e"]

PyPlot.matplotlib.rc("font", family="sans-serif",size=14)
legend_fontsize, plotlabel_fontsize = 12,16
#legend_fontsize, plotlabel_fontsize = 6,8

fig, ax = plt.subplots(5,1,sharex=true,figsize=(18,9))
ax[1].plot(times,1e-5int_APE_control_smooth[1],label = experiment_labels[1])
ax[1].plot(times,1e-5int_APE_homoSfc_smooth[1],label = experiment_labels[2])
ax[1].plot(times,1e-5int_APE_homoRad_smooth[1],label = experiment_labels[3])
ax[1].plot(times,1e-5int_APE_homoRad_homoSfc_smooth[1],label = experiment_labels[4])
ax[1].set_ylabel(string("APE ", L"(J/m^2)",L"x10^5"))
ax[1].yaxis.set_label_coords(-0.05,0.5)
ax[1].yaxis.set_major_formatter(tkr.FormatStrFormatter("%d"))

for variable in 3:7
   ax[2].plot(times,int_APE_control_smooth[variable],label=budget_labels[variable])
end

for variable in 3:7
   ax[3].plot(times,int_APE_homoSfc_smooth[variable])
end
for variable in 3:7
   ax[4].plot(times,int_APE_homoRad_smooth[variable])
end
for variable in 3:7
   ax[5].plot(times,int_APE_homoRad_homoSfc_smooth[variable])
end
for plotid in 1:5
    ax[plotid].annotate(plot_labels[plotid], xy=(0.02, 0.8), xycoords="axes fraction",backgroundcolor="white")
    ax[plotid].grid(b=true,which="both",color="xkcd:gray", linestyle="--",alpha=0.2)
    if plotid > 1
        ax[plotid].text(-0.0, 1.03, experiment_labels[plotid-1], transform=ax[plotid,1].transAxes, size=legend_fontsize)
    end
end

ax[3].set_ylabel(string("APE budget ", L"(W/m^2)"))
ax[3].yaxis.set_label_coords(-0.05,0.25)
ax[5].set_xlabel("Time (days)")

ax[1].legend(loc = (0.1, 1), ncol=4 ,frameon=false,fontsize=legend_fontsize)
ax[2].legend(loc = (0.1, 1), ncol=8 ,frameon=false,fontsize=legend_fontsize)
ax[1].xaxis.set_minor_locator(tkr.MultipleLocator(5))
ax[1].yaxis.set_minor_locator(plt.NullLocator())


text:
using JuMP, Test
using SparseArrays
using LinearAlgebra
using Statistics
using Printf

include("../src/OnePhase.jl")
#Uncomment the below two line of codes when working with HSL is needed
#OnePhase.setUSE_HSL(true)
@show OnePhase.USE_HSL
#OnePhase.loadHSL("../src/linear_system_solvers/")

TEST_MOI = false
#include("CUTEst.jl")
include("problems.jl")
include("kkt_system_solvers.jl")
include("linear_system_solvers.jl")
include("test_moi_nlp_solver.jl")

function unit_tests()
    test_compare_columns()
    test_compute_indicies()
    test_linear_solvers()
    test_kkt_solvers()
end

#function basic_tests(solver)
function basic_tests(options::Dict{String, Any})

    @testset "rosenbrook" begin
       test_rosenbrook1(options)
	   test_rosenbrook2(options)
       test_rosenbrook3(options)
       test_rosenbrook4(options)
    end

    @testset "LP" begin
        test_toy_lp1(options)
        test_toy_lp2(options)
        test_toy_lp3(options)
        test_toy_lp4(options)
        test_toy_lp5(options)
        test_toy_lp6(options)
        test_toy_lp7(options)
    end
    @testset "infeasible" begin
        model = toy_lp_inf1()
	    attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Infeasible

        model = toy_lp_inf2()
        attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Infeasible

        model = circle_nc_inf1()
	    attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Infeasible
    end

    @testset "convex_nlp" begin
        @testset "circle1" begin
            model = circle1()
            attachSolverWithAttributesToJuMPModel(model, options)
            optimize!(model)
            status = MOI.get(model, MOI.TerminationStatus())

            @test status == :Optimal
            check_circle1(model)
        end
        @testset "circle2" begin
            model = circle2()
            attachSolverWithAttributesToJuMPModel(model, options)
            optimize!(model)
            status = MOI.get(model, MOI.TerminationStatus())

            @test status == :Optimal
            check_circle2(model)
        end

        @testset "quad_opt" begin
            model = quad_opt()
	        attachSolverWithAttributesToJuMPModel(model, options)
            optimize!(model)
            status = MOI.get(model, MOI.TerminationStatus())

            if Test.Pass == @test status == :Optimal
                check_quad_opt(model)
            end
        end
    end

    @testset "nonconvex" begin
        model = circle_nc1()
    	attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())

        @test status == :Optimal
        check_circle_nc1(model)

        model = circle_nc2()
	    attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())

        @test status == :Optimal
        check_circle_nc2(model)
    end

    @testset "unbounded_opt_val" begin
        model = lp_unbd()
        attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test :Unbounded == status

        model = circle_nc_unbd()
	    attachSolverWithAttributesToJuMPModel(model, options)
        set_optimizer_attribute(model, "output_level", 0)
        optimize!(model)
	    status = MOI.get(model, MOI.TerminationStatus())
        @test status == :Unbounded

        model = quad_unbd()
	    attachSolverWithAttributesToJuMPModel(model, options)
        optimize!(model)
        status = MOI.get(model, MOI.TerminationStatus())
        @test_broken status == :Unbounded
    end

    @testset "unbounded_feasible_region" begin
        # problems with unbounded feasible region
        test_unbd_feas(options)
    end

    @testset "starting point" begin
        test_starting_point(options,0.5)
        test_starting_point(options,-0.5)
    end

end

function basic_tests()
    #max_it = 100
    max_it = 81
    output_level = 0
    a_norm_penalty = 1e-4
    @testset "basic_tests" begin
        @testset "cholseky linear system solve" begin
	    options = Dict{String, Any}("term!max_it"=>max_it,
	    	"a_norm_penalty"=>a_norm_penalty,
            "output_level"=>output_level,
            "kkt!kkt_solver_type"=>:schur)
	    basic_tests(options)
        end

    if OnePhase.USE_HSL
        @testset "Ma97 linear system solve" begin
	    options = Dict{String, Any}("term!max_it"=>max_it,
	    "a_norm_penalty"=>a_norm_penalty,
            "output_level"=>output_level,
            "kkt!kkt_solver_type"=>:symmetric,
            "kkt!linear_solver_type"=>:HSL)
	    basic_tests(options)
        end

        @testset "Ma97 linear system solve with clever elimination" begin
	    options = Dict{String, Any}("term!max_it"=>max_it,
	    "a_norm_penalty"=>a_norm_penalty,
            "output_level"=>output_level,
            "kkt!kkt_solver_type"=>:clever_symmetric,
            "kkt!linear_solver_type"=>:HSL)
	    basic_tests(options)
        end
    end

        @testset "LDLT julia linear system solve" begin
	    options = Dict{String, Any}("term!max_it"=>max_it,
	    	"a_norm_penalty"=>a_norm_penalty,
            "output_level"=>output_level,
            "kkt!kkt_solver_type"=>:clever_symmetric,
            "kkt!linear_solver_type"=>:julia)
	    basic_tests(options)
        end

    end
end

function moi_nlp_tests()
    @testset "test_moi_nlp_solver" begin
        test_lp1_feasible_MOI()
        test_lp1_optimal_MOI()
        test_nlp1_feasible_MOI()
        test_nlp1_optimal_MOI()
        # test_lp1_feasible_JuMP()
        test_lp1_optimal_JuMP()
        test_nlp1_optimal_JuMP()
    end
end

# lets run the tests!
@time begin
  unit_tests()
  moi_nlp_tests()
  basic_tests()
  #cutest_tests()
  #executeCUTEST_Models()
  #executeCUTEST_Models_benchmark()
end

max_stars_repo_name: ohinder/OnePhase

max_issues_repo_name: ohinder/OnePhase

max_forks_repo_name: ohinder/OnePhase


text:
module Experiment

export trialList2DataFrame
export data_of_run
export mat2jdl

using DataFrames
using MAT
using JLD
using Dates

using ..MyUtils
using ..Neurons

##################################################
#
#
#
##################################################
struct data_of_run
    paradigm::String
    date::Date
    time::DateTime
    sf::Float64
    trial_list::DataFrame
end

##################################################
#
#
#
##################################################
function trialList2DataFrame( data_dict )

    ### orderd according to intended column order in Julia
    ### (f1, f2): f1 -> column index in .mat.TrialList,
    ordered_column_name_pairs = [
        ([ Neurons.recording() for ii in 1:size(data_dict["TrialList"],1) ], "Neurons"),
        ([ [ data_dict["EyeX"][.!isnan.(data_dict["EyeX"][:,trl]),trl], data_dict["EyeY"][.!isnan.(data_dict["EyeY"][:,trl]),trl] ] for trl in 1:size(data_dict["EyeY"],2) ], "EyeXY"),
        (fill([[]],size(data_dict["TrialList"],1)), "Saccades"),
        (fill([[]],size(data_dict["TrialList"],1)), "Saccades_clean"),
        (fill([[]],size(data_dict["TrialList"],1)), "Fixations"),
        (fill([],size(data_dict["TrialList"],1)), "LookedAtTarget"),
        (data_dict["TrialList"][:,9], "TrialStartTime"),
        (data_dict["TrialList"][:,10], "TrialEndTime"),
        (fill(false, size(data_dict["TrialList"],1)), "Valid"),
        (data_dict["TrialList"][:,24].==1, "Hit"),
        (fill(false, size(data_dict["TrialList"],1)), "FA"),
        (fill(false, size(data_dict["TrialList"],1)), "Invalid"),
        (Array{Float64,1}(undef,size(data_dict["TrialList"],1)), "RTfirstFixation"),
        (Array{Float64,1}(undef,size(data_dict["TrialList"],1)), "RTcorrectFixation"),
        (data_dict["TrialList"][:,23].==1, "WasMissed"),                    ### !!! if this row changes, changes the part in mat2jdl() where missed trials are filtered out !!!
        (data_dict["TrialList"][:,22].==1, "WasTimeout"),
        (data_dict["TrialList"][:,16].+1, "ValidTargetIndex"),
        (data_dict["TrialList"][:,12].+1, "CueTargetIndex"),
        (data_dict["TrialList"][:,17].==1, "WasCueValid"),
        # (data_dict["TrialList"][:,18].==1, "WasFixationFailureFixPt"),    ### redundant with WasMissed
        # (data_dict["TrialList"][:,19].==1, "WasFixationFailureTarget"),   ### redundant with WasTimeout
        (data_dict["TrialList"][:,14], "LuminanceChange"),
        (data_dict["TrialList"][:,11], "LuminanceChangeTime"),
        (data_dict["TrialList"][:,21], "PauseDuration"),
        (data_dict["TrialList"][:,13], "ITI"),
        (data_dict["TrialList"][:,15], "NumberOfNonFixFailTrials"),
        (data_dict["TrialList"][:,20].==1, "WasLEDChangeActive")]

    ### -- Now in Behavior.jl --- ###
    # ### extract eye trace to store events in TrialList
    # xy = [[data_dict["EyeX"][:,trl], data_dict["EyeY"][:,trl]] for trl in 1:size(data_dict["EyeY"], 2)];
    # [deleteat!(xy[trl][dim], findall(isnan.(xy[trl][dim]))) for dim in [1,2], trl in 1:size(data_dict["EyeY"], 2)] # remove NaNs
    # sf = data_dict["samplingRate"]

    ### generate dataframe from ordered_column_name_pairs
    trial_list = DataFrame([el[1] for el in ordered_column_name_pairs], [Symbol(el[2]) for el in ordered_column_name_pairs])

    return trial_list
end

##################################################
#
# reads .mat file and generates raw_data stuct
#
##################################################
function mat2jdl(ifi, ofi)
    for (idx, mf_) in enumerate(ifi)

        ### read mat files of runs specified in paths_to_mat_files and store the in data_dicts array
        data_dict = matread(mf_)
        # @show data_dict
        #data_dict_ = deepcopy(data_dict)

        ### remove missed trials in data_dict
        missed_trials_idx = findall(data_dict["TrialList"][:, 23].==1)
        ks = ["TargetY", "VoltX", "TargetX", "EyeT", "TargetTime", "EyeX", "EyeY", "VoltY"]
        for k in ks
            data_dict[k] = data_dict[k][:, setdiff(1:end, missed_trials_idx)]
        end
        data_dict["TrialList"] = data_dict["TrialList"][setdiff(1:end, missed_trials_idx), :]

        ### remove trials with all NaN
        xy = [[data_dict["EyeX"][:,trl], data_dict["EyeY"][:,trl]] for trl in 1:size(data_dict["EyeY"], 2)];
        only_nans = [trl for trl in 1:size(xy, 1) if all(isnan.(xy[trl][1]))] # find all NaN trials
        for k in ks
            data_dict[k] = data_dict[k][:, setdiff(1:end, only_nans)]
        end
        data_dict["TrialList"] = data_dict["TrialList"][setdiff(1:end, only_nans), :];

        print(".mat file loaded: ", idx,"/",size(ifi)[1])

        ### convert trialList from mat file to julia table
        trial_list = trialList2DataFrame(data_dict)
        print(" | trial_list extracted")

        ## get indices of valid and invalid trials
        # valid_trials_idx = findall( trial_list[!, :Valid] )

        raw_data = data_of_run(
            "...",
            Date(mf_[end-26:end-17],"y-m-d"),
            DateTime(mf_[end-15:end-8],"H-M-S"),
            data_dict["samplingRate"],
            trial_list
        )

        ### store processed data to file
        save(joinpath(ofi, mf_[end-41:end-4]*".jld"), "raw_data", raw_data)
        println(" | .jdl saved")

        #return eye_dat
    end
end

end

max_stars_repo_name: maalaria/Newbe

max_issues_repo_name: maalaria/Newbe

max_forks_repo_name: maalaria/Newbe


text:
function OneWayAnova()
    return OneWayAnova(())
end

function anova_f_value(obj::OneWayAnova, arg0::Collection)
    return jcall(obj, "anovaFValue", jdouble, (Collection,), arg0)
end

function anova_p_value(obj::OneWayAnova, arg0::Collection)
    return jcall(obj, "anovaPValue", jdouble, (Collection,), arg0)
end

function anova_p_value(obj::OneWayAnova, arg0::Collection, arg1::jboolean)
    return jcall(obj, "anovaPValue", jdouble, (Collection, jboolean), arg0, arg1)
end

function anova_test(obj::OneWayAnova, arg0::Collection, arg1::jdouble)
    return jcall(obj, "anovaTest", jboolean, (Collection, jdouble), arg0, arg1)
end


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

include("rs-utils.jl")

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

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

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

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

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

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

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

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

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

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

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

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

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

        # retrieve rgb data
        rgb_frame_data = Ptr{UInt8}(rs2_get_frame_data(frame, err))
        check_error(err)
        rgb_image_raw = unsafe_wrap(Array, rgb_frame_data, 3*WIDTH*HEIGHT)

        frame_number = rs2_get_frame_number(frame, err)
        check_error(err)

        frame_timestamp = rs2_get_frame_timestamp(frame, err)
        check_error(err)

        frame_timestamp_domain = rs2_get_frame_timestamp_domain(frame, err)
        check_error(err)
        frame_timestamp_domain_str = unsafe_string(rs2_timestamp_domain_to_string(frame_timestamp_domain))
        frame_metadata_time_of_arrival = rs2_get_frame_metadata(frame, RS2_FRAME_METADATA_TIME_OF_ARRIVAL, err)
        check_error(err)

        println("RGB frame arrived.")
        println("First 10 bytes: $(rgb_image_raw[1:10])")
        println("Frame No.$frame_number")
        println("Timestamp: $frame_timestamp")
        println("Timestamp domain: $frame_timestamp_domain_str")
        println("Time of arrival: $frame_metadata_time_of_arrival")

        rs2_release_frame(frame)
    end
    rs2_release_frame(frames)
end

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

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


text:
# Create an OOI data request URL from input parameters.
#
# @return A URL string that can be input into ooi_submit_request.
#
# @param site An eight character OOI site.
# @param node A five character OOI node.
# @param instrument A twelve character OOI instrument.
#   The third character must be a dash (-)
# @param method Method for data delivery.
# @param stream  Stream for data.
# @param start_date UTC in format YYYY-mm-dd.
# @param start_time UTC in format HH:MM:SS.
# @param stop_date UTC in format YYYY-mm-dd.
# @param stop_time UTC in format HH:MM:SS.

function ooi_create_url(site,node,instrument,method,stream;
                        start_date ="2010-01-01",start_time = "00:00:00",
                        stop_date= "2040-12-31",stop_time="23:59:59")

    url = string("https://", user,":",token,
                 "@ooinet.oceanobservatories.org/api/m2m/12576/sensor/inv/",
                 site,"/",node,"/",instrument,"/",method,"/",stream,
                 "?beginDT=",start_date,"T",start_time,
                 ".000Z&endDT=",stop_date,"T",stop_time,".999Z")

    return url
end  #End of function.

max_stars_repo_name: lstoltz/ooi-data-explorations

max_issues_repo_name: lstoltz/ooi-data-explorations

max_forks_repo_name: lstoltz/ooi-data-explorations


text:
using PyLCM
@pyimport microstrain
@pyimport bot_core

lc = LCM()

function handle_msg(channel, msg_data)
    global lc
    @time msg = microstrain.ins_t[:decode](msg_data)
    @show msg[:utime]

    # msg = microstrain.ins_t()
    nmsg = deepcopy(msg)
    nmsg[:utime] = 12345
    publish(lc, "NEWMSG", nmsg)
end

subscribe(lc, "MICROSTRAIN_IMU", handle_msg)

while true
    handle(lc)
end

function handle_cam(channel, msg_data)
  global lc
  msg = bot_core.image_t[:decode](msg_data)
  @show msg[:width], msg[:height], msg[:pixelformat], msg[:size], typeof(msg[:data]), msg[:data][1]
  ndata = Vector{Int}(msg[:size])
  d = msg[:data]
  for i in 1:5#msg[:size]
    @show convert(Vector{UInt8},d[i]) #ndata[i] = Int(msg[:data][i])
  end
  @show size(ndata)
  publish(lc, "NEWIMAGE", msg)
  nothing
end

subscribe(lc, "CAMLCM_IMAGE", handle_cam)

# for ev in reader
#   ev[:channel] == "MICROSTRAIN_INS" ? println(ev) : nothing
# end

# using PyCall
# pylcm = pyimport("lcm")
# pylcm[:EventLog]("/home/dehann/data/inertialtest/df-backtozero/6point-01.lcm","r")
# for ev in reader
#          ev[:channel] == "MICROSTRAIN_INS" ? println(microstrain.ins_t[:decode](convert(Vector{UInt8}, ev[:data]))) : nothing
#        end

pkgs = [
  "IJulia";
  "Color";
  "Graphs";
  "GraphViz";
  "Gadfly";
  "DataFrames";
  "Distributions";
  "NLsolve";
  "PyCall";
  "PyPlot";
  "HDF5";
  "JLD";
  "JuMP";
  "NLopt";
  "RobotOS"
]

for p in pkgs; Pkg.add(p); end


text:
function almostequal(A,B)

if norm(A-B) < 1e-10
    return true
else
    return false
end

end


text:
# ---
# title: Custom Title
# description: >
#   this is a single line
#   destiption that spans
#   over multiple lines
# ---

# This is the content


text:
import Gmsh: gmsh

function generate_mesh_inclusion(; Lx=1.0, Ly=2.0 ,radius=Lx/40 ,lowres=0.02, highres=0.001, file)
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.add("model")

    xmin, xmax = 0.0, Lx
    ymin, ymax = 0.0, Ly

    p1 = gmsh.model.geo.addPoint(xmin, ymin, 0, lowres)
    p2 = gmsh.model.geo.addPoint(xmax, ymin,  0, lowres)
    p3 = gmsh.model.geo.addPoint(xmax, ymax, 0, lowres)
    p4 = gmsh.model.geo.addPoint(xmin, ymax, 0, lowres)
    l1 = gmsh.model.geo.addLine(p1, p2)
    l2 = gmsh.model.geo.addLine(p2, p3)
    l3 = gmsh.model.geo.addLine(p3, p4)
    l4 = gmsh.model.geo.addLine(p4, p1)
    cl1 = gmsh.model.geo.addCurveLoop([l1, l2, l3, l4])
    # Add two points that will be used as fixed in ux :
    # ptop = gmsh.model.geo.addPoint(0.5, 2, 0, lowres)
    # pbot = gmsh.model.geo.addPoint(0.5, 0, 0, lowres)
    # embed these points in the surface to force the mesh to define nodes at their location.
    # pbot = gmsh.model.geo.addPoint(xmax/2, ymin, 0, lowres)
    # ptop = gmsh.model.geo.addPoint(xmax/2, ymax, 0, lowres)

    #pmid = gmsh.model.geo.addPoint(xmax/2, ymax/2, 0, highres)

    #circ = gmsh.model.occ.addCircle(xmax/2, ymax/2, 0.0, xmax/20)
    pmid = gmsh.model.geo.addPoint(xmax/2, ymax/2, 0, 5*highres)
    pc1 = gmsh.model.geo.addPoint(xmax/2, ymax/2 - radius, 0, 2highres)
    pc2 = gmsh.model.geo.addPoint(xmax/2 + radius, ymax/2, 0, 2highres)
    pc3 = gmsh.model.geo.addPoint(xmax/2, ymax/2 + radius, 0, 2highres)
    pc4 = gmsh.model.geo.addPoint(xmax/2 - radius, ymax/2, 0, 2highres)
    arc1 = gmsh.model.geo.addCircleArc(pc1,pmid,pc2)
    arc2 = gmsh.model.geo.addCircleArc(pc2,pmid,pc3)
    arc3 = gmsh.model.geo.addCircleArc(pc3,pmid,pc4)
    arc4 = gmsh.model.geo.addCircleArc(pc4,pmid,pc1)
    cl2 = gmsh.model.geo.addCurveLoop([arc1,arc2,arc3,arc4])
    s1 = gmsh.model.geo.addPlaneSurface([cl1,cl2])
    s2 = gmsh.model.geo.addPlaneSurface([cl2])
    #gmsh.model.geo.synchronize()
    # gmsh.model.mesh.embed(0, [pbot], 1, l1)
    # gmsh.model.mesh.embed(0, [ptop], 1, l3)
    #gmsh.model.mesh.embed(2, [s2], 2, s1)

    #clamped = gmsh.model.addPhysicalGroup(0, [ptop, pbot])
    bottom = gmsh.model.addPhysicalGroup(1, [l1])
    right = gmsh.model.addPhysicalGroup(1, [l2])
    top = gmsh.model.addPhysicalGroup(1, [l3])
    left = gmsh.model.addPhysicalGroup(1, [l4])
    circle = gmsh.model.addPhysicalGroup(1, [cl2])
    domain1 = gmsh.model.addPhysicalGroup(2, [s1])
    domain2 = gmsh.model.addPhysicalGroup(2, [s2])

    #gmsh.model.setPhysicalName(0, clamped, "clamped")
    gmsh.model.setPhysicalName(1, bottom, "bottom")
    gmsh.model.setPhysicalName(1, right, "right")
    gmsh.model.setPhysicalName(1, top, "top")
    gmsh.model.setPhysicalName(1, left, "left")
    gmsh.model.setPhysicalName(1, circle, "circle")
    gmsh.model.setPhysicalName(2, domain1, "domain1")
    gmsh.model.setPhysicalName(2, domain2, "domain2")
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(2)
    gmsh.write(file)
    gmsh.finalize()

    return nothing
end


text:
struct PSResult{T}
    value::T
end

struct PSTextMagicAPIUserObject

    # data -
    userName::String
    apiKey::String
    label::String

    function PSTextMagicAPIUserObject(userName::String, apiKey::String, label::String)
        this = new(userName,apiKey,label)
    end
end


text:
using Random
using TestAbstractTypes

# Set reproducible random number seed
Random.seed!(0)
testAbstractChar(Char, generator(Char))

# Set reproducible random number seed
Random.seed!(0)
testAbstractString(String, generator(String))

# Set reproducible random number seed
Random.seed!(0)
testAbstractString(SubString, generator(SubString))


text:
function http_get_call_with_url(url::String)::String

    # should we check if this string is formatted as a URL?
    if (occursin("http://",url) == false)
        throw(ArgumentError("url $(url) is not properly formatted"))
    end

    # ok, so we are going to make a HTTP GET call with the URL that was passed in -
    response = HTTP.request("GET",url)

    # ok, so let's check if we are getting a 200 back -
    if (response.status == 200)
        return String(response.body)
    else
        # create an error, and throw it back to the caller -
        throw(ErrorException("http status flag $(response.status) was returned from url $(url)"))
    end
end


text:
export restrictCellCenteredVariables,restrictNodalVariables2,getFWInterp;

# Bilinear "Full Weighting" prolongation
function getFWInterp(n_nodes::Array{Int64,1},geometric::Bool=false)
# n here is the number of NODES!!!
(P1,nc1) = get1DFWInterp(n_nodes[1],geometric);
(P2,nc2) = get1DFWInterp(n_nodes[2],geometric);
if length(n_nodes)==3
	(P3,nc3) = get1DFWInterp(n_nodes[3],geometric);
end
if length(n_nodes)==2
	P = kron(P2,P1);
	nc = [nc1,nc2];
else
	P = kron(P3,kron(P2,P1));
	nc = [nc1,nc2,nc3];
end
return P,nc
end

function get1DFWInterp(n_nodes::Int64,geometric)
# n here is the number of NODES!!!
oddDim = mod(n_nodes,2);
if n_nodes > 2
	halfVec = 0.5*ones(n_nodes-1);
	P = spdiagm(-1=>halfVec,0=>ones(n_nodes),1=>halfVec); #used to be P = spdiagm((halfVec,ones(n_nodes),halfVec),[-1,0,1],n_nodes,n_nodes);
    if oddDim == 1
        P = P[:,1:2:end];
    else
		if geometric
			P = sparse(1.0I,n_nodes,n_nodes);
			println("Warning: getFWInterp(): in geometric mode we stop coarsening because num cells does not divide by two");
		else
			P = P[:,[1:2:end;end]];
			P[end-1:end,end-1:end] = speye(2);
#         	P = P[:,1:2:end];
#         	P[end,end-1:end] = [-0.5,1.5];
		end
    end
else
    P = sparse(1.0I,n_nodes,n_nodes);
end
nc = size(P,2);
return P,nc
end

########################### GEOMETRIC MULTIGRID STUFF ###############################

function restrictCellCenteredVariables(rho::Array,n::Array{Int64})
R,nc = getRestrictionCellCentered(n);
rho_c = R*rho[:];
return rho_c;
## TODO: make this more efficient...
end

export restrictNodalVariables
function restrictNodalVariables(rho::Array,n_nodes::Array{Int64})
P,nc = getFWInterp(n_nodes,true);
rho_c = zeros(eltype(rho),size(P,2));
rho_c[:] = (0.5^length(n_nodes))*(P'*rho[:]);
return rho_c;
## TODO: make this more efficient...
end

function restrictNodalVariables2(rho::Array,n_nodes::Array{Int64})
P,nc = getFWInterp(n_nodes,true);
R1,nc1 = get1DNodeFullWeightRestriction(n_nodes[1]-1);
R2,nc2 = get1DNodeFullWeightRestriction(n_nodes[2]-1);
R = kron(R2,R1);
# rho_c = zeros(eltype(rho),size(P,2));
# println(size(R))
# println(size(rho[:]))

rho_c = R*rho[:];
return rho_c;
## TODO: make this more efficient...
end