text:
# Julia wrapper for header: /Users/dls/jl/openssl-1.0.1e/include/openssl/des.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

@c Ptr{Uint8} DES_options () libcrypto
@c None DES_ecb3_encrypt (Ptr{const_DES_cblock}, Ptr{DES_cblock}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Int32) libcrypto
@c Uint32 DES_cbc_cksum (Ptr{Uint8}, Ptr{DES_cblock}, Int64, Ptr{DES_key_schedule}, Ptr{const_DES_cblock}) libcrypto
@c None DES_cbc_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Int32) libcrypto
@c None DES_ncbc_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Int32) libcrypto
@c None DES_xcbc_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Ptr{const_DES_cblock}, Ptr{const_DES_cblock}, Int32) libcrypto
@c None DES_cfb_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int32, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Int32) libcrypto
@c None DES_ecb_encrypt (Ptr{Uint8}, Ptr{Uint8}, Ptr{DES_key_schedule}, Int32) libcrypto
@c None DES_encrypt1 (Ptr{Uint32}, Ptr{DES_key_schedule}, Int32) libcrypto
@c None DES_encrypt2 (Ptr{Uint32}, Ptr{DES_key_schedule}, Int32) libcrypto
@c None DES_encrypt3 (Ptr{Uint32}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}) libcrypto
@c None DES_decrypt3 (Ptr{Uint32}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}) libcrypto
@c None DES_ede3_cbc_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Int32) libcrypto
@c None DES_ede3_cbcm_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Ptr{DES_cblock}, Int32) libcrypto
@c None DES_ede3_cfb64_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Ptr{Int32}, Int32) libcrypto
@c None DES_ede3_cfb_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int32, Int64, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Int32) libcrypto
@c None DES_ede3_ofb64_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Ptr{Int32}) libcrypto
@c Int32 DES_enc_read (Int32, Ptr{None}, Int32, Ptr{DES_key_schedule}, Ptr{DES_cblock}) libcrypto
@c Int32 DES_enc_write (Int32, Ptr{None}, Int32, Ptr{DES_key_schedule}, Ptr{DES_cblock}) libcrypto
@c Ptr{Uint8} DES_fcrypt (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c Ptr{Uint8} DES_crypt (Ptr{Uint8}, Ptr{Uint8}) libcrypto
@c None DES_ofb_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int32, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}) libcrypto
@c None DES_pcbc_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Int32) libcrypto
@c Uint32 DES_quad_cksum (Ptr{Uint8}, Ptr{DES_cblock}, Int64, Int32, Ptr{DES_cblock}) libcrypto
@c Int32 DES_random_key (Ptr{DES_cblock},) libcrypto
@c None DES_set_odd_parity (Ptr{DES_cblock},) libcrypto
@c Int32 DES_check_key_parity (Ptr{const_DES_cblock},) libcrypto
@c Int32 DES_is_weak_key (Ptr{const_DES_cblock},) libcrypto
@c Int32 DES_set_key (Ptr{const_DES_cblock}, Ptr{DES_key_schedule}) libcrypto
@c Int32 DES_key_sched (Ptr{const_DES_cblock}, Ptr{DES_key_schedule}) libcrypto
@c Int32 DES_set_key_checked (Ptr{const_DES_cblock}, Ptr{DES_key_schedule}) libcrypto
@c None DES_set_key_unchecked (Ptr{const_DES_cblock}, Ptr{DES_key_schedule}) libcrypto
@c None DES_string_to_key (Ptr{Uint8}, Ptr{DES_cblock}) libcrypto
@c None DES_string_to_2keys (Ptr{Uint8}, Ptr{DES_cblock}, Ptr{DES_cblock}) libcrypto
@c None DES_cfb64_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Ptr{Int32}, Int32) libcrypto
@c None DES_ofb64_encrypt (Ptr{Uint8}, Ptr{Uint8}, Int64, Ptr{DES_key_schedule}, Ptr{DES_cblock}, Ptr{Int32}) libcrypto
@c Int32 DES_read_password (Ptr{DES_cblock}, Ptr{Uint8}, Int32) libcrypto
@c Int32 DES_read_2passwords (Ptr{DES_cblock}, Ptr{DES_cblock}, Ptr{Uint8}, Int32) libcrypto


text:
# Runs through all ITMIX glaciers, updates cached files and plots all.
using BITEModel
const BM=BITEModel
plotyes = false
global glacier, gl,gb,pp,pm,pn,pl,hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio
gl_with_errors = []
gl_without_errors = []
glns = 1:1606
t_tot = @elapsed for id in glns
    update_cache = false
    pl_kws = Dict{Symbol,Any}()
    verbose = false
    println(" ")
    println("Penisula glacier $id")
    gid = BM.PeninsulaGlacier(id)
    print("Loading ")
    try
        pl = BM.LoadPara(gid, update_cache=update_cache)
        pl.dataset_opts[BM.ParaData][:pn][:calc_boxcar_M] = false
        print("Loading ... ")
        @time gl,pp,pm,pn,pl = load_glacier(gid,pl)
        print("Elevation bands ... ")
        @time gb,pm = BM.make_bands(gl, pp, pm, pn)

        print("Running fwd-model ")
        @time hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, vol_ratio, gb =
            BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)
        mh = round(Int,sum(hs1d[1:end-1].*gb.ws.*gb.ls)/BM.area(gl))
        mhva = round(Int, BM.volume_area_mean_h(BM.area(gl), isicecap=true))
        println("Mean h: $mh ; Volume-area: $mhva; Glacier area: $(BM.area(gl)/1e6)km^2")
        push!(gl_without_errors, (id, mh, mhva, BM.area(gl)/1e6))
        # println("Plotting")
        # display(BM.plot_run(gb, hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l, pp, pm,
        #                     reuse=false, plot2d=[:h,:iv,:ivlog][1]))
    catch e
        e==InterruptException() && break
        push!(gl_with_errors, (id, e))
        warn("ERROR $e")
        println("")
    end
end

println("\n\nTotal time $(round(t_tot/60,2)) minutes; average per glacier $(round(t_tot/length(glns),3))s.")
println("$(length(gl_with_errors)) out of $(length(glns)) glaciers had errors")

ext: jl


text:
#!/usr/bin/env julia
isdir(Pkg.dir("YAML")) || Pkg.add("YAML")

module ClassifyPackages
using YAML

cd(dirname(@__FILE__))
allnames = filter!(p -> isdir(p) && p != ".git" &&
    p != "METADATA", readdir())

type PackageInfo
    name::String
    hastravis::Bool
    travisyml::Dict
    language::String
    usesppa::Bool
    jlminver::VersionNumber
    jlmaxver::VersionNumber # looking just at package require, not
# metadata require, so not a complete indicator of package deprecation
    lastcommit::DateTime
end

pkgs = Array(PackageInfo, length(allnames))

for (i, p) in enumerate(allnames)
    pkgs[i] = PackageInfo(p, false, Dict(), "", false,
        v"0.0", typemax(VersionNumber), now())
    ymlpath = joinpath(p, ".travis.yml")
    if isfile(ymlpath)
        pkgs[i].hastravis = true
        ymlstring = readstring(ymlpath)
        yml = YAML.load(ymlstring)
        pkgs[i].travisyml = yml
        if haskey(yml, "language")
            pkgs[i].language = yml["language"]
        end
        pkgs[i].usesppa = contains(ymlstring, "ppa:staticfloat")
    end
    requirepath = joinpath(p, "REQUIRE")
    if isfile(requirepath)
        reqs = Pkg.Reqs.parse(requirepath)
        if haskey(reqs, "julia")
            juliareq = reqs["julia"]
            if length(juliareq.intervals) == 1
                interval = juliareq.intervals[1]
                pkgs[i].jlminver = interval.lower
                pkgs[i].jlmaxver = interval.upper
            else
                error("complicated julia requirement for $p: $juliareq")
            end
        end
    end
    cd(p) do
        pkgs[i].lastcommit = Dates.unix2datetime(parse(Int,
            readchomp(`git log -1 --format=%ct`)))
    end
end

#=
# interesting notes:
for f in fieldnames(PackageInfo)
    @eval $f = map(p -> p.$f, pkgs)
end
# or language = getfield.(pkgs, [:language])
unique(language) # julia, cpp, ruby, none
name[language .== "ruby"] # only MsgPackRpcServer
name[(language .== "cpp") & ~usesppa] # only SystemImageBuilder
name[(language .== "") & hastravis] # only Cxx
name[(language .== "julia") & usesppa] # BuildExecutable and Playground for julia-deps patchelf

name[~haskey.(travisyml, ["julia"]) & (language .== "julia")]
# 7 packages: 5 with matrices, CauseMap and LLLplus missing julia entries
=#

end # module

max_stars_repo_name: JuliaTagBot/cloneallpackages

max_issues_repo_name: JuliaTagBot/cloneallpackages

max_forks_repo_name: JuliaTagBot/cloneallpackages


text:
# Romeo LTI utils
#
# Various utility functions.

"""
    kernel_to_signel(k, Δt)

Given a kernel vector `k` compute the corresponding response signal.

This function is an unfortunate necesity of the discretization process, and the consequence
of the way the continuous-time convolutions are implemented using sampled signals. Input
particular, if the response of a system to a given excitation `u` is computed as `g ⋆ u`,
with `⋆` denoting discrete convolution, then the continuous-time impulse response of the
system is not `g`, but `g/Δt`.
"""
kernel_to_signal(k::Array{<:Number}, Δt::Real) = k / Δt

"""
    impulse_resp([T=typeof(1.0)], sys, N, Δt)

Given a system `sys` compute its impulse response vector of length `N`, uniformly
discretized with sample time `Δt`.
"""
impulse_resp(T::Type{<:Number}, sys::SisoLtiSystem, N::Integer, Δt::Real) = kernel_to_signal(kernel(T, sys, N, Δt), Δt)
impulse_resp(sys::SisoLtiSystem, N::Integer, Δt::Real) = impulse_resp(typeof(1.0), sys::SisoLtiSystem, N::Integer, Δt::Real)

"""
    invlaplace([T=typeof(1.0)], sys, N, Δt)

Given a Laplace domain operator (a system) `sys` compute its inverse Laplace transform
as a vector of length `N`, uniformly discretized with sample time `Δt`.
"""
invlaplace(T::Type{<:Number}, sys::SisoLtiSystem, N::Integer, Δt::Real) = impulse_resp(T, sys, N, Δt)
invlaplace(sys::SisoLtiSystem, N::Integer, Δt::Real) = impulse_resp(sys, N, Δt)

"""
    step_resp([T=typeof(1.0)], sys, N, Δt)

Given a system `sys` compute its step response vector of length `N`, uniformly
discretized with sample time `Δt`.
"""
step_resp(T::Type{<:Number}, sys::SisoLtiSystem, N::Integer, Δt::Real) = impulse_resp(T, sys / Diff(1), N, Δt)
step_resp(sys::SisoLtiSystem, N::Integer, Δt::Real) = step_resp(Float64, sys, N, Δt)

"""
    simulate(sys, input, Δt)

Simulate response of `sys` to input signal `input`, uniformly sampled with `Δt`.
"""
function simulate(sys::SisoLtiSystem, input::Array{T}, Δt::Real) where T <: Number
    convolve(kernel(T, sys, length(input), Δt), input)
end

export kernel_to_signal, impulse_resp, step_resp, simulate, invlaplace


text:
module Flux

using Base: tail
using LinearAlgebra, Statistics, Random  # standard lib
using MacroTools, Reexport, ProgressLogging, SpecialFunctions
using MacroTools: @forward

@reexport using NNlib
using MLUtils
import Optimisers: Optimisers, trainable, destructure  # before v0.13, Flux owned these functions

using Zygote, ChainRulesCore
using Zygote: Params, @adjoint, gradient, pullback, @nograd
export gradient

# Pirate error to catch a common mistake. (Internal function `base` because overloading `update!` is more likely to give ambiguities.)
Optimisers.base(dx::Zygote.Grads) = error("Optimisers.jl cannot be used with Zygote.jl's implicit gradients, `Params` & `Grads`")

export Chain, Dense, Maxout, SkipConnection, Parallel,
       RNN, LSTM, GRU, GRUv3,
       SamePad, Conv, CrossCor, ConvTranspose, DepthwiseConv,
       AdaptiveMaxPool, AdaptiveMeanPool, GlobalMaxPool, GlobalMeanPool, MaxPool, MeanPool,
       Dropout, AlphaDropout, LayerNorm, BatchNorm, InstanceNorm, GroupNorm,
       Upsample, PixelShuffle,
       fmap, cpu, gpu, f32, f64,
       testmode!, trainmode!

include("optimise/Optimise.jl")
using .Optimise
using .Optimise: @epochs
using .Optimise: skip
export Descent, ADAM, Momentum, Nesterov, RMSProp,
  ADAGrad, AdaMax, ADADelta, AMSGrad, NADAM, OADAM,
  ADAMW, RADAM, AdaBelief, InvDecay, ExpDecay,
  WeightDecay, ClipValue, ClipNorm

using CUDA
const use_cuda = Ref{Union{Nothing,Bool}}(nothing)

include("utils.jl")
include("onehot.jl")
include("functor.jl")

# Pirate error to catch a common mistake.
Functors.functor(::Type{<:MLUtils.DataLoader}, x) = error("`DataLoader` does not support Functors.jl, thus functions like `Flux.gpu` will not act on its contents.")

include("layers/stateless.jl")
include("layers/basic.jl")
include("layers/conv.jl")
include("layers/recurrent.jl")
include("layers/normalise.jl")
include("layers/upsample.jl")
include("layers/show.jl")

include("loading.jl")

include("outputsize.jl")

include("data/Data.jl")
using .Data

include("losses/Losses.jl")
using .Losses # TODO: stop importing Losses in Flux's namespace in v0.12

include("deprecations.jl")

include("cuda/cuda.jl")

end # module


text:
module GnuplotRecipes

using Gnuplot, Tables, Measurements, InvertedIndices

export bars

"Helper function returning the number of columns in the table."
ncol(table) = length(Tables.columnnames(table))

"Helper function returning the number of rows in the table."
nrow(table) = Tables.rowcount(table)

"""
    bars(table; kwargs...)::Vector{Gnuplot.PlotElement}

Explicit recipe to plot bar graphs from tabular data.
`Measurement`-type data are plot with errorbars style.

# Arguments

- `table`: data to plot. The first column should contains text labels,
  the other columns the plotted values. If the values are of
  `Measurement` type, they will be plotted with errorbars style,
  unless overridden with `hist_style`.
- `box_width=0.8`: the width of the boxes. One means that the boxes
  touch each other.
- `gap::Union{Int64,Nothing}=nothing`: The gap between bar clusters.
  If `nothing`, it is set automatically depending on the number of
  bars in the cluster; to zero for one bar in the cluster, to 1 for
  multiple bars.
- `hist_style=nothing`: histogram style — see Gnuplot documentation.
- `fill_style="solid 1 border -1"`: fill style — see Gnuplot documentation.
- `errorbars=""`: errorbars style — see Gnuplot documentation.
- `label_rot=-30`: label rotation angle; if > 0, align label to right.
- `label_enhanced=false`: whether to apply Gnuplot enhanced formatting to labels.
- `key_enhanced=false`: whether to apply Gnuplot enhanced formatting to data keys.
- `y2cols=[]`: Columns (specified as symbols) which should be plot against *y2* axis.
- `linetypes=1:ncol(table)-1`: Line types (colors) used for different bars
- `xlabelcol=1`: Index of columns containing labels.

# Example

```jldoctest
julia> using Measurements, Gnuplot, DataFrames

julia> table = DataFrame(names=["very long label", "b", "c"],
                         temp=10:12,
                         speed=collect(4:-1:2) .± 1)
3×3 DataFrame
 Row │ names            temp   speed
     │ String           Int64  Measurem…
─────┼───────────────────────────────────
   1 │ very long label     10    4.0±1.0
   2 │ b                   11    3.0±1.0
   3 │ c                   12    2.0±1.0

julia> @gp bars(table)

```
"""
function bars(table;
              box_width = 0.8,
              gap::Union{Int64,Nothing} = nothing,
              hist_style = nothing,
              fill_style = "solid 1 border -1",
              errorbars = "",
              label_rot = -30,
              label_enhanced = false,
              key_enhanced = false,
              y2cols = [],
              linetypes = 1:ncol(table)-1,
              xlabelcol = 1,
              )::Gnuplot.PlotElement
    nr = nrow(table)

    foreach(y2cols) do col
        @assert col in propertynames(table)
    end

    axes(group) = (group in y2cols) ? "axes x1y2 " : ""

    if hist_style === nothing
        hist_style = (any(map(eltype.(Tables.columns(table))) do type type <: Measurement end)
                      ? "errorbars" : "clustered")
    end
    eb = hist_style == "errorbars"

    if gap === nothing
        gap = ncol(table) == 2 ? 0 : 1
    end
    gap_str = (hist_style ∈ ["clustered", "errorbars"]) ? "gap $gap" : ""
    hist_style == "columnstacked" && @warn "columnstacked is not well supported"

    function gpusing(i)
        cols = [i+1, "xticlabels(1)"]
        eb && insert!(cols, 2, i+ncol(table))
        join(cols, ":")
    end

    label_column = Tables.columns(table)[xlabelcol]
    value_columns = Tables.columns(table)[Not(xlabelcol)]
    group_names = Tables.columnnames(value_columns)

    data=Gnuplot.DatasetText(
        String.(label_column), # labels
        [Measurements.value.(column) for column in value_columns]...,
        [eltype(column) <: Measurement ?
            Measurements.uncertainty.(column) :
            fill(0.0, nr)
         for column in value_columns if eb]...,
    )
    Gnuplot.PlotElement(
        xr=[-0.5, nr - 0.5],
        cmds=vcat(["set grid ytics y2tics",
                   "set style data histogram",
                   "set style histogram $hist_style $gap_str",
                   "set boxwidth $box_width",
                   "set style fill $fill_style",
                   """set xtics rotate by $label_rot $(label_rot > 0 ? "right" : "") $(label_enhanced ? "" : "no")enhanced""",
                   ],
                  !isempty(errorbars) ? ["set errorbars $errorbars"] : String[],
                  !isempty(y2cols) ? ["set ytics nomirror", "set y2tics"] : String[],
                  ),
        data=data,
        plot=[ "using $(gpusing(i)) $(axes(group_name))" *
            """title '$(String(group_name))' $(key_enhanced ? "" : "no")enhanced lt $(linetypes[i])"""
               for (i, group_name) in enumerate(group_names)]
    )
end

end


text:
using OrthogonalPolynomialsQuasi, FillArrays, BandedMatrices, ContinuumArrays, QuasiArrays, LazyArrays, FastGaussQuadrature, Test
import OrthogonalPolynomialsQuasi: recurrencecoefficients, basis

@testset "Jacobi" begin
    @testset "JacobiWeight" begin
        a,b = 0.1,0.2
        w = JacobiWeight(a,b)
        @test w.^2 == JacobiWeight(2a,2b)
        @test sqrt.(w) == JacobiWeight(a/2,b/2)
        @test JacobiWeight(0.2,0.3) .* w == JacobiWeight(a+0.2,b+0.3)
        @test LegendreWeight() .* w == w .* LegendreWeight() == w
        @test ChebyshevWeight() .* w == w .* ChebyshevWeight() == JacobiWeight(a-1/2,b-1/2)
    end

    @testset "basics" begin
        @test Legendre() == Jacobi(0,0)
        @test Jacobi(0,0) == Legendre()
    end

    @testset "basis" begin
        b,a = 0.1,0.2
        P = Jacobi(a,b)
        @test P[0.1,2] ≈ 0.16499999999999998
        P = Jacobi(b,a)
        @test P[-0.1,2] ≈ -0.16499999999999998

        A,B,C = recurrencecoefficients(P)
        X = jacobimatrix(P)
        @testset "recurrence coefficient and jacobimatrix" begin
            @test 1/A[1] ≈ X[2,1]
            @test -B[1]/A[1] ≈ X[1,1]
            @test C[2]/A[2] ≈ X[1,2]
            @test 1/A[2] ≈ X[3,2]
            @test -B[2]/A[2] ≈ X[2,2]
            @test C[3]/A[3] ≈ X[2,3]

            @test A[1] ≈ 1/X[2,1]
            @test B[1] ≈ -X[1,1]/X[2,1]
            @test C[2] ≈ X[1,2]/X[3,2]
            @test A[2] ≈ 1/X[3,2]
            @test B[2] ≈ -X[2,2]/X[3,2]
            @test C[3] ≈ X[2,3]/X[4,3]
        end
    end

    @testset "orthogonality" begin
        a,b = 0.1,0.2
        x,w = gaussjacobi(3,a,b)
        P = Jacobi(a,b)

        M = P[x,1:3]'Diagonal(w)*P[x,1:3]
        @test M ≈ Diagonal(M)
        x,w = gaussradau(3,a,b)
        M = P[x,1:3]'Diagonal(w)*P[x,1:3]
        @test M ≈ Diagonal(M)

        w = JacobiWeight(a,b)
        w_2 = sqrt.(w)
        @test w_2 .^ 2 == w
        @test (P' * (w .* P))[1:3,1:3] ≈ (P' * (w .* P))[1:3,1:3] ≈ ((w_2 .* P)'*(w_2 .* P))[1:3,1:3] ≈ M

        @test (Jacobi(0,0)'Jacobi(0,0))[1:10,1:10] ≈ (Legendre()'Legendre())[1:10,1:10]
    end

    @testset "operators" begin
        a,b = 0.1,0.2
        S = Jacobi(a,b)

        @testset "Jacobi" begin
            x = 0.1
            @test S[x,1] === 1.0
            X = jacobimatrix(S)
            @test X[1,1] ≈ (b^2-a^2)/((a+b)*(a+b+2))
            @test X[2,1] ≈ 2/(a+b+2)
            @test S[x,2] ≈ 0.065
            @test S[x,10] ≈ 0.22071099583604945

            w = JacobiWeight(a,b)
            @test w[x] ≈ (1-x)^a * (1+x)^b
            @test OrthogonalPolynomial(w) == S
            wS = w.*S
            @test wS == WeightedJacobi(a,b) == WeightedJacobi{Float64}(a,b)
            @test wS[0.1,1] ≈ w[0.1]
            @test wS[0.1,1:2] ≈ w[0.1] .* S[0.1,1:2]

            w_A = WeightedJacobi(-1/2,0)
            w_B =  WeightedJacobi(1/2,0)

            u = w_A * [1 ; 2; zeros(∞)]
            v = w_B * [1 ; 2; zeros(∞)]
            @test basis(u + v) == w_A
            @test (u+v)[0.1] == u[0.1] + v[0.1]
        end

        @testset "Clenshaw" begin
            x = axes(S,1)
            a = S * (S \ exp.(x))
            A = S \ (a .* S)
            @test (S * (A * (S \ a)))[0.1] ≈ exp(0.2)
            a = S * [1; Zeros(∞)]
            A = S \ (a .* S)
            @test A[1:10,1:10] == I
        end

        @testset "Derivative" begin
            a,b,c = 0.1,0.2,0.3
            S = Jacobi(a,b)
            x = axes(S,1)
            D = Derivative(x)
            u = S \ exp.(x)
            x̃ = 0.1
            @test ((D*S) * u)[x̃] ≈ exp(x̃)
            @test (D * (JacobiWeight(0,b) .* S) * u)[x̃] ≈ exp(x̃) * (1+x̃)^(b-1) * (1+b+x̃)
            @test (D * (JacobiWeight(a,0) .* S) * u)[x̃] ≈ exp(x̃) * (1-x̃)^(a-1) * (1-a-x̃)
            @test (D * (JacobiWeight(a,b) .* S) * u)[x̃] ≈ -exp(x̃)*(1-x̃)^(-1+a)*(1+x̃)^(-1+b)*(a*(1+x̃)+(-1+x̃)*(1+b+x̃))
            @test (D * (JacobiWeight(0,c) .* S) * u)[x̃] ≈ exp(x̃) * (1+x̃)^(c-1) * (1+c+x̃)
            @test (D * (JacobiWeight(c,0) .* S) * u)[x̃] ≈ exp(x̃) * (1-x̃)^(c-1) * (1-c-x̃)
            @test (D * (JacobiWeight(0,0) .* S) * u)[x̃] ≈ exp(x̃)
        end
    end

    @testset "functions" begin
        b,a = 0.1,0.2
        P = Jacobi(a,b)
        D = Derivative(axes(P,1))

        f = P*Vcat(randn(10), Zeros(∞))
        @test (Jacobi(a,b+1) * (Jacobi(a,b+1)\f))[0.1] ≈ f[0.1]
        h = 0.0000001
        @test (D*f)[0.1] ≈ (f[0.1+h]-f[0.1])/h atol=100h

        (D*(JacobiWeight(a,b) .* f))
    end

    @testset "expansions" begin
        P = Jacobi(1/2,0.)
        x = axes(P,1)
        @test (P * (P \ exp.(x)))[0.1] ≈ exp(0.1)

        wP = WeightedJacobi(1/2,0.)
        f = @.(sqrt(1 - x) * exp(x))
        @test wP[0.1,1:100]'*(wP[:,1:100] \ f) ≈ sqrt(1-0.1) * exp(0.1)
        @test (wP * (wP \ f))[0.1] ≈ sqrt(1-0.1) * exp(0.1)

        P̃ = P[affine(Inclusion(0..1), x), :]
        x̃ = axes(P̃, 1)
        @test (P̃ * (P̃ \ exp.(x̃)))[0.1] ≈ exp(0.1)
        wP̃ = wP[affine(Inclusion(0..1), x), :]
        f̃ = @.(sqrt(1 - x̃) * exp(x̃))
        @test wP̃[0.1,1:100]'*(wP̃[:,1:100] \ f̃) ≈ sqrt(1-0.1) * exp(0.1)
        @test (wP̃ * (wP̃ \ f̃))[0.1] ≈ sqrt(1-0.1) * exp(0.1)
    end

    @testset "trivial weight" begin
        S = JacobiWeight(0.0,0.0) .* Jacobi(0.0,0.0)
        @test S == S
        @test Legendre() == S
        @test Legendre()\S isa Eye
    end

    @testset "Jacobi integer" begin
        S = Jacobi(true,true)
        D = Derivative(axes(S,1))
        P = Legendre()

        @test pinv(pinv(S)) === S
        @test P\P === pinv(P)*P === Eye(∞)

        Bi = pinv(Jacobi(2,2))
        @test Bi isa QuasiArrays.PInvQuasiMatrix

        A = Jacobi(2,2) \ (D*S)
        @test typeof(A) == typeof(pinv(Jacobi(2,2))*(D*S))
        @test A isa BandedMatrix
        @test bandwidths(A) == (-1,1)
        @test size(A) == (∞,∞)
        @test A[1:10,1:10] == diagm(1 => 2:0.5:6)

        M = @inferred(D*S)
        @test M isa MulQuasiMatrix
        @test M.args[1] == Jacobi(2,2)
        @test M.args[2][1:10,1:10] == A[1:10,1:10]
    end

    @testset "Weighted Jacobi integer" begin
        S = Jacobi(true,true)
        w̃ = JacobiWeight(false,true)
        A = Jacobi(true,false)\(w̃ .* S)
        @test A isa BandedMatrix
        @test size(A) == (∞,∞)
        @test A[1:10,1:10] ≈ (Jacobi(1.0,0.0) \ (JacobiWeight(0.0,1.0) .* Jacobi(1.0,1.0)))[1:10,1:10]

        w̄ = JacobiWeight(true,false)
        A = Jacobi(false,true)\(w̄.*S)
        @test A isa BandedMatrix
        @test size(A) == (∞,∞)
        @test A[1:10,1:10] ≈ (Jacobi(0.0,1.0) \ (JacobiWeight(1.0,0.0) .* Jacobi(1.0,1.0)))[1:10,1:10]

        P = Legendre()
        w̄ = JacobiWeight(true,false)
        @test_broken P \ (w̃ .* Jacobi(true,false))
        w̄ = JacobiWeight(false,true)
        @test (P \ (w̃ .* Jacobi(false,true)))[1:10,1:10] == diagm(0 => ones(10), -1 => ones(9))

        w = JacobiWeight(true,true)
        A,B = (P'P),P\(w.*S)

        M = Mul(A,B)
        @test M[1,1] == 4/3

        M = ApplyMatrix{Float64}(*,A,B)
        M̃ = M[1:10,1:10]
        @test M̃ isa BandedMatrix
        @test bandwidths(M̃) == (2,0)

        @test A*B isa MulMatrix
        @test bandwidths(A*B) == bandwidths(B)

        A,B,C = (P\(w.*S))',(P'P),P\(w.*S)
        M = ApplyArray(*,A,B,C)
        @test bandwidths(M) == (2,2)
        @test M[1,1] ≈  1+1/15
        M = A*B*C
        @test bandwidths(M) == (2,2)
        @test M[1,1] ≈  1+1/15

        S = Jacobi(1.0,1.0)
        w = JacobiWeight(1.0,1.0)
        wS = w .* S

        W = Diagonal(w)
        @test W[0.1,0.2] ≈ 0.0
    end

    @testset "Jacobi and Chebyshev" begin
        T = ChebyshevT()
        U = ChebyshevU()
        JT = Jacobi(T)
        JU = Jacobi(U)

        @testset "recurrence degenerecies" begin
            A,B,C = recurrencecoefficients(JT)
            @test A[1] == 0.5
            @test B[1] == 0.0
        end

        @test JT[0.1,1:4] ≈ [1.0,0.05,-0.3675,-0.0925]

        @test ((T \ JT) * (JT \ T))[1:10,1:10] ≈ Eye(10)
        @test ((U \ JU) * (JU \ U))[1:10,1:10] ≈ Eye(10)

        @test T[0.1,1:10]' ≈ JT[0.1,1:10]' * (JT \ T)[1:10,1:10]
        @test U[0.1,1:10]' ≈ JU[0.1,1:10]' * (JU \ U)[1:10,1:10]

        w = JacobiWeight(1,1)
        @test (w .* U)[0.1,1:10] ≈ (T * (T \ (w .* U)))[0.1,1:10]
        @test (w .* U)[0.1,1:10] ≈ (JT * (JT \ (w .* U)))[0.1,1:10]
        @test (w .* JU)[0.1,1:10] ≈ (T * (T \ (w .* JU)))[0.1,1:10]
        @test (w .* JU)[0.1,1:10] ≈ (JT * (JT \ (w .* JU)))[0.1,1:10]
    end

    @testset "Jacobi-Chebyshev-Ultraspherical transforms" begin
        @test Jacobi(0.0,0.0) \ Legendre() == Eye(∞)
        @test ((Ultraspherical(3/2) \ Jacobi(1,1))*(Jacobi(1,1) \ Ultraspherical(3/2)))[1:10,1:10] ≈ Eye(10)
        f = Jacobi(0.0,0.0)*[[1,2,3]; zeros(∞)]
        g = (Legendre() \ f) - f.args[2]
        @test_skip norm(g) ≤ 1E-15
        @test_broken (Legendre() \ f) == f.args[2]
        @test (Legendre() \ f)[1:10] ≈ f.args[2][1:10]
        f = Jacobi(1.0,1.0)*[[1,2,3]; zeros(∞)]
        g = Ultraspherical(3/2)*(Ultraspherical(3/2)\f)
        @test f[0.1] ≈ g[0.1]

        @testset "Chebyshev-Legendre" begin
            T = Chebyshev()
            P = Legendre()
            @test T[:,Base.OneTo(5)] \ P[:,Base.OneTo(5)] == (T\P)[1:5,1:5]

            x = axes(P,1)
            u = P * (P \ exp.(x))
            @test u[0.1] ≈ exp(0.1)

            P = Legendre{BigFloat}()
            x = axes(P,1)
            u = P * (P \ exp.(x))
            @test u[BigFloat(1)/10] ≈ exp(BigFloat(1)/10)
        end
    end

    @testset "hcat" begin
        L = LinearSpline(range(-1,1;length=2))
        S = JacobiWeight(1.0,1.0) .* Jacobi(1.0,1.0)
        P = apply(hcat,L,S)
        @test P isa ApplyQuasiArray
        @test axes(P) == axes(S)
        V = view(P,0.1,1:10)
        @test all(arguments(V) .≈ [L[0.1,:], S[0.1,1:8]])
        @test P[0.1,1:10] == [L[0.1,:]; S[0.1,1:8]]
        D = Derivative(axes(P,1))
        # applied(*,D,P) |> typeof
        # MemoryLayout(typeof(D))
    end

    @testset "Jacobi Clenshaw" begin
        P = Jacobi(0.1,0.2)
        x = axes(P,1)
        a = P * (P \ exp.(x))
        @test a[0.1] ≈ exp(0.1)
        M = P \ (a .* P);
        u = [randn(1000); zeros(∞)];
        @test (P * (M*u))[0.1] ≈ (P*u)[0.1]*exp(0.1)
    end

    @testset "mapped" begin
        R = jacobi(0,1/2,0..1) \ jacobi(0,-1/2,0..1)
        R̃ = Jacobi(0,1/2) \ Jacobi(0,-1/2)
        @test R[1:10,1:10] == R̃[1:10,1:10]
    end

    @testset "Christoffel–Darboux" begin
        a,b = 0.1,0.2
        P = Jacobi(a,b)
        X = P\ (axes(P,1) .* P)
        Mi = inv(P'*(JacobiWeight(a,b) .* P))
        x,y = 0.1,0.2
        n = 10
        Pn = Diagonal([Ones(n); Zeros(∞)])
        Min = Pn * Mi
        @test norm((X*Min - Min*X')[1:n,1:n]) ≤ 1E-13
        β = X[n,n+1]*Mi[n+1,n+1]
        @test (x-y) * P[x,1:n]'Mi[1:n,1:n]*P[y,1:n] ≈ P[x,n:n+1]' * (X*Min - Min*X')[n:n+1,n:n+1] * P[y,n:n+1] ≈ P[x,n:n+1]' * [0 -β; β 0] * P[y,n:n+1]

        @testset "extrapolation" begin
            x,y = 0.1,3.4
            @test (x-y) * P[x,1:n]'Mi[1:n,1:n]*Base.unsafe_getindex(P,y,1:n) ≈ P[x,n:n+1]' * [0 -β; β 0] * Base.unsafe_getindex(P,y,n:n+1)
        end
    end
end


text:
using RDatasets, DataFrames, Gadfly

set_default_plot_size(5inch, 5inch/golden);

plot(dataset("datasets","iris"), x="SepalWidth", y="SepalLength", color="Species")

max_stars_repo_path: Section 4/Julia+Iris.jl

max_issues_repo_path: Section 4/Julia+Iris.jl

max_forks_repo_path: Section 4/Julia+Iris.jl


text:
# ---
# title: 85. Maximal Rectangle
# id: problem85
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Hash Table, Dynamic Programming, Stack
# link: <https://leetcode.com/problems/maximal-rectangle/description/>
# hidden: true
# ---
#
# Given a `rows x cols` binary `matrix` filled with `0`'s and `1`'s, find the
# largest rectangle containing only `1`'s and return _its area_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/14/maximal.jpg)
#
#
#
#     Input: matrix = [["1","0","1","0","0"],["1","0","1","1","1"],["1","1","1","1","1"],["1","0","0","1","0"]]
#     Output: 6
#     Explanation: The maximal rectangle is shown in the above picture.
#
#
# **Example 2:**
#
#
#
#     Input: matrix = []
#     Output: 0
#
#
# **Example 3:**
#
#
#
#     Input: matrix = [["0"]]
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: matrix = [["1"]]
#     Output: 1
#
#
# **Example 5:**
#
#
#
#     Input: matrix = [["0","0"]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `rows == matrix.length`
#   * `cols == matrix.length`
#   * `0 <= row, cols <= 200`
#   * `matrix[i][j]` is `'0'` or `'1'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# TODO: add check to ensure minimum seg value is 1
inorder(x, y) = (min(x, y), max(x, y))
interleave(x, y) = map((x, y) -> (x, y), x, y)

function collapse{T}(x::Array{T, 3})
    # Collapse a 3 dimensional array to a 2 dimensional array of vectors
    nrows, ncols, depth = size(x)
    result = [T[] for i=1:nrows, j=1:ncols]
    for k = 1:depth
        for j = 1:ncols
            for i = 1:nrows
                push!(result[i, j], x[i, j, k])
            end
        end
    end

    result
end

function neighborpairs{T}(spmap::Matrix{T})
    neighbordict = Dict{Int,Set{Int}}()
    for j = 1:size(spmap, 2)
        for i = 1:size(spmap, 1)
            if j < size(spmap, 2)
                left = spmap[i,j]
                right = spmap[i,j+1]
                if left != right
                    minind = min(left, right)
                    maxind = max(left, right)

                    if !haskey(neighbordict, minind)
                        neighbordict[minind] = Set{Int}()
                    end

                    push!(neighbordict[minind], maxind)
                end
            end
            if i < size(spmap, 1)
                top = spmap[i,j]
                bottom = spmap[i+1,j]
                if top != bottom
                    minind = min(top, bottom)
                    maxind = max(top, bottom)

                    if !haskey(neighbordict, minind)
                        neighbordict[minind] = Set{Int}()
                    end

                    push!(neighbordict[minind], maxind)
                end
            end
        end
    end

    rval = NTuple{2,Int}[]
    for key in sort(collect(keys(neighbordict)))
        for val in sort(collect(neighbordict[key]))
            push!(rval, (key, val))
        end
    end
    rval
end

function spdist{T<:Integer}(segmap::Matrix{T}, sp1indices::Vector{Int}, sp2indices::Vector{Int})
    # Return the distance between the pixels in sp1indices and sp2indices, as defined by segmap.
    sp1seg = segmap[sp1indices]
    sp2seg = segmap[sp2indices]

    maxsegval = maximum([maximum(sp1seg), maximum(sp2seg)])

    _, sp1segfreq = hist(sp1seg, 0:maxsegval)
    _, sp2segfreq = hist(sp2seg, 0:maxsegval)

    acc = 0
    for x in unique(sp1seg)
        acc += sp1segfreq[x] * sp2segfreq[x]
    end

    return 1 - (acc / (length(sp1indices) * length(sp2indices)))
end

function spdist{T<:Integer}(spmap::Matrix{T}, segmap::Matrix{T})
    spindices = Dict{Int, Vector{Int}}()
    for (index, val) in enumerate(spmap)
        if !haskey(spindices, val)
            spindices[val] = [index]
        else
            push!(spindices[val], index)
        end
    end

    ((T, T) => Float64)[(i, j) => spdist(segmap, spindices[i], spindices[j]) for (i, j) in neighborpairs(spmap)]
end

function seg2boundary{T<:Integer}(seg::Matrix{T}; format=:imagesize)
    # Ported by Jake Snell <jsnell@cs.toronto.edu> from MATLAB code due to Pablo Arbelaez <arbelaez@eecs.berkeley.edu>
    (tx, ty) = size(seg)

    boundary = convert(BitArray, zeros(2*tx+1, 2*ty+1))
    hedges = hcat(!(seg[:, 1:end-1] .== seg[:, 2:end]), zeros(Bool, (tx,1)))
    vedges = vcat(!(seg[1:end-1, :] .== seg[2:end, :]), zeros(Bool, (1,ty)))

    boundary[3:2:end, 2:2:end] = vedges
    boundary[2:2:end, 3:2:end] = hedges
    boundary[3:2:end-1, 3:2:end-1] = hedges[1:end-1, 1:end-1] | hedges[2:end, 1:end-1] |
                                     vedges[1:end-1, 1:end-1] | vedges[1:end-1, 2:end]

    boundary[1, :] = boundary[2, :]
    boundary[:, 1] = boundary[:, 2]
    boundary[end, :] = boundary[end-1, :]
    boundary[:, end] = boundary[:, end-1]

    if format == :imagesize
        boundary = boundary[3:2:end, 3:2:end]
    end
    boundary
end

function seg2boundaryindex{T<:Integer}(seg::Matrix{T})
    (tx, ty) = size(seg)

    hedges = hcat(interleave(seg[:, 1:end-1], seg[:, 2:end]), fill((0,0), (tx, 1)))
    vedges = vcat(interleave(seg[1:end-1, :], seg[2:end, :]), fill((0,0), (1, ty)))

    ul = collapse(cat(3, hedges[1:end-1, 1:end-1], hedges[2:end, 1:end-1],
    vedges[1:end-1, 1:end-1], vedges[1:end-1, 2:end]))
    ur = collapse(reshape(vedges[1:end-1, end], (tx-1, 1, 1)))
    bl = collapse(reshape(hedges[end, 1:end-1], (1, ty-1, 1)))
    br = collapse(cat(3, hedges[end, end], vedges[end, end]))

    unfiltered = [ul ur; bl br]

    Vector{(Int, Int)}[unique(map(ys -> inorder(ys[1], ys[2]), filter(zs -> zs[1] != zs[2], unfiltered[i, j]))) for i = 1:tx, j = 1:ty]
end

function seg2boundarystrength{T1<:Integer,T2<:FloatingPoint}(spmat::Matrix{T1}, strengths::Dict{NTuple{2,T1},T2})
    nrows, ncols = size(spmat)

    rval = zeros(T2, 2nrows+1, 2ncols+1)
    hedges = hcat(!(spmat[:, 1:end-1] .== spmat[:, 2:end]), zeros(T2, nrows, 1))
    vedges = vcat(!(spmat[1:end-1, :] .== spmat[2:end, :]), zeros(T2, 1, ncols))

    hedges = zeros(T2, nrows, ncols)
    for j = 1:ncols-1
        for i = 1:nrows
            if spmat[i, j] != spmat[i, j+1]
                hedges[i, j] = strengths[spmat[i, j], spmat[i, j+1]]
            end
        end
    end

    vedges = zeros(T2, nrows, ncols)
    for j = 1:ncols
        for i = 1:nrows-1
            if spmat[i, j] != spmat[i+1, j]
                vedges[i, j] = strengths[spmat[i, j], spmat[i+1, j]]
            end
        end
    end

    rval[3:2:end, 2:2:end] = vedges
    rval[2:2:end, 3:2:end] = hedges
    rval[3:2:end-1, 3:2:end-1] = max(hedges[1:end-1, 1:end-1],
                                     hedges[2:end, 1:end-1],
                                     vedges[1:end-1, 1:end-1],
                                     vedges[1:end-1, 2:end])
    rval[1, :] = rval[2, :]
    rval[:, 1] = rval[:, 2]
    rval[end, :] = rval[end-1, :]
    rval[:, end] = rval[:, end-1]

    rval[3:2:end, 3:2:end]
end


text:
module Jig

export @runtest,
       @context,
       jtest

macro runtest(pkg, file...)
  for f in file
    print_with_color(:blue, string(f))
    print_with_color(:blue, ":")
    println("")
    include(Pkg.dir("$pkg/test/$f.jl"))
  end
end

macro context(s)
  print("      ")
  print_with_color(:magenta,s)
  print("  ")
end

function jtest(ex...)
  fails = 0
  for e in ex
    if eval(:($e))
      print_with_color(:green, ".")
    else
      print_with_color(:red, "x")
      fails += 1
    end
  end
  println("")
  fails
end

end


text:
using Base.Test
using Requests: json
using HavenOnDemand

try
    set_api_key(ENV["HOD_API_KEY"])
catch ex
    if isa(ex, KeyError)
        info("You have to setup a `HOD_API_KEY` environment variable first.")
    end
end

# `analyze_sentiment`
@test_throws HODException analyze_sentiment()

try
    # needs at least one valid kwarg
    analyze_sentiment(non = :valid, keyword = :arguments)
catch ex
    if isa(ex, HODException)
        response = json(ex.response)
        @test response["error"]  == 4015
        @test response["reason"] == "Missing required parameter(s)"
    end
end

positive_response = analyze_sentiment(text = "I love all zombies!")
@test isa(positive_response, Dict{AbstractString, Any})

@test positive_response["aggregate"]["sentiment"] == "positive"
@test isempty(positive_response["negative"])
@test positive_response["positive"][1]["topic"] == "all zombies"
@test positive_response["positive"][1]["score"] > 0

negative_response = analyze_sentiment(
    text    = "I hate all plants!",    # ...some zombie
    this_is = :ignored    # only non valid kwargs are ignored
)

@test negative_response["aggregate"]["sentiment"] == "negative"
@test isempty(negative_response["positive"])
@test negative_response["negative"][1]["topic"] == "all plants"
@test negative_response["negative"][1]["score"] < 0


text:
using PyCall

@pyimport modeller
@pyimport modeller.scripts as scripts
@pyimport _modeller

mod_lib = _modeller.mod_libdir_get()

modeller.log[:verbose]()    # request verbose output
env = modeller.environ()
env[:libs][:topology][:read](file=string(mod_lib, "/top_heav.lib")) # read topology
env[:libs][:parameters][:read](file=string(mod_lib, "/par.lib")) # read parameters

# read model file
mdl = scripts.complete_pdb(env, "TvLDH.B99990002.pdb")

# Assess with DOPE:
s = modeller.selection(mdl)   # all atom selection
s[:assess_dope](output="ENERGY_PROFILE NO_REPORT", file="TvLDH.profile",
              normalize_profile=true, smoothing_window=15)


text:
export
    opt  # least square solver

"""
    solve the least square problem

arguments:

  x,
  u,
  B,

"""
function opt(x, u, B)
  set_default_solver(SCSSolver(verbose=0))
  A = x'
  b = u'
  m, n = size(A)

  w = Variable(n)
  problem = minimize(sumsquares(A * w - b), [sumsquares(w) <= B^2])

  solve!(problem)

  problem.status == :Optimal && print(".")

  w.value'
end

max_stars_repo_name: innerlee/Halfspace

max_issues_repo_name: innerlee/Halfspace

max_forks_repo_name: innerlee/Halfspace


text:
using MPI
using Test
using StaticArrays
using Logging, Printf

using CLIMA
using CLIMA.LinearSolvers
using CLIMA.GeneralizedConjugateResidualSolver
using CLIMA.GeneralizedMinimalResidualSolver
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.MPIStateArrays
using CLIMA.VariableTemplates
using CLIMA.DGmethods
import CLIMA.DGmethods: BalanceLaw, vars_aux, vars_state, vars_gradient,
                        vars_diffusive, flux_nondiffusive!, flux_diffusive!,
                        source!, boundary_state!,
                        gradvariables!,
                        diffusive!, init_aux!, init_state!,
                        LocalGeometry

import CLIMA.DGmethods.NumericalFluxes: NumericalFluxDiffusive,
                                        numerical_flux_diffusive!

@static if haspkg("CuArrays")
  using CUDAdrv
  using CUDAnative
  using CuArrays
  CuArrays.allowscalar(false)
  const ArrayTypes = (CuArray, )
else
  const ArrayTypes = (Array, )
end

if !@isdefined integration_testing
  const integration_testing =
    parse(Bool, lowercase(get(ENV,"JULIA_CLIMA_INTEGRATION_TESTING","false")))
end

struct PoissonModel{dim} <: BalanceLaw end

vars_aux(::PoissonModel,T) = @vars(rhs_ϕ::T)
vars_state(::PoissonModel, T) = @vars(ϕ::T)
vars_gradient(::PoissonModel, T) = @vars(ϕ::T)
vars_diffusive(::PoissonModel, T) = @vars(∇ϕ::SVector{3, T})

boundary_state!(nf, bl::PoissonModel, _...) = nothing

function flux_nondiffusive!(::PoissonModel, flux::Grad, state::Vars,
                            auxstate::Vars, t::Real)
  nothing
end

function flux_diffusive!(::PoissonModel, flux::Grad, state::Vars,
                         diffusive::Vars, auxstate::Vars, t::Real)
  flux.ϕ = diffusive.∇ϕ
end

struct PenaltyNumFluxDiffusive <: NumericalFluxDiffusive end
function numerical_flux_diffusive!(::PenaltyNumFluxDiffusive,
                                   bl::PoissonModel, F::MArray, nM,
                                   QM, QVM, auxM,
                                   QP, QVP, auxP,
                                   t)
  numerical_flux_diffusive!(CentralNumericalFluxDiffusive(),
                            bl, F, nM,
                            QM, QVM, auxM,
                            QP, QVP, auxP,
                            t)
  FT = eltype(F)
  nstate = length(F)
  tau = FT(1.0)
  for s = 1:nstate
    @inbounds F[s] -= tau * (QM[1] - QP[1])
  end
end

function gradvariables!(::PoissonModel, transformstate::Vars, state::Vars, auxstate::Vars, t::Real)
  transformstate.ϕ = state.ϕ
end

function diffusive!(::PoissonModel, diffusive::Vars,
                    ∇transform::Grad, state::Vars, auxstate::Vars, t::Real)
  diffusive.∇ϕ = ∇transform.ϕ
end

source!(::PoissonModel, source::Vars, state::Vars, aux::Vars, t::Real) = nothing

# note, that the code assumes solutions with zero mean
sol1d(x) = sin(2pi * x) ^ 4 - 3 / 8
dxx_sol1d(x) = -16pi ^ 2 * sin(2pi * x) ^ 2 * (sin(2pi * x) ^ 2 - 3cos(2pi * x) ^ 2)

function init_aux!(::PoissonModel{dim}, aux::Vars, g::LocalGeometry) where dim
  aux.rhs_ϕ = 0
  @inbounds for d = 1:dim
    x1 = g.coord[d]
    x2 = g.coord[1 + mod(d, dim)]
    x3 = g.coord[1 + mod(d + 1, dim)]
    x23 = SVector(x2, x3)
    aux.rhs_ϕ -= dxx_sol1d(x1) * prod(sol1d, view(x23, 1:dim-1))
  end
end

function init_state!(::PoissonModel{dim}, state::Vars, aux::Vars, coords, t) where dim
  state.ϕ = prod(sol1d, view(coords, 1:dim))
end

function run(mpicomm, ArrayType, FT, dim, polynomialorder, brickrange, periodicity, linmethod)

  topology = BrickTopology(mpicomm, brickrange, periodicity=periodicity)
  grid = DiscontinuousSpectralElementGrid(topology,
                                          polynomialorder = polynomialorder,
                                          FloatType = FT,
                                          DeviceArray = ArrayType)
  dg = DGModel(PoissonModel{dim}(),
               grid,
               CentralNumericalFluxNonDiffusive(),
               PenaltyNumFluxDiffusive(),
               CentralGradPenalty())

  Q = init_ode_state(dg, FT(0))
  Qrhs = dg.auxstate
  Qexact = init_ode_state(dg, FT(0))

  linearoperator!(y, x) = dg(y, x, nothing, 0; increment = false)

  linearsolver = linmethod(Q)

  iters = linearsolve!(linearoperator!, Q, Qrhs, linearsolver)

  error = euclidean_distance(Q, Qexact)

  @info @sprintf """Finished
  error = %.16e
  iters = %d
  """ error iters
  error, iters
end

let
  MPI.Initialized() || MPI.Init()
  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))
  @static if haspkg("CUDAnative")
    device!(MPI.Comm_rank(mpicomm) % length(devices()))
  end

  polynomialorder = 4
  base_num_elem = 4
  tol = 1e-9

  linmethods = (b -> GeneralizedConjugateResidual(3, b, tol),
                b -> GeneralizedMinimalResidual(7, b, tol)
               )

  expected_result = Array{Float64}(undef, 2, 2, 3) # method, dim-1, lvl

  # GeneralizedConjugateResidual
  expected_result[1, 1, 1] = 5.0540243616448058e-02
  expected_result[1, 1, 2] = 1.4802275366044011e-03
  expected_result[1, 1, 3] = 3.3852821775121401e-05
  expected_result[1, 2, 1] = 1.4957957657736219e-02
  expected_result[1, 2, 2] = 4.7282369781541172e-04
  expected_result[1, 2, 3] = 1.4697449643351771e-05

  # GeneralizedMinimalResidual
  expected_result[2, 1, 1] = 5.0540243587512981e-02
  expected_result[2, 1, 2] = 1.4802275409186211e-03
  expected_result[2, 1, 3] = 3.3852820667079927e-05
  expected_result[2, 2, 1] = 1.4957957659220951e-02
  expected_result[2, 2, 2] = 4.7282369895963614e-04
  expected_result[2, 2, 3] = 1.4697449516628483e-05

  lvls = integration_testing ? size(expected_result)[end] : 1

  for ArrayType in ArrayTypes, (m, linmethod) in enumerate(linmethods), FT in (Float64,)
    result = Array{Tuple{FT, Int}}(undef, lvls)
    for dim = 2:3

      for l = 1:lvls
        Ne = ntuple(d -> 2 ^ (l - 1) * base_num_elem, dim)
        brickrange = ntuple(d -> range(FT(0), length = Ne[d], stop = 1), dim)
        periodicity = ntuple(d -> true, dim)

        @info (ArrayType, FT, m, dim)
        result[l] = run(mpicomm, ArrayType, FT, dim,
                        polynomialorder, brickrange, periodicity, linmethod)

        @test isapprox(result[l][1], FT(expected_result[m, dim-1, l]), rtol = sqrt(tol))
      end

      if integration_testing
        @info begin
          msg = ""
          for l = 1:lvls-1
            rate = log2(result[l][1]) - log2(result[l + 1][1])
            msg *= @sprintf("\n  rate for level %d = %e\n", l, rate)
          end
          msg
        end
      end
    end
  end
end

nothing

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
const liquids = Vector{Symbol}()

const DefaultLiquid = Liquid()
const Water = Liquid()
export DefaultLiquid, Water
push!(liquids,:Water)

@liquid Seawater 15.6u"°C" 1030 1.20e-3  7.34e-2 1.77e3  2.34e9
@liquid Glycerin 20u"°C" 1260  1.50e0  6.33e-2 1.4e-2 4.52e9
@liquid Gasoline 15.6u"°C" 680  3.1e-4 2.2e-2  5.5e4 1.3e9

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


text:
using HardSphereDynamics
include("makie_visualization.jl")

fluid = HardSphereFluid(1, 10, 0.03)  # create hard spheres in 1D
positions, velocities, times = evolve!(fluid, δt, final_time)

positions2 =  HardSphereDynamics.to_3D(positions)
velocities2 = HardSphereDynamics.to_3D(velocities)

visualize_3d(positions2, velocities2, [ball.r for ball in fluid.balls], sleep_step=0.01)

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
include("Yfind.jl")

# setprecision(100)
# Q=.2
# P=.3
#
#
# Energy=.25
# H=(2*Energy)^2
#
# Y=Yfind(Q,P,H)
#
# u0=zeros(4)
#
# u0[1]=0 #X
# u0[2]=P #P
# u0[3]=Q #Q
# u0[4]=Y #Y
#
# TEST1_0=H_test(u0)-H

Q=BigFloat(".2")
P=BigFloat(".3")

Energy=BigFloat(".25")
H=(2*Energy)^2

Y=Yfind(Q,P,H)

u0=zeros(BigFloat ,4)

u0[1]=0 #X
u0[2]=P #P
u0[3]=Q #Q
u0[4]=Y #Y

TEST1_0=H_test(u0)-H


text:
const Theme = Scene

# this should probably be part of @default or actually replace it.
# but a seperate macro for a node conversion seemed easier right now
# So this just works like @default, but for individual nodes with concrete values -
# it's there to build themes etc.

function matchnode(assignment)
    @capture(assignment, a_ = b_) || error("needs to be an assignment. Found: $assignment")
    if isa(b, Expr)
        if b.head == :block
            rhs = []
            for elem in b.args
                Base.is_linenumber(elem) && continue
                expr, name = matchnode(elem)
                push!(rhs, :($(QuoteNode(name)) => $expr))
            end
            expr = quote
                $(esc(a)) = Scene($(rhs...))
                $(esc(a))
            end
            return expr, a
        else
            if @capture(b,
                    (f_(args__)) |
                    (args__::f_)
                )
                return :(to_node($(esc.(args)...), x-> $(esc(f))(Makie.current_backend[], x))), a
            end
        end
    end
    return esc(b), a
end

macro theme(assignment)
    expr, name = matchnode(assignment)
    expr
end

# BIG TODO: make themes + defaults + conversion function @default + @theme macro suck less!
# Right now there is a lot of annoying and unecessary dublication going on

function default_theme(scene)
    q1 = qrotation(Vec3f0(1, 0, 0), -0.5f0*pi)
    q2 = qrotation(Vec3f0(0, 0, 1), 1f0*pi)
    tickrotations = (
        qrotation(Vec3f0(0,0,1), -1.5pi),
        q2,
        qmul(qmul(q2, q1), qrotation(Vec3f0(0, 1, 0), 1pi))
    )

    tickalign = (
        (:hcenter, :left), # x axis
        (:right, :vcenter), # y axis
        (:right, :vcenter), # z axis
    )
    dark_text = RGBAf0(0.0, 0.0, 0.0, 0.4)

    axisnames = map(x-> ("$x Axis", 0.1, dark_text, Vec4f0(0,0,0,1), (:center, :bottom)), (:X, :Y, :Z))
    showticks = ntuple(i-> true, 3)
    tickfont3d = ntuple(i-> (0.1, RGBAf0(0.5, 0.5, 0.5, 0.6), tickrotations[i], tickalign[i]), 3)
    tickfont2d = ntuple(i-> (0.1, RGBAf0(0.5, 0.5, 0.5, 0.6), Vec4f0(0,0,0,1), tickalign[i]), 2)
    showaxis = ntuple(i-> true, 3)
    showgrid = ntuple(i-> true, 3)

    scalefuncs = ntuple(i-> identity, 3)
    gridthickness = ntuple(x-> 1f0, 3)
    colors = UniqueColorIter(:Set1)
    meshrotation = Vec4f0(0, 0, 0, 1)
    light = Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), Vec3f0(20,20,20)]
    @theme theme = begin

        rotation = to_rotation(Vec4f0(0, 0, 0, 1))
        scale = to_scale(Vec3f0(1))
        offset = to_offset(Vec3f0(0))

        camera = to_camera(:auto)
        visible = to_bool(true)
        show = to_bool(true)

        light = to_static_vec(light)

        drawover = to_bool(false)

        color = colors
        linewidth = to_float(1)
        colormap = to_colormap(:YlGnBu)
        colornorm = nothing # nothing for calculating it from intensity

        surface = begin
            image = nothing
        end

        contour = begin
            levels = to_float(5)
            fillrange = false
        end

        scatter = begin
            marker = to_spritemarker(Circle)
            markersize = to_markersize2d(0.1)
            strokecolor = to_color(RGBA(0, 0, 0, 0))
            strokewidth = to_float(0)
            glowcolor = to_color(RGBA(0, 0, 0, 0))
            glowwidth = to_float(0)
            rotations = to_rotations(Billboard())
        end

        meshscatter = begin
            marker = to_mesh(Sphere(Point3f0(0), 0.1f0))
            markersize = to_markersize3d(1)
            rotations = meshrotation
        end

        lines = begin
            linestyle = to_linestyle(nothing)
        end

        image = begin
            spatialorder = :yx
        end

        mesh = begin
            shading = true
            attribute_id = nothing
            indices = nothing
        end

        axis = begin
            axisnames = axisnames
            visible = true

            showticks = showticks
            tickfont2d = tickfont2d
            tickfont3d = tickfont3d
            showaxis = showaxis
            showgrid = showgrid

            scalefuncs = scalefuncs
            gridcolors = to_color(ntuple(x-> RGBAf0(0.5, 0.5, 0.5, 0.4), 3))
            gridthickness = gridthickness
            axiscolors = to_color(ntuple(x-> dark_text, 3))
        end

        heatmap = begin
            linewidth = to_float(0)
            levels = to_float(1f0)
        end
        volume = begin
            algorithm = to_volume_algorithm(GLVisualize.MaximumIntensityProjection)
            absorption = to_float(1f0)
            isovalue = to_float(0.5f0)
            isorange = to_float(0.01f0)
        end

        legend = begin
            backgroundcolor = to_color(:white)
            strokecolor = to_color(RGBA(0.3, 0.3, 0.3, 0.9))
            strokewidth = to_float(1)
            position = to_position((0, 1))
            gap = to_float(20)
            textgap = to_float(15)
            labelwidth = to_float(20)
            padding = to_float(10)
            outerpadding = to_markersize2d(10)
            align = to_textalign((:left, :hcenter))
            rotation = to_rotation(Vec4f0(0, 0, 0, 1))
            textcolor = to_color(:black)
            textsize = to_float(16)
            markersize = to_markersize2d(5)
            linepattern = to_positions(Point2f0[(0, 0), (1, 0.0)])
            scatterpattern = to_positions(Point2f0[(0.5, 0.0)])
        end

        color_legend = begin
            width = to_markersize2d((20, 200))
            backgroundcolor = to_color(:white)
            strokecolor = to_color(RGBA(0.3, 0.3, 0.3, 0.9))
            strokewidth = to_float(1)
            position = to_position((1, 1))
            textgap = to_float(15)
            padding = to_markersize2d(10)
            outerpadding = to_markersize2d(10)
            align = Makie.to_textalign((:left, :hcenter))
            rotation = to_rotation(Vec4f0(0, 0, 0, 1))
            textcolor = to_color(:black)
            textsize = to_float(16)
        end
    end
    scene[:theme] = theme
end


text:
using Pkg; Pkg.activate(".")
using Temporal

# define universe and gather data
assets = ["DOGE-BUSD"]
function datasource(asset::String)::TS
    savedata_path = joinpath("data", "$asset.csv")
    return Temporal.tsread(savedata_path, indextype=UInt64, format="yyyy-mm-dd HH:MM:SS")

end

df = datasource(assets[end])

save_path(asset::String) = joinpath("/home/bane/projects/julia/Temporal.jl", "data", "$asset.csv")

file = save_path(assets[1])
eol = '\n'
header = true
dlm = ','
indextype = UInt64
format="yyyy-mm-dd"

csv = Vector{String}(split(read(file, String), eol))
if csv[end] == ""
    pop!(csv)  # remove final blank line
end
if header
    fields = Vector{String}(split(popfirst!(csv), dlm)[2:end])
    k = length(fields)
    n = length(csv)
else
    # k = matchcount(csv[1], dlm)
    k = sum([csv[1][i] == dlm for i in 1:length(csv[1])])
    n = length(csv)
    fields = autocol(1:k)
end
# Fill data
arr = zeros(Float64, (n,k))
if indextype != UInt64
    idx = fill("", n)::Vector{String}
else
    idx = zeros(UInt64, n)::Vector{UInt64}
end
for i = 1:n
    s = Vector{String}(split(csv[i], dlm))
    idx[i] = indextype != UInt64 ? popfirst!(s) : parse(UInt64, popfirst!(s))
    s[s.==""] .= "NaN"
    for j in 1:length(s)
        arr[i,j] = parse(Float64, s[j])
    end
end
    return TS(arr, indextype.(idx, format), fields)


text:
#
# files.jl --
#
# Implement reading/writing of OI-FITS data from/to FITS files.
#
#------------------------------------------------------------------------------

# Read a column from an OI-FITS table.
function read_column(ff::FITSFile, colnum::Integer, multiplier::Integer)
    # Minimal check.
    fits_assert_open(ff)

    # Get the type and the dimensions of the data stored in the column.
    (eqcoltype, repcnt, width) = fits_get_eqcoltype(ff, colnum)
    dims = fits_read_tdim(ff, colnum)
    nrows = fits_get_num_rows(ff)

    # Allocate the array and read the column contents.
    T = eqcoltype_to_type(eqcoltype) # FIXME: improve type stability
    if T <: String
        # Column contains an array of strings.  Strip the leading dimension
        # which is the maximum length of each strings.  On return trailing
        # spaces are removed (they are insignificant according to the FITS
        # norm) and the result converted to a regular string.
        if length(dims) == 1
            dims = nrows
        else
            dims[1:end-1] = dims[2:end]
            dims[end] = nrows
        end
        data = Array{T}(undef, dims...)
        fits_read_col(ff, colnum, 1, 1, data)
        return map(str -> String(rstrip(c -> c == ' ', str)), data)
    elseif T === Nothing
        error("unsupported column data")
    else
        # Column contains numerical data.
        if length(dims) == 1 && dims[1] == multiplier == 1
            # Result will be a simple vector.
            resize!(dims, 1)
            dims[1] = nrows
        else
            # Result will be a multi-dimensional array.
            push!(dims, nrows)
        end
        data = Array{T}(undef, dims...)
        fits_read_col(ff, colnum, 1, 1, data)
        return data
    end
end

function get_file_handle(hdu::HDU)
    fh = hdu.fitsfile
    fits_assert_open(fh)
    fits_movabs_hdu(fh, hdu.ext)
    return fh
end

"""
    OIFITS.read_table(hdu, sel=everything)

reads the FITS table in FITS Header Data Units `hdu` and returns a dictionary
whose keys are the column names and whose values are the column data.  By
default all columns are returned but a selector `sel` may be specified to
select a subset of columns, it must be a callable object such that `sel(key)`
yields `true` if column named `key` (in uppercase letters and trailing spaces
removed) is to be returned.  Units of returned columns, if any, are available
under the name `"\$(key).units"`.

"""
function read_table(hdu::Union{TableHDU,ASCIITableHDU},
                    sel = everything)
    fh = get_file_handle(hdu)
    hdr = read_header(hdu)
    data = Dict{String,Any}()
    ncols = get_integer(hdr, "TFIELDS", 0)
    for k in 1:ncols
        key = fix_name(get_string(hdr, "TTYPE$k", ""))
        if sel(key)
            if haskey(data, key)
                warn("Duplicate column name: \"", key, "\"")
                continue
            end
            data[key] = read_column(fh, k)
            units = strip(get_string(hdr, "TUNIT$k", ""))
            if length(units) > 0
                data[key*".units"] = units
            end
        end
    end
    return data
end

everything(args...) = true

"""
    OIFITS.read_datablock([T = Float64,] hdu; quiet=false)

reads the OI-FITS data in FITS Header Data Units `hdu` and returns an instance
of one of the `OIDataBlock` sub-types or `nothing` if `hdu` does not contain
OI-FITS data.  Optional argument `T` is the floating-point type to use for the
returned object.  If keyword `quiet` is `true`, no warning messages are
printed.

If keyword `quiet` is `true`, no warning messages are printed.

"""
read_datablock(hdu::HDU; kwds...) = read_datablock(Float64, hdu; kwds...),

# OI-FITS data-blocks are stored as FITS binary tables, hence returns nothing
# for any other HDU type.
read_datablock(::Type{<:AbstractFloat}, hdu::HDU; kwds...) = nothing

function read_datablock(::Type{T},
                        hdu::TableHDU;
                        quiet::Bool=false) where {T<:AbstractFloat}
    # Read the header of the binary table and check extension name.
    local extname::String
    hdr = read_header(hdu)
    let val = get_value(hdr, "EXTNAME", nothing)
        if ! isa(val, AbstractString)
            quiet || warn(val === nothing
                          ? "missing keyword \"EXTNAME\""
                          : "EXTNAME value is not a string")
            return nothing
        end
        extname = fix_name(val)
    end
    startswith(extname, "OI_") || return nothing
    if ! (extname in EXTNAMES)
        quiet || warn("unknown OI-FITS extension \"", extname, "\"")
        return nothing
    end

    # Get revision number.
    local revn::Int
    let val = get_value(hdr, "OI_REVN", nothing)
        if ! isa(val, Integer)
            quiet || warn(val === nothing
                          ? "missing keyword \"OI_REVN\""
                          : "\"OI_REVN\" value is not an integer")
            return nothing
        end
        revn = val
    end
    if revn < 1
        quiet || warn("invalid \"OI_REVN\" value: ", revn)
        return nothing
    end

    # Get format definition.
    local defn::Builder.DataBlockDefinition
    let val = Builder.get_definition(extname, revn, nothing)
        if val === nothing
            quiet || warn("unknown OI-FITS extension \"", extname,
                          "\" revision ", revn)
            return nothing
        end
        defn = val
    end

    # So far so good, make a dictionary of the column numbers of the table
    # indexed by the column names.
    cols = Dict{String,Int}()
    for k in 1:get_integer(hdr, "TFIELDS", 0)
        ttype = get_string(hdr, "TTYPE$k")
        cols[fix_name(ttype)] = k
    end

    # Read the columns into a new data-block instance.
    _read_datablock(get_datablock_type(T, extname), defn, hdu, hdr, cols)
end

function _read_datablock(::Type{T},
                         defn::Builder.DataBlockDefinition,
                         hdu::TableHDU,
                         hdr::FITSHeader,
                         cols::Dict{String,Int}) where {T<:OIDataBlock}
    obj = T()
    obj.revn = defn.revn
    nerrs = 0
    fh = get_file_handle(hdu)
    for field in defn.fields
        spec = defn.spec[field]
        name = spec.name
        if spec.iskeyword
            let val = get_value(hdr, name, nothing)
                if val === nothing
                    if spec.isoptional == false
                        warn("missing keyword \"", name,
                             "\" in OI-FITS extension ", defn.extname)
                        nerrs += 1
                    end
                else
                    _setproperty!(obj, spec.symb, val)
                end
            end
        else
            col = get(cols, name, 0)
            if col < 1
                if spec.isoptional == false
                    warn("missing column \"", name,
                         "\" in OI-FITS extension ", defn.extname)
                    nerrs += 1
                end
            else
                let val = read_column(fh, col, spec.multiplier)
                    _setproperty!(obj, spec.symb, val)
                end
            end
        end
    end
    nerrs > 0 && error("bad OI-FITS extension ", defn.extname)
    Builder.check_contents(obj)
    return obj
end

"""
    OIFITS.load([T = Float64,] src) -> master

reads all OI-FITS data-blocks from FITS file `src` and return an instance of
`OIMaster`.  Argument `src` can be a file name of a FITS handle.  Optional
argument `T` is the floating-point type to use for the returned object.  If
keyword `quiet` is `true`, no warning messages are printed.

"""
load(src::Union{AbstractString,FITS}; kwds...) = load(Float64, src; kwds...)

function load(::Type{T}, src::Union{AbstractString,FITS};
              kwds...) where {T<:AbstractFloat}
    load!(OIMaster{T}(), src; kwds...)
end

"""
    OIFITS.load!(dst, src) -> dst

stores all OI-FITS data-blocks from FITS file `src` in `OIMaster` instance
`dst` and return `dst`.  Argument `src` can be a file name of a FITS handle.
If keyword `quiet` is `true`, no warning messages are printed.

"""
load!(dst::OIMaster, src::AbstractString; kwds...) =
    load!(dst, FITS(src, "r"); kwds...)

function load!(dst::OIMaster{T}, src::FITS; quiet::Bool=true) where {T}
    # Read all contents, skipping first HDU.
    for hdu in 2:length(src)
        let db = read_datablock(T, src[hdu], quiet=quiet)
            if db === nothing
                quiet || println("skipping HDU ", hdu, " (no OI-FITS data)")
            else
                quiet || println("reading OI-FITS ", db.extname, " in HDU ", hdu)
                Builder._push!(dst, db)
            end
        end
    end
    Builder._update_links!(dst)
end


text:
@inline function kernel(a, b, c)
    f(x) = 1 - 1.895 * x^2
    ϵ = 0.1
    return (1 - ϵ) * f(b) + ϵ / 2 * (f(a) + f(c))
end

stencil!(us::AbstractMatrix; kwargs...) = stencil!(kernel, us; kwargs...)

function example_prob(; nx::Integer = 100, nt::Integer = 1000)
    us = zeros(nx, nt)
    us[:, begin] .= range(-0.5, 0.5; length = nx) .+ 0.01
    return us
end

function example(; nx::Integer = 100, nt::Integer = 1000, kwargs...)
    us = example_prob(; nx = nx, nt = nt)
    return stencil!(us; smallsize = 2^10, basesize = Inf, timewindow = Inf, kwargs...)
end

as_basesize_int64(x) = max(2, Int64(x))
function as_basesize_int64(x::AbstractFloat)
    if x == Inf
        return typemax(Int64)
    end
    error("expected an integer or `Inf`; got: ", x)
end

function stencil_parameters(; smallsize, basesize, timewindow)
    return (
        smallsize = as_basesize_int64(smallsize),
        basesize = as_basesize_int64(basesize),
        timewindow = as_basesize_int64(timewindow),
    )
end

const DUMMY_VALUE = 12345

function stencil!(f::F, us::AbstractMatrix; kwargs...) where {F}
    @ifdebug fill!(@view(us[begin+1:end-1, begin+1:end]), DUMMY_VALUE)
    smallsize, basesize, timewindow = stencil_parameters(; kwargs...)
    xs = firstindex(us, 1)+1:lastindex(us, 1)-1
    for tf in firstindex(us, 2)+1:timewindow:lastindex(us, 2)
        tl = Int(min(tf + Int128(timewindow) - 1, lastindex(us, 2)))
        trapezoid_par!(f, us, tf:tl, xs, smallsize, basesize)
    end
    return us
end

function trapezoid_par!(
    f::F,
    us::AbstractMatrix,
    ts::UnitRange,
    xs::UnitRange,
    smallsize::Integer,
    basesize::Integer,
) where {F}
    if length(xs) ≤ basesize
        dx0 = first(xs) > firstindex(us, 1) + 1 ? 1 : 0
        dx1 = last(xs) < lastindex(us, 1) - 1 ? -1 : 0
        trapezoid_seq!(f, us, ts, xs, dx0, dx1, smallsize)
    else
        xm = first(xs) + length(xs) ÷ 2 - 1
        @sync begin
            @spawn trapezoid_par!(f, us, ts, xm+1:last(xs), smallsize, basesize)
            trapezoid_par!(f, us, ts, first(xs):xm, smallsize, basesize)
        end
        # Compute the "inverted trapezoid" in between the DAC'ed cases above
        let ts = first(ts)+1:last(ts),  # "t=1" is already computed
            xs = xm:xm+1,             # "t=2" has two uncomputed states
            dx0 = -1,                   # expand to left
            dx1 = 1
            # expand to right
            trapezoid_basecase!(f, us, ts, xs, dx0, dx1)
        end
    end
    return
end

function trapezoid_seq!(
    f::F,
    us::AbstractMatrix,
    ts::UnitRange,
    xs::UnitRange,
    dx0,
    dx1,
    smallsize::Integer,
) where {F}
    @ifdebug @check dx0 ∈ (-1, 0, 1) && dx1 ∈ (-1, 0, 1)
    if max(length(xs), length(ts)) ≤ smallsize
        trapezoid_basecase!(f, us, ts, xs, dx0, dx1)
    else
        doublewidth = 2 * length(xs) + (dx1 - dx0) * length(ts)
        if doublewidth ≥ 4 * length(ts)  #  width ≥ 2 height  ⇒  space cut
            xm = first(xs) + length(xs) ÷ 2 - 1
            xsl = first(xs):xm
            xsr = xm+1:last(xs)
            trapezoid_seq!(f, us, ts, xsl, dx0, -1, smallsize)
            trapezoid_seq!(f, us, ts, xsr, -1, dx1, smallsize)
        else  # width < 2 height  ⇒  time cut
            half = length(ts) ÷ 2
            tm = first(ts) + half - 1
            tsl = first(ts):tm
            tsr = tm+1:last(ts)
            x0 = first(xs) + dx0 * half
            x1 = last(xs) + dx1 * half
            trapezoid_seq!(f, us, tsl, xs, dx0, dx1, smallsize)
            trapezoid_seq!(f, us, tsr, x0:x1, dx0, dx1, smallsize)
        end
    end
    return
end

function trapezoid_basecase!(f, us, ts, xs, dx0, dx1)
    @ifdebug @check dx0 ∈ (-1, 0, 1) && dx1 ∈ (-1, 0, 1)
    x0 = first(xs)
    x1 = last(xs)
    for t in ts
        # @simd ivdep for x in x0:x1
        for x in x0:x1
            @ifdebug assert_us_not_written(us, x, t)
            @inbounds us[x, t] = f(us[x-1, t-1], us[x, t-1], us[x+1, t-1])
        end
        x0 += dx0
        x1 += dx1
    end
end

function assert_us_not_written(us, x, t)
    u = us[x, t]
    if u != DUMMY_VALUE
        error("us[x, t] (x=$x, t=$t) already has value $u")
    end
end


text:
export page, row, cell
import Stipple

"""
    `function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                    channel::String = Genie.config.webchannels_default_route , head_content::String = "", kwargs...)`

Generates the HTML code corresponding to an SPA (a single page application), defining the root element of the Vue app.

### Example

```julia
julia> page(:elemid, [
        span("Hello", @text(:greeting))
        ])
"<!DOCTYPE html>\n<html><head><title></title><meta name=\"viewport\" content=\"width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui\" /></head><body class style><link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><div id=elemid><span v-text='greeting'>Hello</span></div><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/underscore-min.js\"></script><script src=\"/js/vue.js\"></script><script src=\"/js/quasar.umd.min.js\"></script>\n<script src=\"/js/apexcharts.min.js\"></script><script src=\"/js/vue-apexcharts.min.js\"></script><script src=\"/js/stipplecore.js\" defer></script><script src=\"/js/vue_filters.js\" defer></script></body></html>"
```
"""
function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "container", style::String = "",
              channel::String = Genie.config.webchannels_default_route, head_content::String = "",
              prepend::Union{S,Vector} = "", append::Union{T,Vector} = [],
              core_theme::Bool = true, kwargs...)::ParsedHTMLString where {S<:AbstractString,T<:AbstractString}
  layout(
    [
      join(prepend)
      Genie.Renderer.Html.div(id = elemid, args...; class = class, kwargs...)
      join(append)
    ],
    partial = partial, title = title, style = style, head_content = head_content, channel = channel,
    core_theme = core_theme)
end
function page(model::T, args...; kwargs...)::ParsedHTMLString where {T<:Stipple.ReactiveModel}
  hasproperty(model, :channel) || throw(Stipple.MissingPropertyException(:channel, model))

  page(vm(model), args...; channel = model.channel, kwargs...)
end

"""
    `function row(args...; kwargs...)`

Creates a `div` HTML element with a CSS class named `row`. This works with Stipple's Twitter Bootstrap to create the
responsive CSS grid of the web page. The `row` function creates rows which should include `cell`s.

### Example

```julia
julia> row(span("Hello"))
"<div class=\"row\"><span>Hello</span></div>"
```
"""
function row(args...; kwargs...)
  kwargs = NamedTuple(Dict{Symbol,Any}(kwargs...), :class, "row")

  Genie.Renderer.Html.div(args...; kwargs...)
end

"""
    `function cell(args...; size::Int=0, kwargs...)`

Creates a `div` HTML element with CSS classes named `col col-12` and `col-sm-$size`.
This works with Stipple's Twitter Bootstrap to create the responsive CSS grid of the web page. The `cell`s should be
included within `row`s.

### Example

```julia
julia> row(cell(size=2, span("Hello")))
"<div class=\"row\"><div class=\"col col-12 col-sm-2\"><span>Hello</span></div></div>"
```
"""
function cell(args...; size::Int=0, kwargs...)
  kwargs = NamedTuple(Dict{Symbol,Any}(kwargs...), :class, "col col-12 col-sm$(size > 0 ? "-$size" : "")")

  Genie.Renderer.Html.div(args...; kwargs...)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

function apply_collision!(u, dt, collision_op,
                          mesh::AbstractMesh{3}, equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      u_node = get_node_vars(u, equations, dg, i, j, k, element)
      update = collision_op(u_node, dt, equations)
      add_to_node_vars!(u, update, equations, dg, i, j, k, element)
    end
  end

  return nothing
end

end # @muladd


text:
# dummy module for testing assignment in parcel
module E
    struct ET
        x::Int
    end
end


text:
module NumericPredicates

import Base: IEEEFloat,
             iszero, isone,
             isinteger, iseven, isodd,
             isreal, isfinite, isinf, isnan, issubnormal

export negone, two, half,
       isnegone, istwo, ishalf,
       ispos, isneg, isnonneg, isnonpos,
       isfloat, isrational, iscomplex, iscomplexint, iscomplexfloat,
       isfractional, ismixed, isextended, nextint, prevint,

 using IEEEFloats

function nextint(x::T) where {T<:Union{Signed,Unsigned}}
     x === typemax(T) && throw(ErrorException("OverflowError($x)"))
     x += one(T)
     return x
end

function prevint(x::T) where {T<:Union{Signed,Unsigned}}
     x === typemin(T) && throw(ErrorException("UnderflowError($x)"))
     x -= one(T)
     return x
end

function nextint(x::T) where {T<:IEEEFloat}
     x = floor(x)
     x === floatintmax(T) && throw(ErrorException("OverflowError($x)"))
     x += one(T)
     return x
end

function prevint(x::T) where {T<:IEEEFloat}
     x = ceil(x)
     x === floatintmin(T) && throw(ErrorException("UnderflowError($x)"))
     x -= one(T)
     return x
end

function nextint(x::T) where {T<:Real}
     x = floor(x)
     z = x + one(T)
     z === x && throw(ErrorException("OverflowError($x)"))
     return z
end

function prevint(x::T) where {T<:Real}
     x = ceil(x)
     z = x - one(T)
     z === x && throw(ErrorException("UnderflowError($x)"))
     return z
end

isfractional(x::T) where {T<:Union{Signed,Unsigned}} = false
isfractional(x::T) where {T<:Real} = -one(T) < x < one(T)
ismixed(x::T) where {T<:Union{Signed,Unsigned}} = false
ismixed(x::T) where {T<:AbstractFloat} = abs(x) > one(T) && !isinteger(x)

end # module NumericPredicates


text:
# dummycamera

dummycamTimer = nothing
exposurefactor = camSettings.exposureTime/12.0 #dummy setting to simulate ideal 12.0 time

# Typically imported functions, made to act like t dummy camera
function getimage!(cam::String,image::Array{UInt8};normalize::Bool=false)
    image .= round.(UInt8,clamp.(rand(UInt8,size(image,1),size(image,2)).*exposurefactor,UInt8(0),UInt8(255)))
    wait(dummycamTimer)
end
function exposure!(cam::String)
    exposure = 12.0 #dummy value
    updateExposureFactor(exposure,camSettings.gain) #dummy brightness change
    return exposure
end
function exposure!(cam::String,exposure::Real)
    updateExposureFactor(exposure,camSettings.gain) #dummy brightness change
end
function gain!(cam::String)
    gain = 0.0 #dummy value
    updateExposureFactor(camSettings.exposureTime,gain) #dummy brightness change
    return gain
end
function gain!(cam::String,gain::Real)
    updateExposureFactor(camSettings.exposureTime,gain) #dummy brightness change
end
function start!(cam::String)
    global dummycamTimer
    dummycamTimer = Timer(0.0,interval=1/camSettings.acquisitionFramerate)
end
function stop!(cam::String)
    global dummycamTimer
    close(dummycamTimer)
    dummycamTimer = nothing
end
function framerate!(cam::String,framerate::Real)
    global dummycamTimer
    dummycamTimer = Timer(0.0,interval=1/camSettings.acquisitionFramerate)
end
function imagedims!(cam::String,dims)
    global camImage
    camImage = Array{UInt8}(undef,dims...)
end
function offsetdims!(cam::String, dims)
    # do something with dims
end
function sensordims(cam::String)
    return (2048,1536)
end
function isrunning(cam::String)
    return camRunning
end

# Helper functions for dummy behaviour
function updateExposureFactor(exposure,gain)
    global exposurefactor
    exposurefactor = (exposure/12.0)*(1.0 + gain)
end

# Main functions
function runCamera()
    global dummycamTimer
    global perfGrabFramerate
    global cam, camImage, camImageFrameBuffer
    global camSettingsLimits

    # Initialize dummy camera
    cam = "dummycam"

    camSettingsLimitsUpdater!(cam,camSettingsLimits)

    start!(cam)

    camImage = Array{UInt8}(undef,camSettings.width,camSettings.height)

    perfGrabTime = time()
    while gui_open
        if isrunning(cam)
            #cim_id, cim_timestamp, cim_exposure = getimage!(camera, previewimage)
            try
                getimage!(cam,camImage,normalize=false)
            catch err
                if err isa EOFError
                    # Do nothing. This happens if the camera is stopped before the camRunning
                    # bool is set due to async
                else
                    rethrow()
                end
            end
            # Loop timing
            perfGrabFramerate = 1/(time() - perfGrabTime)
            perfGrabTime = time()
        end
        yield()
    end
    stop!(cam)
end


text:
@testset "43.multiply-strings.jl" begin
    @test trap([0,1,0,2,1,0,1,3,2,1,2,1]) == 6
    @test trap([4,2,0,3,2,5]) == 9
end


text:
export L1ball

struct L1ball{T <: Real} <: ConvexSet
    center::Vector{T}
    radius::T
end

"""

    L1ball(c, r)

A ball with center `c` and radius `r`.
"""

function project!(s::L1ball{T}, v::Vector{T}, y::Vector{T}) where {T <: Real}
    n = length(y)
    z = sum(abs, y-s.center)
   if z <= s.radius
      v = y
   else
      w = project(Projections.Simplex(s.radius),y)
      sgn = zeros(n)
      d = y-s.center
      for i in 1:n
            if d[i] > 0
                sgn[i] = 1
            elseif d[i] == 0
                sgn[i] = 0
            else
                sgn[i] = -1
            end
      end
        v = s.center .+ sgn .* w
   end
end

function project(s::L1ball{<:Real}, y::Vector{<:Real})
    project!(s, similar(y), y)
end


text:
function extrema(point_array::Vector{<: AbstractPoint}, is_less_x, is_less_y, large = [1e100 0 ; 0 1e100])
    # is_less_x(p1, p2) corresponds to the first dimension output, X

    L = length(point_array)
    if L == 1
        return (point_array[1], point_array[1], point_array[1], point_array[1])
    end

    PointType = eltype(point_array)

    # min_x = PointType(+Inf, 0)
    # max_x = PointType(-Inf, 0)
    # min_y = PointType(0, +Inf)
    # max_y = PointType(0, -Inf)

    # large = 1e300
    min_x = PointType(+large[1,1], +large[1,2])
    max_x = PointType(-large[1,1], -large[1,2])
    min_y = PointType(+large[2,1], +large[2,1])
    max_y = PointType(large[2,2], -large[2,2])

    is_odd = L % 2 == 1
    Lc = is_odd ? L-1 : L

    for kk = 1:2:Lc

        p1 = point_array[kk]
        p2 = point_array[kk+1]

        if is_less_x(p1, p2) # p1.x < p2.x
            if is_less_x(p1, min_x) # min_x.x > p1.x
                min_x = p1
            end
            if is_less_x(max_x, p2) # max_x.x < p2.x
                max_x = p2
            end
        else
            if is_less_x(p2, min_x) # min_x.x > p2.x
                min_x = p2
            end
            if is_less_x(max_x, p1) # max_x.x < p1.x
                max_x = p1
            end
        end

        if is_less_y(p1, p2) # p1.y < p2.y
            if is_less_y(p1, min_y) # min_y.y > p1.y
                min_y = p1
            end
            if is_less_y(max_y, p2) # max_y.y < p2.y
                max_y = p2
            end
        else
            if is_less_y(p2, min_y) # min_y.y > p2.y
                min_y = p2
            end
            if is_less_y(max_y, p1) # max_y.y < p1.y
                max_y = p1
            end
        end
    end
    if is_odd
        p1 = point_array[end]
        # check for the last element too
        if is_less_x(p1, min_x) # min_x.x > p1.x
            min_x = p1
        end
        if  is_less_x(max_x, p1) # max_x.x < p1.x
            max_x = p1
        end
        if is_less_y(p1, min_y) # min_y.y > p1.y
            min_y = p1
        end
        if is_less_y(max_y, p1) # max_y.y < p1.y
            max_y = p1
        end
    end
    return (min_x, max_x, min_y, max_y)
end

function extrema(point_array::Vector{<: AbstractPoint})
    L = length(point_array)
    if L == 1
        return (point_array[1], point_array[1], point_array[1], point_array[1])
    end

    PointType = eltype(point_array)
    min_x = PointType(+Inf, 0)
    max_x = PointType(-Inf, 0)
    min_y = PointType(0, +Inf)
    max_y = PointType(0, -Inf)

    is_odd = L % 2 == 1
    Lc = is_odd ? L-1 : L

    for kk = 1:2:Lc

        p1 = point_array[kk]
        p2 = point_array[kk+1]

        if p1.x < p2.x
            if min_x.x > p1.x
                min_x = p1
            end
            if max_x.x < p2.x
                max_x = p2
            end
        else
            if min_x.x > p2.x
                min_x = p2
            end
            if max_x.x < p1.x
                max_x = p1
            end
        end

        if p1.y < p2.y
            if min_y.y > p1.y
                min_y = p1
            end
            if max_y.y < p2.y
                max_y = p2
            end
        else
            if min_y.y > p2.y
                min_y = p2
            end
            if max_y.y < p1.y
                max_y = p1
            end
        end
    end
    if is_odd
        p1 = point_array[end]
        # check for the last element too
        if min_x.x > p1.x
            min_x = p1
        end
        if max_x.x < p1.x
            max_x = p1
        end
        if min_y.y > p1.y
            min_y = p1
        end
        if max_y.y < p1.y
            max_y = p1
        end
    end
    return (min_x, max_x, min_y, max_y)
end

function skew_extrema(point_array::Vector{<: AbstractPoint})
    geom.extrema( point_array,
                  (p1, p2) -> (p1.x + p1.y) < (p2.x + p2.y),
                  (p1, p2) -> (p1.y - p2.y) < (p1.x - p2.x) )
end

function extrema_y(point_array::Vector{<: AbstractPoint})
    L = length(point_array)
    if L == 1
        return (point_array[1], point_array[1], 1, 1)
    end

    PointType = eltype(point_array)
    min_y = PointType(0, +Inf)
    max_y = PointType(0, -Inf)
    idx_min = 0
    idx_max = 0

    is_odd = L % 2 == 1
    Lc = is_odd ? L-1 : L

    for kk = 1:2:Lc

        p1 = point_array[kk]
        p2 = point_array[kk+1]

        if p1.y < p2.y
            if min_y.y > p1.y
                min_y = p1
                idx_min = kk
            end
            if max_y.y < p2.y
                max_y = p2
                idx_max = kk+1
            end
        else
            if min_y.y > p2.y
                min_y = p2
                idx_min = kk+1
            end
            if max_y.y < p1.y
                max_y = p1
                idx_max = kk
            end
        end
    end
    if is_odd
        p1 = point_array[end]
        # check for the last element too
        if min_y.y > p1.y
            min_y = p1
            idx_min = L
        end
        if max_y.y < p1.y
            max_y = p1
            idx_max = L
        end
    end
    return (min_y, max_y, idx_min, idx_max)
end

function extrema_x(point_array::Vector{<: AbstractPoint})
    L = length(point_array)
    if L == 1
        return (point_array[1], point_array[1], 1, 1)
    end

    PointType = eltype(point_array)
    min_x = PointType(+Inf, 0)
    max_x = PointType(-Inf, 0)

    is_odd = L % 2 == 1
    Lc = is_odd ? L-1 : L

    for kk = 1:2:Lc

        p1 = point_array[kk]
        p2 = point_array[kk+1]

        if p1.x < p2.x
            if min_x.x > p1.x
                min_x = p1
                idx_min = kk
            end
            if max_x.x < p2.x
                max_x = p2
                idx_max = kk
            end
        else
            if min_x.x > p2.x
                min_x = p2
                idx_min = kk
            end
            if max_x.x < p1.x
                max_x = p1
                idx_max = kk
            end
        end
    end
    if is_odd
        p1 = point_array[end]
        # check for the last element too
        if min_x.x > p1.x
            min_x = p1
            idx_min = L
        end
        if max_x.x < p1.x
            max_x = p1
            idx_max = L
        end
    end
    return (min_x, max_x, idx_min, idx_max)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using MdfReader
const m = MdfReader
path = "test\\testfiles\\test.mdf"
# mdf = MDF(path)

##

##
# abstract type bitstypes end
# struct bittype <: bitstypes end
# struct nonbittype <: bitstypes end
# # nameing is stolen from StructIO.jl
# bitstypes(a::DataType) = isbitstype(a) ? bittype : nonbittype

##
a = m.BlockHeader("HD",100)

##

# function Base.read(io,a::Union{Type{BLOCK},Type{MdfReader.BlockHeader}},

a = HDBLOCK()
pointer(a::BLOCK) = convert(Ptr{typeof(a)},pointer_from_objref(Ref{a}))
p = pointer(a)
padding = Base.padding(a)