text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.
# This file includes parts from Combinatorics.jl in modified form.

# Version of: UTC 2020-01-30 22:26:36
# 323c3f30-43a7-11ea-3bfd-cb8fb16dde07

# Do not edit this file, it is generated from the modules and will be overwritten!
# Edit the modules in the src directory and build this file with BuildSequences.jl!

tstdir = realpath(joinpath(dirname(@__FILE__)))
srcdir = joinpath(dirname(tstdir), "src")
tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)
srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)
module perftests
using IntegerSequences, DataStructures, Dates, InteractiveUtils
InteractiveUtils.versioninfo()
start = Dates.now()
# +++ Abundant.jl +++
GC.gc()
println("\nTEST: L005101(5000)")
@time L005101(5000)
# +++ AltPartitions.jl +++
# +++ AndreNumbers.jl +++
GC.gc()
println("\nTEST: (for m ∈ 1:20, n ∈ 0:100")
@time (for m ∈ 1:20, n ∈ 0:100
André(m, n)
end)
# +++ BellNumbers.jl +++
GC.gc()
BellNumberList(5)
println("\nTEST: BellNumberList(1000)")
@time BellNumberList(1000)
println("\nTEST: BellTriangle(100, [1 for _ ∈ 1:100])")
@time BellTriangle(100, [1 for _ ∈ 1:100])
# +++ BernoulliNumbers.jl +++
GC.gc()
println("\nTEST: (for m ∈ 1:10, n ∈ 0:499 BernoulliInt(m, n) end)")
@time (for m ∈ 1:10, n ∈ 0:499 BernoulliInt(m, n) end)
println("\nTEST: (for m ∈ 1:10 BernoulliIntList(m, 500) end)")
@time (for m ∈ 1:10 BernoulliIntList(m, 500) end)
GC.gc()
println("\nTEST: (for n ∈ 0:1000 Bernoulli(n) end)")
@time (for n ∈ 0:1000 Bernoulli(n) end)
println("\nTEST: BernoulliList(1000)")
@time BernoulliList(1000)
println("\nTEST: (for n ∈ 0:1000 V027641(n) end)")
@time (for n ∈ 0:1000 V027641(n) end)
GC.gc()
println("\nTEST: (for n ∈ 0:10000 V195441(n) end)")
@time (for n ∈ 0:10000 V195441(n) end)
# +++ BinaryInteger.jl +++
println("\nTEST: [V295513(k) for k ∈ 0:100000]")
@time [V295513(k) for k ∈ 0:100000]
# +++ BinaryQF.jl +++
GC.gc()
println("\nTEST: L034017(10000)")
@time L034017(10000)
# +++ BooleanOperations.jl +++
# +++ CantorMachines.jl +++
# +++ CarmichaelNumbers.jl +++
println("\nTEST: L = L002997(30)")
@time L = L002997(30)
# +++ ClausenNumbers.jl +++
GC.gc()
println("\nTEST: ClausenNumberList(10000)")
@time ClausenNumberList(10000)
println("\nTEST: (for n ∈ 0:10000")
@time (for n ∈ 0:10000
ClausenNumber(n)
end)
# +++ CombinationsIterator.jl +++
# +++ Compositions.jl +++
GC.gc()
println("\nTEST: L097805(500)")
@time L097805(500)
# +++ Counts.jl +++
println("\nTEST: [V257993(n) for n ∈ 1:10000]")
@time [V257993(n) for n ∈ 1:10000]
GC.gc()
println("\nTEST: PrimePiList(10000)")
@time PrimePiList(10000)
# +++ CyclotomicBinaryForms.jl +++
println("\nTEST: L299214(100)")
@time L299214(100)
# +++ DedekindEta.jl +++
println("\nTEST: PartitionNumberList(10000)")
@time PartitionNumberList(10000)
println("\nTEST: L000731(10000)")
@time L000731(10000)
println("\nTEST: RamanujanTauList(10000)")
@time RamanujanTauList(10000)
# +++ DelehamDelta.jl +++
GC.gc()
println("\nTEST: T225478(100)")
@time T225478(100)
println("\nTEST: T055883(100)")
@time T055883(100)
# +++ EulerTransforms.jl +++
println("\nTEST: (for n ∈ 0:150 V107895(n) end)")
@time (for n ∈ 0:150 V107895(n) end)
# +++ Fibonacci.jl +++
GC.gc()
println("\nTEST: (for n ∈ 1:1000")
@time (for n ∈ 1:1000
V000045(n)
end)
println("\nTEST: (for fib ∈ I000045(10000)")
@time (for fib ∈ I000045(10000)
end)
println("\nTEST: L000045(10000)")
@time L000045(10000)
# +++ FigurativeNumbers.jl +++
println("\nTEST: (for n ∈ 1:1000")
@time (for n ∈ 1:1000
V000326(n)
end)
# +++ GaussFactorials.jl +++
# +++ GeneralizedBinomial.jl +++
println("\nTEST: (for n ∈ 0:10000")
@time (for n ∈ 0:10000
Binomial(2n, n)
end)
println("\nTEST: (for n ∈ -100:100, k ∈ -100:100")
@time (for n ∈ -100:100, k ∈ -100:100
Binomial(n, k)
end)
println("\nTEST: (for k ∈ -10000:10000")
@time (for k ∈ -10000:10000
Binomial(-5, k)
end)
# +++ GeneralizedFibonacci.jl +++
println("\nTEST: (for n ∈ 0:1000 Multinacci(3, n) end)")
@time (for n ∈ 0:1000 Multinacci(3, n) end)
# +++ HighlyAbundant.jl +++
GC.gc()
println("\nTEST: (for n ∈ 1:100 V034885(n) end)")
@time (for n ∈ 1:100 V034885(n) end)
# +++ Hyper1F1.jl +++
GC.gc()
println("\nTEST: (for n ∈ 0:150")
@time (for n ∈ 0:150
V000262(n)
end)
# +++ IntPartitions.jl +++
println("\n-- performance test")
println("\nTEST: for n ∈ 1:50, p ∈ IntegerPartitions(n, byNumPart) p end")
@time for n ∈ 1:50, p ∈ IntegerPartitions(n, byNumPart) p end
println("\nTEST: for n ∈ 1:50, p ∈ IntegerPartitions(n, byMaxPart) p end")
@time for n ∈ 1:50, p ∈ IntegerPartitions(n, byMaxPart) p end
# +++ JacobiTheta.jl +++
println("\nTEST: L000143(100000)")
@time L000143(100000)
println("\nTEST: L035016(100000)")
@time L035016(100000)
# +++ Kolakoski.jl +++
println("\nTEST: I000002(100000)")
@time I000002(100000)
println("\nTEST: KolakoskiList(100000)")
@time KolakoskiList(100000)
# +++ LahNumbers.jl +++
println("\nTEST: Lah(1000)")
@time Lah(1000)
println("\nTEST: RowSum(LahTriangle(1000))")
@time RowSum(LahTriangle(1000))
# +++ LandauConstants.jl +++
"V001511" |> println
println("\nTEST: for n ∈ 0:100 V001511(n) end")
@time for n ∈ 0:100 V001511(n) end
"V163590" |> println
println("\nTEST: for n ∈ 0:100 V163590(n) end")
@time for n ∈ 0:100 V163590(n) end
"V327492" |> println
println("\nTEST: for n ∈ 0:100 V327492(n) end")
@time for n ∈ 0:100 V327492(n) end
"V005187" |> println
println("\nTEST: for n ∈ 0:100 V005187(n) end")
@time for n ∈ 0:100 V005187(n) end
"V001790" |> println
println("\nTEST: for n ∈ 0:100 V001790(n) end")
@time for n ∈ 0:100 V001790(n) end
"V277233" |> println
println("\nTEST: for n ∈ 0:100 V277233(n) end")
@time for n ∈ 0:100 V277233(n) end
"V327495" |> println
println("\nTEST: for n ∈ 0:100 V327495(n) end")
@time for n ∈ 0:100 V327495(n) end
# +++ NarayanaCows.jl +++
println("\nTEST: (for cow ∈ NarayanasCows(10000)")
@time (for cow ∈ NarayanasCows(10000)
end)
# +++ NodesAndEdges.jl +++
println("\nTEST: [V000055(n) for n in 0:99]")
@time [V000055(n) for n in 0:99]
# +++ NumberTheory.jl +++
println("\nTEST: [Divisors(n) for n ∈ 1:10000]")
@time [Divisors(n) for n ∈ 1:10000]
println("\nTEST: [Radical(n)  for n ∈ 1:10000]")
@time [Radical(n)  for n ∈ 1:10000]
# +++ OEISUtils.jl +++
# +++ OrthoPolynomials.jl +++
GC.gc()
println("\nTEST: T111062(500)")
@time T111062(500)
GC.gc()
println("\nTEST: T066325(500)")
@time T066325(500)
GC.gc()
println("\nTEST: T053120(500)")
@time T053120(500)
# +++ Polynomials.jl +++
# +++ PrimesIterator.jl +++
println("\nTEST: primepi = PrimePi(1000000)")
@time primepi = PrimePi(1000000)
# +++ PrimeSwingFactorial.jl +++
GC.gc()
println("\nTEST: Base.factorial(BigInt(1000000))")
@time Base.factorial(BigInt(1000000))
println("\nTEST: PSfactorial(1000000)")
@time PSfactorial(1000000)
# +++ Products.jl +++
GC.gc()
println("\nTEST: (for n ∈ 1:10000")
@time (for n ∈ 1:10000
F!(n)
end)
GC.gc()
println("\nTEST: (for n ∈ 1:1000")
@time (for n ∈ 1:1000
V000407(n)
end)
GC.gc()
println("\nTEST: (for n ∈ 1:200, k ∈ 1:200")
@time (for n ∈ 1:200, k ∈ 1:200
RisingFactorial(n, k)
end)
# +++ QueensProblems.jl +++
println("\nTEST: L319284(15)")
@time L319284(15)
# +++ RationalTrees.jl +++
GC.gc()
println("\nTEST: EuclidTree(100)")
@time EuclidTree(100)
GC.gc()
println("\nTEST: [SchinzelSierpinskiEncoding(l) for l ∈ EuclidTree(100)]")
@time [SchinzelSierpinskiEncoding(l) for l ∈ EuclidTree(100)]
# +++ RecordSearch.jl +++
# +++ RiordanSquares.jl +++
# +++ SelfConvolutive.jl +++
GC.gc()
println("\nTEST: L005411(500)")
@time L005411(500)
# +++ SeqTests.jl +++
# +++ SeqUtils.jl +++
# +++ SeriesExpansion.jl +++
println("\nTEST: L068875(1000)")
@time L068875(1000)
# +++ SetPartitions.jl +++
# +++ SpigotPi.jl +++
println("\nTEST: Pi(1600)")
@time Pi(1600)
# +++ StirlingLahNumbers.jl +++
GC.gc()
println("\nTEST: (for row ∈ T132393(500)")
@time (for row ∈ T132393(500)
end)
GC.gc()
println("\nTEST: (for row ∈ T048993(500)")
@time (for row ∈ T048993(500)
end)
GC.gc()
println("\nTEST: (for row ∈ T271703(500)")
@time (for row ∈ T271703(500)
end)
# +++ SwingFactorial.jl +++
GC.gc()
n = 1000000
println("\nTEST: Base.factorial(BigInt(n))")
@time Base.factorial(BigInt(n))
println("\nTEST: Sfactorial(n)")
@time Sfactorial(n)
println("n = $n")
# +++ Triangles.jl +++
# +++ UlamNumbers.jl +++
GC.gc()
for n ∈ 6:10
println("\nTEST: UlamList(2^n)")
@time UlamList(2^n)
end
# +++ ZumkellerNumbers.jl +++
GC.gc()
println("\nTEST: (for n ∈ 1:2000")
@time (for n ∈ 1:2000
isZumkeller(n)
end)
println("\nTEST: L083207(500)")
@time L083207(500)
stop = Dates.now()
tdiff = stop - start
println("\nJulia version: " * string(VERSION) )
println(start)
println("Total test time: ", tdiff)
end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
#!/usr/bin/env julia

function reversort!(case, N, L)
    c = @inbounds sum(1:N - 1) do i
        @views j = i + argmin(L[i:end]) - 1
        @views L[i:j] .= Iterators.reverse(L[i:j])
        j - i + 1
    end
    println("Case #$case: $c")
end

#function cost(case, N, L)
#    c = sum(((i, x),) -> abs(x - i + 1), enumerate(L))
#    println("Case #$case: $c")
#end

function main()
    T = parse(Int, readline())
    for case in 1:T
        N = parse(Int, readline())
        L = map(x -> parse(Int, x), split(readline()))
        reversort!(case, N, L)
    end
end

isinteractive() || main()

ext: jl

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
using Atom
using Base.Test

include("./eval.jl")
include("./utils.jl")


text:
for (op,cfunc) in ((:factorial,:arb_fac_ui), (:doublefactorial,:arb_doublefac_ui))
  @eval begin
    function ($op){P}(x::ArbFloat{P})
      signbit(x) && throw(ErrorException("Domain Error: argument is negative"))
      y = trunc(UInt, x)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, UInt, Int), &z, y, P)
      z
    end
  end
end

function doublefactorial{R<:Real}(xx::R)
   P = precision(ArbFloat)
   x = convert(ArbFloat{P},xx)
   doublefactorial(x)
end

for (op,cfunc) in ((:risingfactorial,:arb_rising),)
  @eval begin
    function ($op){P}(x::ArbFloat{P}, y::ArbFloat{P})
      signbit(x) && throw(ErrorException("Domain Error: argument is negative"))
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    function ($op){R<:Real,P}(xx::R, y::ArbFloat{P}, prec::Int=P)
      x = convert(ArbFloat{P},xx)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    function ($op){R<:Real,P}(x::ArbFloat{P}, yy::R, prec::Int=P)
      y = convert(ArbFloat{P},yy)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    function ($op){R1<:Real,R2<:Real}(xx::R1, yy::R2)
      P = precision(ArbFloat)
      x = convert(ArbFloat{P},xx)
      y = convert(ArbFloat{P},yy)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
  end
end

for (op,cfunc) in ((:agm, :arb_agm), (:polylog, :arb_polylog))
  @eval begin
    function ($op){P}(x::ArbFloat{P}, y::ArbFloat{P}, prec::Int=P)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    function ($op){R<:Real,P}(xx::R, y::ArbFloat{P}, prec::Int=P)
      x = convert(ArbFloat{P},xx)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    function ($op){R<:Real,P}(x::ArbFloat{P}, yy::R, prec::Int=P)
      y = convert(ArbFloat{P},yy)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
    function ($op){R1<:Real,R2<:Real}(xx::R1, yy::R2)
      P = precision(ArbFloat)
      x = convert(ArbFloat{P},xx)
      y = convert(ArbFloat{P},yy)
      z = ArbFloat{P}()
      ccall(@libarb($cfunc), Void, (Ptr{ArbFloat}, Ptr{ArbFloat}, Ptr{ArbFloat}, Int), &z, &x, &y, P)
      z
    end
  end
end


text:
using Plots
import JSON

function plot_anderson(folder; tol=1e-10, iteration=40, ymax=1e2, letter=nothing)
    mixing = "None"
    isfile(joinpath(folder, mixing * ".optimal_damping.json")) || (mixing = "Kerker")

    energies = Dict(m => Float64.(open(JSON.parse,
                                       joinpath(folder, "anderson_$m.json"))["energies"])
                    for m in (0, 5, 10))
    mintotal = minimum(minimum, values(energies))
    mintotal -= tol / 100  # To avoid zero error (and thus -Inf on a log scale)

    errors = Dict{Int, Vector{Float64}}()
    for m in sort(collect(keys(energies)))
        values = abs.(energies[m] .- mintotal)
        idx = findfirst(values .< 1e-16)
        isnothing(idx) || (values = values[1:idx-1])
        isempty(values) && continue
        errors[m] = values
    end

    maxiter = 0
    ms = [:x, :+, :utriangle, :dtriangle]
    p = plot(legend=:topright, grid=false)
    for (i, m) in enumerate(sort(collect(keys(errors))))
        maxiter = max(maxiter, length(errors[m]))
        plot!(p, errors[m]; label="Anderson m=$m", m=ms[i], markerstrokecolor=i,
              linewidth=2)
    end

    fn = joinpath(folder, "$mixing.optimal_damping.json")
    κ = open(JSON.parse, fn)["condition_number"]
    rate_simple = (κ - 1) / (κ + 1)                # Convergence rate for simple mixing
    rate_anderson = (sqrt(κ) - 1) / (sqrt(κ) + 1)  # Convergence rate for anderson

    # Square because energy converges quadratically (as opposed to density residual)
    x = collect(1:0.1:maxiter)
    plot!(p, x, errors[0][1]  * rate_simple.^2x, label="Simple rate", ls=:dash,
          color=:black, linewidth=2)
    plot!(p, x, errors[10][1] * rate_anderson.^2x, label="Anderson rate", ls=:dot,
          color=:black, linewidth=2)

    xlims!(p, (0, iteration))
    ylims!(p, (tol, ymax))
    xaxis!(p, "Iteration")
    yaxis!(p, "Total energy absolute error", :log10)
    if !isnothing(letter)
        plot!(p, [1.1, 8], [18, 18], label="", line=(0.5, :white), lw=15)
        annotate!(p, 1.2, 1.8e1, Plots.text(letter, 10, :left))
    end

    p
end


text:
using NCDatasets
using Formatting
using ArgParse
using Statistics
using JSON

include("constants.jl")
include("CESMReader.jl")

using .CESMReader

function mreplace(x)
    return replace(x, missing=>NaN)
end

function integrate(x, dydx)
    y = copy(dydx) * 0.0

    for i = 2:length(x)
        y[i] = y[i-1] + (dydx[i-1] + dydx[i]) * (x[i] - x[i-1]) / 2.0
    end

    return y
end

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file-prefix"
            help = "Data filename prefix including folder and path until the timestamp. File extension `nc` is assumed."
            arg_type = String
            required = true

        "--data-file-timestamp-form"
            help = "Data filename timestamp form. Either `YEAR` or `YEAR_MONTH`."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--beg-year"
            help = "Year of begin."
            arg_type = Int64
            required = true

        "--end-year"
            help = "Year of end."
            arg_type = Int64
            required = true

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["domain-file"], "r") do ds
    global mask = ds["mask"][:] |> mreplace
    global lat = ds["yc"][1, :] |> mreplace

    global lat_weight = Re * cos.(deg2rad.(lat)) * 2π
    global y          = Re * deg2rad.(lat)

    mask[mask.!=0] .= 1
end

let

    if parsed["data-file-timestamp-form"] == "YEAR"
        filename_format = format("{:s}{{:04d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR
    elseif parsed["data-file-timestamp-form"] == "YEAR_MONTH"
        filename_format = format("{:s}{{:04d}}-{{:02d}}.nc", joinpath(parsed["data-file-prefix"]))
        form = :YEAR_MONTH
    end

    fh = FileHandler(filename_format=filename_format, form=form)

    beg_t = (parsed["beg-year"] - 1) * 12 + 1
    end_t = (parsed["end-year"] - 1) * 12 + 12

    FSNT, FLNT, FSNS, FLNS, SHFLX, LHFLX = getData(fh, ["FSNT", "FLNT", "FSNS", "FLNS", "SHFLX", "LHFLX"], (parsed["beg-year"], parsed["end-year"]), (:, :))

    global FFLX_TOA = ( - mean(FSNT, dims=(1, ))  + mean(FLNT,  dims=(1, )) )[1, :, :]
    global FFLX_SFC = ( - mean(FSNS, dims=(1, ))  + mean(FLNS,  dims=(1, )) )[1, :, :]
    global HFLX_SFC = (   mean(SHFLX, dims=(1, )) + mean(LHFLX, dims=(1, )) )[1, :, :]

    global (Ny, Nt) = size(FFLX_TOA)

    if mod(Nt, 12) != 0
        throw(ErrorException("Time should be a multiple of 12."))
    end

    global nyears = Int64(Nt / 12)
    println(format("We got {:d} years of data.", nyears))
end

IAET       = zeros(Float64, Ny, Nt)
EFLX_CONV  = zeros(Float64, Ny, Nt)

IAET_AM    = zeros(Float64, Ny, nyears)

for t = 1:Nt

    FFLX_TOA[:, t] .*= lat_weight
    FFLX_SFC[:, t] .*= lat_weight
    HFLX_SFC[:, t] .*= lat_weight
    EFLX_CONV[:, t] = - ( FFLX_TOA[:, t] - FFLX_SFC[:, t] - HFLX_SFC[:, t] )
    IAET[:, t] = integrate(y, EFLX_CONV[:, t])

end

for j = 1:Ny, t = 1:nyears
    IAET_AM[j, t] = mean(IAET[j, (t-1)*12+1:t*12])
end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "Ny", Ny)
    defDim(ds, "time", Inf)

    for (varname, vardata, vardim, attrib) in [
        ("FFLX_TOA",  FFLX_TOA,  ("Ny", "time"), Dict()),
        ("FFLX_SFC",  FFLX_SFC,  ("Ny", "time"), Dict()),
        ("HFLX_SFC",  HFLX_SFC,  ("Ny", "time"), Dict()),
        ("EFLX_CONV", EFLX_CONV, ("Ny", "time"), Dict()),
        ("IAET",      IAET,      ("Ny", "time"), Dict()),
        ("FFLX_TOA_mean",  mean(FFLX_TOA, dims=(2,))[:, 1],  ("Ny", ), Dict()),
        ("FFLX_SFC_mean",  mean(FFLX_SFC, dims=(2,))[:, 1],  ("Ny", ), Dict()),
        ("HFLX_SFC_mean",  mean(HFLX_SFC, dims=(2,))[:, 1],  ("Ny", ), Dict()),
        ("EFLX_CONV_mean", mean(EFLX_CONV, dims=(2,))[:, 1], ("Ny", ), Dict()),
        ("IAET_AM",        mean(IAET_AM, dims=(2,))[:, 1],   ("Ny", ), Dict()),
        ("IAET_AMSTD",     std( IAET_AM, dims=(2,))[:, 1],   ("Ny", ), Dict()),
    ]

        println("Doing var: ", varname)

        var = defVar(ds, varname, Float64, vardim)
        var.attrib["_FillValue"] = 1e20

        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
multiply.jl defines methods for the *(x, y) function
"""

# --- Exports

import Base.:(*)

# --- Method definitions

# ------ Specializations involving an AbstractMultivector instance

# M::AbstractMultivector, B::AbstractBlade
# B::AbstractBlade, M::AbstractMultivector
*(M::AbstractMultivector, B::AbstractBlade) =
    Multivector(map(C -> C * B, blades(M)))

*(B::AbstractBlade, M::AbstractMultivector) =
    Multivector(map(C -> B * C, blades(M)))

# M::AbstractMultivector, B::AbstractScalar,
# B::AbstractScalar, M::AbstractMultivector
*(M::AbstractMultivector, B::AbstractScalar) = B * M
*(B::AbstractScalar, M::AbstractMultivector) = value(B) * M

# M::AbstractMultivector, B::One
# B::One, M::AbstractMultivector
*(M::AbstractMultivector, B::One) = M
*(B::One, M::AbstractMultivector) = M

# M::AbstractMultivector, B::Zero
# B::Zero, M::AbstractMultivector
*(M::AbstractMultivector, B::Zero) = B
*(B::Zero, M::AbstractMultivector) = B

# M::AbstractMultivector, x::Real
# x::Real, M::AbstractMultivector
*(M::AbstractMultivector, x::Real) = x * M
*(x::Real, M::AbstractMultivector) = Multivector(map(B -> x * B, blades(M)))

# M::AbstractMultivector, v::Vector
# v::Vector, M::AbstractMultivector
*(M::AbstractMultivector, v::Vector) = Multivector(map(B -> B * v, blades(M)))
*(v::Vector, M::AbstractMultivector) = Multivector(map(B -> v * B, blades(M)))

# ------ Specializations involving an AbstractBlade instance

# B::AbstractBlade, C::AbstractScalar
# B::AbstractScalar, B::AbstractBlade
*(B::AbstractBlade, C::AbstractScalar) =
    Blade(B, volume=volume(B) * value(C), copy_basis=false)
*(C::AbstractScalar, B::AbstractBlade) = B * C

# B::AbstractBlade, C::One
# B::One, C::AbstractBlade
*(B::AbstractBlade, C::One) = B
*(B::One, C::AbstractBlade) = C

# B::AbstractBlade, C::Zero
# B::Zero, C::AbstractBlade
*(B::AbstractBlade, C::Zero) = C
*(B::Zero, C::AbstractBlade) = B

# B::AbstractBlade, x::Real
# x::Real, B::AbstractBlade
*(B::AbstractBlade, x::Real) = Blade(B, volume=volume(B) * x)
*(x::Real, B::AbstractBlade) = B * x

# ------ Specializations involving a Blade instance

# B::Blade, C::Blade
function *(B::Blade, C::Blade)
    # Check arguments
    assert_dim_equal(B, C)

    # Compute geometric product
    M = C
    for i in grade(B):-1:1
        M = basis(B)[:, i] * M
    end
    M = volume(B) * M
end

# B::Blade, C::Pseudoscalar
# B::Pseudoscalar, C::Blade
*(B::Blade, C::Pseudoscalar) = B ⋅ C
*(B::Pseudoscalar, C::Blade) = zero(B)

# B::Blade, v::Vector
# v::Vector, B::Blade
function *(B::Blade, v::Vector{<:Real})
    # --- Check arguments

    # Check that B and v are from the same real vector space
    assert_dim_equal(v, B)

    # Check that grade(B) == 1. B * v is not defined when grade(B) > 1
    if grade(B) > 1
        error("Geometric product B * v is not defined when grade(B) > 1")
    end

    # Compute geometric product
    B_dot_v = B ⋅ v
    B_wedge_v = B ∧ v

    if B_dot_v == zero(B_dot_v)
        return B_wedge_v
    elseif B_wedge_v == zero(B_wedge_v)
        return B_dot_v
    end

    Multivector([B_dot_v, B_wedge_v])
end

function *(v::Vector{<:Real}, B::Blade)
    # Check arguments
    assert_dim_equal(v, B)

    # Compute geometric product
    v_dot_B = v ⋅ B
    v_wedge_B = v ∧ B

    if v_dot_B == zero(v_dot_B)
        return v_wedge_B
    elseif v_wedge_B == zero(v_wedge_B)
        return v_dot_B
    end

    Multivector([v_dot_B, v_wedge_B])
end

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
*(B::Pseudoscalar, C::Pseudoscalar) = contract_left(B, C)

# B::Pseudoscalar, C::AbstractScalar
# B::AbstractScalar, C::Pseudoscalar
*(B::Pseudoscalar, C::AbstractScalar) =
    Pseudoscalar(B, value=value(B) * value(C))

*(B::AbstractScalar, C::Pseudoscalar) = C * B

# B::Pseudoscalar, C::One
# B::One, C::Pseudoscalar
*(B::Pseudoscalar, C::One) = B
*(B::One, C::Pseudoscalar) = C

# B::Pseudoscalar, C::Zero
# B::Zero, C::Pseudoscalar
*(B::Pseudoscalar, C::Zero) = C
*(B::Zero, C::Pseudoscalar) = B

# B::Pseudoscalar, x::Real
# x::Real, B::Pseudoscalar
*(B::Pseudoscalar, x::Real) = Pseudoscalar(B, value=x * value(B))
*(x::Real, C::Pseudoscalar) = C * x

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::One
# B::One, C::AbstractScalar
*(B::AbstractScalar, C::One) = B
*(B::One, C::AbstractScalar) = C

# B::AbstractScalar, C::Zero
# B::Zero, C::AbstractScalar
*(B::AbstractScalar, C::Zero) = C
*(B::Zero, C::AbstractScalar) = B

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
*(B::AbstractScalar, x::Real) = Scalar{typeof(value(B))}(value(B) * x)
*(x::Real, B::AbstractScalar) = B * x

# B::AbstractScalar, v::Vector
# v::Vector, B::AbstractScalar
*(v::Vector{<:Real}, B::AbstractScalar) = Blade(value(B) * v)
*(B::AbstractScalar, v::Vector{<:Real}) = v * B

# ------ Specializations involving a Scalar instance

# B::Scalar, C::Scalar
*(B::Scalar, C::Scalar) = Scalar{typeof(value(B))}(value(B) * value(C))

# ------ Specializations involving a One instance

# B::One, C::One
*(B::One, C::One) = B

# B::One, C::Zero
# B::Zero, C::One
*(B::One, C::Zero) = C
*(B::Zero, C::One) = B

# ------ Specializations involving a Zero instance

# B::Zero, C::Zero
*(B::Zero, C::Zero) = B


text:
include("../src/ZZDiffusionBridge.jl")

using Random

function BM(ξ, L)
    u = 0.0
    v = randn()
    T = 1.0
    x =fs_expansion(ξ, u, v, L, T)
    t = range(0.0, stop=1, length=2<<(L) + 1)
    a = plot(leg=false)
    plot!(a, t, x, label = "random realization of X⁵", alpha = 0.5, color = :red)
    xaxis!(a, "t")
    scatter!(a, t, x, color= :black, markersize=2)
    return a
end

L =5
plot(BM(randn(2<<L-1), L))

function FSbasis(L)
    a = plot()
    plot!(a, 0:0.001:1, 0:0.001:1, label = "first $L levels basis functions", alpha = 0.5)
    for i in 0:L
        for k in 0:2^(i)-1
            plot!(a, 0:0.001:1, [Λ(t, i , k, 1.0) for t in 0:0.001:1], alpha = 0.5, leg=false)
            xaxis!(a, "t")
        end
    end
    return a
end

L = 6
Random.seed!(4)
ξ = randn(2^(L+1) - 1)
plot( FSbasis(L), BM(ξ, L), layout = (1,2))
savefig("../../output/bm_6levels.pdf")

using LinearAlgebra
using Plots
using Random

function zomming()
    dt = 10^(-4)
    t = 0.0:dt:100.0
    d = length(t)
    ν = 2.0
    μ = 0.0
    σ = 1.0
    N = randn(d-1)
    x = zeros(d)
    for i = 1:d-1
        x[i+1] = x[i] + ν*(μ -x[i])*dt  + σ*sqrt(dt)*N[i]
    end
    return (t, x)
end

function runall()
    tt, xx = zomming()
    f = length(tt)
    i = 1
    dt = 1000
    a = []
    for rep in 0:6
        id = i:dt:f
        push!(a, plot(tt[id], xx[id], color = :darkred, legend = false, xaxis = "t"))
        #, label = "zoom X$(2^rep)", legend =:bottomleft ))
        i = Int(ceil((i + f)*0.5))
        dt =  Int(ceil(dt/2))
    end
    return a
end

Random.seed!(5)
a = runall()

plot(plot(a[1], a[3], a[5], a[7], size = (750,500),
       layout = (2,2)))
savefig("../../output/bmzooming.pdf")

struct Normal
    μ::Float64
    σ::Float64
    v::Float64
    function Normal(μ, σ)
        new(μ, σ, inv(σ))
    end
end

function ZigZagSampler(N::Normal, ξ0::Float64, T::Float64)
    t = 0.0 ; θ = 1; ξ = ξ0
    Ξ = [ξ]
    Δ = [t]
    while (t<T)
        τ = sqrt.(-2*N.σ*log(rand()) + max(0, θ*(ξ - N.μ))^2) - θ*(ξ - N.μ)
        t = t + τ
        ξ = ξ + θ*τ
        push!(Ξ, ξ)
        push!(Δ, t)
        θ = -θ
    end
    return(Δ, Ξ)
end
function vel(x)
    θ = fill(0, length(x[1]))
    θ[1] = 1
    for i in 2:length(x[1])
        θ[i] = -θ[i-1]
    end
    return θ
end

Random.seed!(0)
x = ZigZagSampler(Normal(0.0,1.0), 0.0, 40.0)
θ = vel(x)
a1 = plot(x[1], x[2], leg= false,xaxis = "t", yaxis = "\\xi")
hline!(a1, [0.0], alpha = 0.4, color = :violet)
a2 =  plot(x[1], θ, linetype=:steppost, leg= false, xaxis = "t", yaxis = "\\theta")

plot( a1, a2, layout = (1,2), size = (500, 200))
savefig("../../output/zigzag1.pdf")

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
@testset "piecewisepolynomial" begin
    D = 2
    v1 = rand(D)
    v2 = rand(D)
    maha = Matrix{Float64}(I, D, D)
    degree = 3

    k = PiecewisePolynomialKernel(; degree=degree, dim=D)
    k2 = PiecewisePolynomialKernel{degree}(D)
    k3 = @test_deprecated PiecewisePolynomialKernel{degree}(maha)
    k4 = @test_deprecated PiecewisePolynomialKernel(; degree=degree, maha=maha)
    k5 = @test_deprecated PiecewisePolynomialKernel(; v=degree, dim=D)
    k6 = @test_deprecated PiecewisePolynomialKernel(; v=degree, maha=maha)

    @test k2(v1, v2) == k(v1, v2)
    @test k3(v1, v2) ≈ k(v1, v2)
    @test k4(v1, v2) ≈ k(v1, v2)
    @test k5(v1, v2) ≈ k(v1, v2)
    @test k6(v1, v2) ≈ k(v1, v2)

    @test_throws ErrorException PiecewisePolynomialKernel{4}(maha)
    @test_throws ErrorException PiecewisePolynomialKernel{4}(D)
    @test_throws ErrorException PiecewisePolynomialKernel{degree}(-1)
    @test_throws ErrorException PiecewisePolynomialKernel()
    @test_throws ErrorException PiecewisePolynomialKernel(; degree=degree)

    # default degree
    @test PiecewisePolynomialKernel(; dim=D) isa PiecewisePolynomialKernel{0}

    @test repr(k) ==
          "Piecewise Polynomial Kernel (degree = $(degree), ⌊dim/2⌋ = $(div(D, 2)))"

    # Standardised tests.
    TestUtils.test_interface(k, ColVecs{Float64}; dim_in=2)
    TestUtils.test_interface(k, RowVecs{Float64}; dim_in=2)
    test_ADs(() -> PiecewisePolynomialKernel{degree}(D))

    test_params(k, ())
end


text:
# test1dreadwrite.jl: Writing of 1D curves to Silo files

dbfile = Silo.siloopen(joinpath(dirname(@__FILE__),"test1dreadwrite.silo"), "w")
xin = collect(linspace(0, 1, 10))
yin = collect(linspace(1, 2, 10))

Silo.silowrite(dbfile, "Test_Curve", xin, yin)
# @show dbfile
Silo.siloclose(dbfile)
# @show dbfile

dbfile = Silo.siloopen(joinpath(dirname(@__FILE__),"test1dreadwrite.silo"), "r")
# @show Silo.names(dbfile)
x,y = Silo.siloread(dbfile, "Test_Curve")
Silo.siloclose(dbfile)

@test_approx_eq x xin
@test_approx_eq y yin

# @show x
# @show y


text:
"""
    SeisGain(d ; <keyword arguments>)

Gain a group of traces. Input and output are 2D.

# Arguments
- `d::Array{Real,2}`: two dimensional data.

# Keyword arguments
- `dt::Real=0.004`: sampling interval in secs.
- `kind::AbstractString="time"`: if kind="time", gain = t.^a . * exp(-bt);
                         if kind="agc", automatic gain control is applied.
- `coef::Vector{Real}=[2.0,0.0]`: if kind="time", coef = [a,b];
                                   if kind="agc", coef = [agc_gate]
- `normal::Int=0`: `normal=0` no normalization; `normal=1` normalize each trace by
                  amplitude; `normal=2` normalize each trace by rms value/

# Example
```julia
julia> using PyPlot
julia> d = SeisHypEvents();
       dout = SeisGain(d, kind="agc", coef=[0.05]);
       SeisPlotTX([d dout]);
```

Credits:  Aaron Staton
Updates: Juan I. Sabbione, Mauricio D. Sacchi, Fernanda Carozzi

"""
function SeisGain(d::Array{Td,2}; dt::Real=0.004, kind::AbstractString="time",
                    coef::Vector{Tp}=[2.0,0.0], normal::Int=0) where {Td<:Real,Tp<:Real}

    nt = size(d,1)
    nx = size(d,2)
    dout = zero(d)

    if kind == "time"   # Geometrical spreading-like gain

        a,b = coef
        t = collect(0:1:nt-1)*dt

        tgain = [(t[i]^a)*exp(b*t[i]) for i in 1:nt]
        dout = [d[i,k]*tgain[i] for i in 1:nt, k in 1:nx]

    end

    if kind=="agc"   #AGC

        L = floor(Int,round(coef[1]/(2dt)))
        h = triang(2L+1)

        for k = 1:nx
            e = [d[i,k]^2 for i in 1:nt]
            c = conv(promote(e,h)...)[L+1:nt+L]
            nc = length(c)
            rms = [sqrt(abs(c[i])) for i in 1:nc]
            epsi = 1.e-10*maximum(rms)
	    if epsi == 0
	       epsi = 0.0001
	    end
            op = [rms[i]/(rms[i]^2 + epsi) for i in 1:nc]
            dout[:,k] = d[:,k] .* op
        end
    end

    if normal==1     # Normalize by amplitude

        for k = 1:nx
            amax = maximum([abs(d[i,k]) for i in 1:nt])
	    if amax == 0
	       amax= 0.0001
	    end
            dout[:,k] = dout[:,k] ./ amax
        end

    end

    if normal==2;    # Normalize by rms

        for k = 1:nx
            amax = [sqrt(sum(d[i,k]^2)/nt) for i in 1:nt];
	    if amax ==0
	       amax = 0.0001
	    end

	    dout[:,k] = dout[:,k] ./ amax;
        end

    end

    return dout

end

"""
    SeisGain(in,out,parameters; <keyword arguments>)

    Gain a group of traces. Input and output are file names.

# Arguments
- `in::String`: Input file - Seis format
- `out::String`: Output file - Seis format.
- `parameters` : list of the keyword arguments for the function SeisGain.

# Keyword arguments
- `group="gather"` : Options are all, some or gather
- `key=["imx","imy"]` : Defines type of gather
- `itrace=1` : Initial trace number
- `ntrace=10000` : Total number of traces to process at once

# Example
```
julia> param = Dict(:dt=>0.01,:kind=>"time",:coef=>[2.0,0.0],:normal=>0)
julia> SeisGain(filein,fileout, param,group="all")
```

"""
function SeisGain(in::String,out::String,parameters;group="gather",key=["imx","imy"],itrace=1,ntrace=10000)

	kind = get(parameters,:kind,"time")
	coef = get(parameters,:coef,[2.0,0.0])
	normal = get(parameters,:normal,0)
	dt = get(parameters,:dt,0.004)

	ext = SeisMain.ReadTextHeader(in);
    	d1 = ext.d1
    	parameters = Dict(:dt=>dt,:kind=>kind,:coef=>coef,:normal=>normal)

	SeisProcessFile(in,out,[SeisGain],[parameters];group=group,key=key,itrace=itrace,ntrace=ntrace)
end

function triang(n::Integer)
    convert(Array{Float32,1},[1 - abs((k - (n-1)/2))/(n/2) for k in 0:(n-1)])
end


text:
isdefined(Base, :__precompile__) && __precompile__()

#__precompile__()

module StatLib

using Compat, DataArrays, GLM, DataFrames, Distributions, NLopt, Showoff, StatsBase, DataStructures, StatsFuns, JuMP, DBAPI, JavaCall, JDBC
using JSON, Requests, HttpParser

import DataArrays
import DBAPI, JavaCall, JDBC
import StatsBase: coef, coeftable, df, deviance, fit!, fitted, loglikelihood, model_response, nobs, vcov
import DataFrames: @~

using JSON, Compat, DataArrays, GLM, DataFrames, Distributions, NLopt, Showoff, StatsBase, DataStructures, StatsFuns, JuMP, MixedModels, NamedArrays, NLsolve
#using JSON, Requests, HttpParser

import DataArrays
#import DBAPI, JavaCall, JDBC
import StatsBase: coef, coeftable, df, deviance, fit!, fitted, loglikelihood, model_response, nobs, vcov
import Base: cond, std
import Distributions: Bernoulli, Binomial, Poisson, Gamma
import GLM: LogitLink, LogLink, InverseLink
import DataFrames: @~

export
       @~,
       rest

import Base: ==, *

include("rest.jl")

end # module

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
module Spirit
using Sockets
using Distributed

const Nullable{T} = Union{T,Nothing}
include("./constants.jl")
include("./HTTPError.jl")
include("./URI.jl")
include("./Headers.jl")
include("./SizedStream.jl")
include("./ChunkedTransferDecodeStream.jl")
include("./Connection.jl")
include("./Continuation.jl")
include("./Pipeline.jl")
include("./Request.jl")
# include("./Response.jl")
include("./Util.jl")
include("./HTTPServer.jl")
end # module Spirit


text:
using Documenter, DistStat

makedocs(
    format = Documenter.HTML(),
    sitename = "DistStat.jl",
    authors = "Seyoon Ko",
    clean = true,
    debug = true,
    pages = [
        "index.md"
    ]
)

deploydocs(
    repo   = "github.com/kose-y/DistStat.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
# Computes the dot product of a single column of a, specified by aColumnIdx, with the vector b.
# The number of elements used in the dot product determined by the size(A)[1].
# Note: bIdx is the last element of b used in the dot product.
function unsafedot( a::Matrix, aColIdx::Integer, b::Vector, bLastIdx::Integer )
    aLen     = size(a)[1]
    bBaseIdx = bLastIdx - aLen
    dotprod  = a[ 1, aColIdx ] * b[ bBaseIdx + 1 ]
    @simd for i in 2:aLen
        @inbounds dotprod += a[ i, aColIdx ] * b[ bBaseIdx + i ]
    end

    return dotprod
end

function unsafedot{T}( a::Matrix, aColIdx::Integer, b::Vector{T}, c::Vector{T}, cLastIdx::Integer )
    aLen = size(a)[1]
    bLen = length(b)
    bLen == aLen-1  || error( "length(b) must equal to length(a)[1] - 1" )
    cLastIdx < aLen || error( "cLastIdx but be < length(a)")

    dotprod = a[ 1, aColIdx] * b[ cLastIdx ]
    @simd for i in 2:aLen-cLastIdx
        @inbounds dotprod += a[ i, aColIdx ] * b[ i+cLastIdx-1 ]
    end
    @simd for i in 1:cLastIdx
        @inbounds dotprod += a[ aLen-cLastIdx+i, aColIdx ] * c[ i ]
    end

    return dotprod
end

function unsafedot( a::Vector, b::Vector, bLastIdx::Integer )
    aLen     = length(a)
    bBaseIdx = bLastIdx - aLen
    dotprod  = a[ 1 ] * b[ bBaseIdx + 1 ]
    @simd for i in 2:aLen
        @inbounds dotprod += a[ i ] * b[ bBaseIdx + i ]
    end

    return dotprod
end

function unsafedot{T}( a::Vector, b::Vector{T}, c::Vector{T}, cLastIdx::Integer )
    aLen    = length(a)
    dotprod = zero(a[1]*b[1])
    @simd for i in 1:aLen-cLastIdx
        @inbounds dotprod += a[ i ] * b[ i+cLastIdx-1 ]
    end
    @simd for i in 1:cLastIdx
        @inbounds dotprod += a[ aLen-cLastIdx+i ] * c[ i ]
    end

    return dotprod
end

# Shifts b into the end a.
# a = [ a, b ][1:length(a)]
function shiftin!{T}( a::Vector{T}, b::Vector{T} )
    aLen = length( a )
    bLen = length( b )

    if bLen >= aLen
        copy!( a, 1, b, bLen - aLen + 1, aLen )
    else

        for i in 1:aLen-bLen
            @inbounds a[i] = a[i+bLen]
        end
        bIdx = 1
        for i in aLen-bLen+1:aLen
            @inbounds a[i] = b[bIdx]
            bIdx += 1
        end
    end

    return a
end

function polyfit( y::AbstractVector, polyorder::Integer )
  A = [ x^p for x in 1:length(y), p = 0:polyorder ]
  Poly(A \ y)
end


text:
#
println("Dep) Test of the AutoIonization and PlasmaShift modules with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-AutoIonization-Plasma.sum")
setDefaults("method: continuum, Galerkin")
setDefaults("method: normalization, pure sine")   ## setDefaults("method: normalization, pure Coulomb")    setDefaults("method: normalization, pure sine")

grid = JAC.Radial.Grid(false)
grid = JAC.Radial.Grid(grid, rnt = 2.0e-6, h = 5.0e-2, hp = 1.0e-2, NoPoints = 900)
wa   = Atomic.Computation(Atomic.Computation(), name="xx", grid=grid, nuclearModel=Nuclear.Model(10.),
                          initialConfigs=[Configuration("1s 2s^2 2p^6")],
                          finalConfigs  =[Configuration("1s^2 2s^2 2p^4"), Configuration("1s^2 2s 2p^5"), Configuration("1s^2 2p^6")],
                          process = AugerInPlasma(),
                          processSettings = PlasmaShift.AugerSettings(PlasmaShift.DebyeHueckel(), 0.25, 0., 0, true, true, Tuple{Int64,Int64}[(1,0)]) )

wb = perform(wa)
setDefaults("print summary: close", "")


text:
include("testjacobi.jl")
include("testgauss.jl")


text:
module HampelFilter

using LinearAlgebra
using Statistics
export hampel
"""
    hampel(x, k, [,thr])

    Simple Hampel Filter for Julia lang.

    # Arguments
    - `x::n-elements Vector{Float64}`: time series data.
    - `k::Int`: the step size of a filter.
    - `thr::Float64`: Float threshold value (default: 3)

    # Examples

    ```jldoctest
        julia> x = [1, 2, 3, 10, 3, 2];
        julia> hampel(x,3)
        6-element Vector{Float64}:
        1.0
        2.0
        3.0
        2.5
        3.0
        2.0
    ```
"""

function hampel(x,k,thr=3)

    # Exceptions
    if length(x)==0
        throw(DomainError(x, "hampel() is not defined for zero-length vector"))
    elseif ~isinteger(k) || k<1
        throw(DomainError(k, "k should be positive integer."))
    elseif thr<=0
        throw(DomainError(thr, "thr should be positive."))
    end

    arrSize = length(x)
    output_x = Float64[]
    for i in 1:arrSize
        if i<2k
            kernel = collect(1:i)
        elseif i+2k>arrSize
            kernel = collect(i:arrSize)
        else
            kernel = collect(i:i+2k)
        end
        med = median(x[kernel])
        std = 1.4826*median(abs.(x[kernel].-med))
        if abs(x[i]-med)>thr*std
            push!(output_x,med)
        else
            push!(output_x,x[i])
        end
    end
    output_x
end
end


text:
using JuliaParser
using JuliaParser.LineNumbers
using Base.Test

code = "f() = 1 + 2 - 3\n"
@test SourceFile(code)[1:2] == Vector{UInt8}[Vector{UInt8}("f() = 1 + 2 - 3"), Vector{UInt8}("")]
data2 = Vector{UInt8}[Vector{UInt8}(code),[]]
LineBreaking(UInt64(1),SourceFile(code),data2)[7] = 'x'
@test String(data2[1]) == "f() = x + 2 - 3\n"
@test length(SourceFile("f() = 1 + 2 - 3")) == 1


text:
"""
initializations of Housedorff kernelalso clear up functions
"""
module HKernelInits
using CUDA, ..CUDAGpuUtils, Logging,StaticArrays
using ..HFUtils
export memoryAllocations,clearMainShmem,clearPadding,clearHalfOfPadding

"""
initialize constants
"""
function memoryAllocations()
    #for storing results
    #shmemGold,shmemSegm,shmemSum = createAndInitializeShmem(datablockDim, threadIdxX())
    resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
    #we need to set the values of the shared memory to 0 as it is not guaranteed by CUDA to be so
    clearMainShmem(resShmem,32)
    clearPadding(resShmem,32)# we separately clear padding
    #coordinates of data in main array
    #we will use this to establish weather we should mark  the data block as empty or full ...
    isMaskFull= zeros(MVector{1,Bool})
    isMaskEmpty= ones(MVector{1,Bool})
    #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
    locArr= zeros(MVector{32,Bool})

    return(resShmem,isMaskFull, isMaskEmpty,locArr)

end#initializeMemory

end#HKernelInits


text:
include("hits_tests.jl")
include("tools_tests.jl")
include("io_tests.jl")
include("calib_tests.jl")
include("math_tests.jl")
include("mc_tests.jl")

println("All tests complete")


text:
using Test
using DandelionWebSockets.WebSocketUpgrade
using DandelionWebSockets.WebSocketUpgrade: BadHTTPResponse, InvalidHTTPResponse, Request

todata(xs...) = codeunits(join(xs))

@testset "WebSocketUpgrade       " begin
    @testset "Requirement 3.1-1" begin
        # ## 3.1-1 MUST
        # The version of an HTTP message is indicated by an HTTP-Version field
        # in the first line of the message.
        # HTTP-Version   = "HTTP" "/" 1*DIGIT "." 1*DIGIT

        @testset "HTTP Version; The HTTP-Version field is not the first line; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "HTTP/1.1    101 Switch Protocols\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "HTTP Version; Version is HTTP/1.1; Major is 1 and minor is 1" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.major == 1
            @test response.httpversion.minor == 1
        end

        @testset "HTTP Version; Version is HTTP/12.34; Major is 12 and minor is 34" begin
            # Arrange
            responsetext = todata(
                "HTTP/12.34 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.major == 12
            @test response.httpversion.minor == 34
        end

        @testset "HTTP Version; The major version is not digits; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "HTTP/1abc2.34 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "HTTP Version; The minor version is not digits; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "HTTP/12.3abc4 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "HTTP Version; The line does not start with HTTP; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "NOTHTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "HTTP Version; The period is whitespace; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "HTTP/1 1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end
    end

    @testset "Requirement 3.1-2" begin
        @testset "HTTP Version; Version numbers are separate; Major is 1 and minor is 1" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.major == 1
            @test response.httpversion.minor == 1
        end
    end

    @testset "Requirement 3.1-3" begin
        @testset "HTTP Version; Version is HTTP/12.34; Major is 12 and minor is 34" begin
            # Arrange
            responsetext = todata(
                "HTTP/12.34 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.major == 12
            @test response.httpversion.minor == 34
        end
    end

    @testset "Requirement 3.1-4" begin
        @testset "HTTP Version; Version is HTTP/012.1; Major is 12" begin
            # Arrange
            responsetext = todata(
                "HTTP/012.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.major == 12
        end

        @testset "HTTP Version; Version is HTTP/000012.1; Major is 12" begin
            # Arrange
            responsetext = todata(
                "HTTP/000012.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.major == 12
        end

        @testset "HTTP Version; Version is HTTP/1.0000012; Minor is 12" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.0000012 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.httpversion.minor == 12
        end
    end

    @testset "Requirement 4.1-2" begin
        @testset "ResponseParser; A complete HTTP response without excess; The parser has a complete HTTP response" begin
            # Arrange
            response = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()

            # Act
            dataread(parser, response)

            # Assert
            @test hascompleteresponse(parser)
        end

        @testset "ResponseParser; No empty line; The parser does not have a complete response" begin
            # Arrange
            response = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
            )

            parser = ResponseParser()

            # Act
            dataread(parser, response)

            # Assert
            @test !hascompleteresponse(parser)
        end

        @testset "ResponseParser; A complete response with some excess data; The parser has a complete response" begin
            # Arrange
            response = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
                "Some excess data"
            )

            parser = ResponseParser()

            # Act
            dataread(parser, response)

            # Assert
            @test hascompleteresponse(parser)
        end

        @testset "ResponseParser; Complete response in two parts; The parser has a complete response" begin
            # Arrange
            response1 = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
            )
            response2 = todata(
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()

            # Act
            dataread(parser, response1)
            dataread(parser, response2)

            # Assert
            @test hascompleteresponse(parser)
        end

        @testset "ResponseParser; Complete response in two parts with excess; The parser has a complete response" begin
            # Arrange
            response1 = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
            )
            response2 = todata(
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
                "Some excess data"
            )

            parser = ResponseParser()

            # Act
            dataread(parser, response1)
            dataread(parser, response2)

            # Assert
            @test hascompleteresponse(parser)
        end
    end

    @testset "Requirement 6-1" begin
        # Requirement 6-1
        # After receiving and interpreting a request message, a server responds
        # with an HTTP response message.
        #
        #     Response      = Status-Line               ; Section 6.1
        #                     *(( general-header        ; Section 4.5
        #                     | response-header        ; Section 6.2
        #                     | entity-header ) CRLF)  ; Section 7.1
        #                     CRLF
        #                     [ message-body ]          ; Section 7.2
        # This section tests the validity of the Response grammar, including reading all headers.
        # Finding the boundary is tested by requirement 4.1-2.

        @testset "Header fields; First field is Date; Response has the field Date with the right date" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test findheaderfield(response, "Date") == "Sun, 06 Nov 1998 08:49:37 GMT"
        end

        @testset "Header fields; ETag field has value xyzzy; Response has the field ETag with value xyzzy" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "ETag: xyzzy\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test findheaderfield(response, "ETag") == "xyzzy"
        end

        @testset "Header fields; No ETag field; findheaderfield returns nothing" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test findheaderfield(response, "ETag") == nothing
        end

        @testset "Header fields; X-MyHeader with value Foo; Response has the field X-MyHeader with value Foo" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "ETag: xyzzy\r\n",
                "X-MyHeader: Foo\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test findheaderfield(response, "X-MyHeader") == "Foo"
        end
    end

    @testset "Requirement 6.1-1" begin
        # Requirement 6.1-1
        # The first line of a Response message is the Status-Line, consisting
        # of the protocol version followed by a numeric status code and its
        # associated textual phrase, with each element separated by SP
        # characters. No CR or LF is allowed except in the final CRLF sequence.
        #
        #     Status-Line = HTTP-Version SP Status-Code SP Reason-Phrase CRLF
        #
        # This section tests the validation of the Status-Line as a whole.
        # Requirement 6.1-2 tests the Status-Code.
        # Requirement 6.1-3 tests the Reason-Phrase.
        # Requirement 3.1-1 tests the HTTP Version field.
        # Requirement 4.1-2 tests the boundary between the message header and the body.
        #
        # NOTE: Requirement 3.1-1 deals only with the HTTP-Version field, both for requests and
        #       responses. The tests however deal only with responses. This should perhaps be
        #       restructured later on.
        @testset "Validation; Malformed Status Line; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "Malformed HTTP header",
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "Validation; Status code is not an integer; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 ABC OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "Validation; Reason Phrase is missing; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        @testset "Validation; HTTP Version is missing; BadHTTPResponse is thrown" begin
            # Arrange
            responsetext = todata(
                "200 OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            @test_throws BadHTTPResponse parseresponse(parser)
        end

        # I believe it is syntactically correct to have zero headers in the response.
        @testset "Validation; Status is 200 and no headers; Response is ok and status is 200" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.status == 200
        end
    end

    @testset "Requirement 6.1-2" begin
        @testset "ResponseParser; Response with status code 101; Parse result has status 101" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.status == 101
        end

        @testset "ResponseParser; Response with status code 200; Parse result has status 200" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "Date: Sun, 06 Nov 1994 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.status == 200
        end

        @testset "ResponseParser; Status code is 101 and year is 200x; Parse result has status 101" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 2008 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.status == 101
        end

        @testset "Status Line; Multiple spaces before status code 101; Status code is 101" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1    101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.status == 101
        end

        # TODO: Test that the status is exactly three digits
    end

    @testset "Requirement 6.1-3" begin
        @testset "Status Line; Reason Phrase is 'Switch Protocols' in header; Response has Reason Phrase 'Switch Protocols'" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 101 Switch Protocols\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.reasonphrase == "Switch Protocols"
        end

        @testset "Status Line; Reason Phrase is 'OK' in header; Response has Reason Phrase 'OK'" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200 OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.reasonphrase == "OK"
        end

        @testset "Status line; Multiple spaces before Reason Phrase; Spaces are trimmed" begin
            # Arrange
            responsetext = todata(
                "HTTP/1.1 200        OK\r\n",
                "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
                "\r\n",
            )

            parser = ResponseParser()
            dataread(parser, responsetext)

            # Act
            response = parseresponse(parser)

            # Assert
            @test response.reasonphrase == "OK"
        end
    end

    @testset "Requirement 6.2-1" begin
    end

    @testset "Requirement 14.10-3" begin
        # TODO Mark some requirements 14.10-* as inapplicable.
    end

    @testset "Requirement 14.42-5" begin
        # TODO Test that an Upgrade field can have more than one token
        # TODO Test that the order between tokens is maintained

        # @testset "Upgrade field; No Connection field; InvalidHTTPResponse is thrown" begin
        #     # Arrange
        #     responsetext = todata(
        #         "HTTP/1.1 200        OK\r\n",
        #         "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
        #         "Upgrade: websocket\r\n",
        #         "\r\n",
        #     )

        #     parser = ResponseParser()
        #     dataread(parser, responsetext)

        #     # Act
        #     @test_throws InvalidHTTPResponse parseresponse(parser)
        # end

        # @testset "Upgrade field; Connection field with value upgrade; InvalidHTTPResponse is not thrown" begin
        #     # Arrange
        #     responsetext = todata(
        #         "HTTP/1.1 200        OK\r\n",
        #         "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
        #         "Upgrade: websocket\r\n",
        #         "Connection: upgrade\r\n",
        #         "\r\n",
        #     )

        #     parser = ResponseParser()
        #     dataread(parser, responsetext)

        #     # Act
        #     # Tests that is does not throw an InvalidHTTPResponse
        #     parseresponse(parser)
        # end

        # @testset "Upgrade field; Connection field does not have token 'upgrade'; InvalidHTTPResponse is thrown" begin
        #     # Arrange
        #     responsetext = todata(
        #         "HTTP/1.1 200        OK\r\n",
        #         "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
        #         "Upgrade: websocket\r\n",
        #         "Connection: keep-alive\r\n",
        #         "\r\n",
        #     )

        #     parser = ResponseParser()
        #     dataread(parser, responsetext)

        #     # Act
        #     @test_throws InvalidHTTPResponse parseresponse(parser)
        # end

        # @testset "Upgrade field; Connection field with tokens upgrade and keep-alive; InvalidHTTPResponse is not thrown" begin
        #     # Arrange
        #     responsetext = todata(
        #         "HTTP/1.1 200        OK\r\n",
        #         "Date: Sun, 06 Nov 1998 08:49:37 GMT\r\n",
        #         "Upgrade: websocket\r\n",
        #         "Connection: keep-alive, upgrade\r\n",
        #         "\r\n",
        #     )

        #     parser = ResponseParser()
        #     dataread(parser, responsetext)

        #     # Act
        #     # Tests that is does not throw an InvalidHTTPResponse
        #     parseresponse(parser)
        # end

        # TODO Even when the Connection field has more than one token?
        # TODO Fail when the Connection field is missing the upgrade token
    end

    @testset "Good enough Request serialization" begin
        @testset "Request serialization; abs_path is /some/path; HTTP request has correct abs_path" begin
            # Arrange
            request = Request("example.org", "/some/path", [])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"^GET /some/path HTTP/1.1", String(take!(io))) != nothing
        end

        @testset "Request serialization; Path is /; Serialized path is /" begin
            # Arrange
            request = Request("example.org", "/", [])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"^GET / HTTP/1.1", String(take!(io))) !== nothing
        end

        @testset "Request serialization; Path is /; Line is ended by \\r\\n" begin
            # Arrange
            request = Request("example.org", "/", [])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"^GET / HTTP/1.1\r\n", String(take!(io))) !== nothing
        end

        @testset "Request serialization; Connection header is Upgrade; Header is listed in the request" begin
            # Arrange
            request = Request("example.org", "/", ["Connection" => "Upgrade"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"Connection: Upgrade\r\n", String(take!(io))) !== nothing
        end

        @testset "Request serialization; Sec-WebSocket-Key is abc; Header is listed in the request" begin
            # Arrange
            request = Request("example.org", "/", ["Sec-WebSocket-Key" => "abc"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"Sec-WebSocket-Key: abc\r\n", String(take!(io))) !== nothing
        end

        @testset "Request serialization; Sec-WebSocket-Key is def; Header is listed in the request" begin
            # Arrange
            request = Request("example.org", "/", ["Sec-WebSocket-Key" => "def"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"Sec-WebSocket-Key: def\r\n", String(take!(io))) !== nothing
        end

        @testset "Request serialization; Sec-WebSocket-Key is def; Sec-WebSocket-Key is not abc" begin
            # Arrange
            request = Request("example.org", "/", ["Sec-WebSocket-Key" => "def"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"Sec-WebSocket-Key: abc\r\n", String(take!(io))) === nothing
        end

        @testset "Request serialization; Sec-WebSocket-Key is def; Request ends with \\r\\n\\r\\n" begin
            # Arrange
            request = Request("example.org", "/", ["Sec-WebSocket-Key" => "def"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"\r\n\r\n$", String(take!(io))) !== nothing
        end

        @testset "Request serialization; Sec-WebSocket-Key is def; Request _only_ ends in two newlines" begin
            # Arrange
            request = Request("example.org", "/", ["Sec-WebSocket-Key" => "def"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"[^\r\n]\r\n\r\n$", String(take!(io))) !== nothing
        end

        @testset "Request serialization; No Host header provided; Host is written" begin
            # Arrange
            request = Request("example.org", "/", ["Sec-WebSocket-Key" => "def"])
            io = IOBuffer()

            # Act
            write(io, request)

            # Assert
            @test match(r"Host: example.org\r\n", String(take!(io))) !== nothing
        end
    end
end


text:
function toeplitz(c, r)
    @assert c[1] == r[1]
    nc = length(c)
    nr = length(r)
    A = zeros(eltype(c), nc, nr)
    A[:, 1] = c
    A[1, :] = r
    for i = 2:nr
        A[2:end, i] = A[1:end-1, i-1]
    end
    A
end

@inline obslength(y::AbstractMatrix) = size(y, 1)
@inline obslength(y::AbstractVector) = length(y[1])

# @inline Base.oftype(x::T, y::T) where T = y
@inline Base.oftype(x::AbstractMatrix, y::Vector{<:Vector}) = reduce(hcat, y)
@inline Base.oftype(x::AbstractMatrix, y::Vector{<:Number}) = reshape(y, 1, :)
@inline Base.oftype(x::Vector{<:Vector}, y::Matrix) = [y[:, i] for i = 1:size(y, 2)]
@inline Base.oftype(x::Vector{<:Vector}, y::Vector{<:Number}) =
    [[y[i]] for i in eachindex(y)]
@inline function Base.oftype(x::Vector{<:Number}, y::Matrix)
    size(y, 1) == 1 || size(y, 2) == 1 && return vec(y)
    throw(ArgumentError("Cannot convert a matrix with both dimensions greater than 1 to a vector."))
end

@inline function Base.oftype(x::Vector{<:Number}, y::Vector{<:Vector})
    size(y, 1) == 1 ||
        size(y, 2) == 1 ||
        throw(ArgumentError("Cannot convert a matrix with both dimensions greater than 1 to a vector."))
    reduce(vcat, y)
end

@inline Base.oftype(::Type{Matrix}, y::Vector{<:Number}) = reshape(y, 1, :)
@inline Base.oftype(::Type{Matrix}, y::Vector{<:Vector}) = reduce(hcat, y)
@inline Base.oftype(::Type{Matrix}, y::Matrix) = y
@inline Base.oftype(::Type{Matrix}, y::AbstractMatrix) = Matrix(y)

@inline time1(y::Vector) = y
@inline time1(y::Vector{<:Vector}) = reduce(hcat, y)'
@inline function time1(y::Matrix)
    if size(y, 1) == 1
        return vec(y)
    end
    Matrix(y')
end
@inline function time1(y::AbstractMatrix)
    if size(y, 1) == 1
        return vec(y)
    end
    Matrix(y')
end

@inline time2(y::Vector) = oftype(Matrix, y)
@inline time2(y::Vector{<:Vector}) = oftype(Matrix, y)
@inline time2(y::Matrix) = oftype(Matrix, y)
@inline time2(y::AbstractMatrix) = oftype(Matrix, y)

function Base.oftype(x::Number, y::AbstractArray)
    length(y) == 1 || throw(ArgumentError("Can only autoconvert one-element arrays"))
    y[]
end

rms(x::AbstractVector) = sqrt(mean(abs2, x))
sse(x::AbstractVector) = x ⋅ x
mse(x::AbstractVector) = sse(x) / length(x)

rms(x::AbstractMatrix) = sqrt.(mean(abs2.(x), dims = 2))[:]
sse(x::AbstractMatrix) = sum(abs2, x, dims = 2)[:]
modelfit(y, yh) = 100 * (1 .- rms(y .- yh) ./ rms(y .- mean(y, dims = 2)))
modelfit(y::T, yh::T) where {T<:AbstractVector} =
    100 * (1 .- rms(y .- yh) ./ rms(y .- mean(y)))
aic(x::AbstractVector, d) = log(sse(x)) .+ 2d / size(x, 2)
const nrmse = modelfit


text:
function foo(x)
    if x == 1 || x == 100
        return "hello"
    else
        return "goodbye"
    end
end


text:
abstract type Ode <: DiffEqFlux.FastLayer end

include("ctrnn0.jl")
include("ctrnn1.jl")
include("affine.jl")

""" Required as a FastLayer, define a paramsum field """
paramlength(f::Ode) = f.paramsum

""" Required as a FastLayer """
initial_params(f::Ode) = f.initial_params()

""" convert xu concatenation to x,u in 1d case """
@inline function extract_inputs(f::O, xu::AbstractArray{<:T,1}) where {T} where {O<:Ode}
    x = view(xu, 1:f.state_dim)
    u = view(xu, f.state_dim+1:f.state_dim+f.input_dim)
    return x, u
end

""" convert xu concatenation to x,u in 2d case """
@inline function extract_inputs(f::O, xu::AbstractArray{<:T,2}) where {T} where {O<:Ode}
    x = view(xu, 1:f.state_dim, :)
    u = view(xu, f.state_dim+1:f.state_dim+f.input_dim, :)
    return x, u
end

""" convert xu concatenation to x,u in 3d case """
@inline function extract_inputs(f::O, xu::AbstractArray{<:T,3}) where {T} where {O<:Ode}
    x = view(xu, 1:f.state_dim, :, :)
    u = view(xu, f.state_dim+1:f.state_dim+f.input_dim, :, :)
    return x, u
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
function baselinerotate(dimred)
    function rotate(Fₕ, Fₗ)
        m = project(Array(Fₕ), dimred)
        function antinan(F)
            # G = deepcopy(F)
            # G[mapslices(all, isnan.(F), dims=2)[:], :] .= 0.0
            if any(mapslices(all, isnan.(F), dims=2)[:])
                @warn "The test data has some all-NaN features. Replacing with epsilon noise."
                F[mapslices(all, isnan.(F), dims=2)[:], :] .= eps().*randn(sum(mapslices(all, isnan.(F), dims=2)), size(F, 2))
            end
            return F
        end
        F -> embed(m, Array(antinan(F)))
    end
end
export baselinerotate

""""""
function bPCA(Fₕ, Fₗ=zeros(size(Fₕ));
                    transform=identity,
                    prescale=identity,
                    prefilter=identity,
                    rotate=identity,
                    postscale=identity,
                    postfilter=identity)

    transform = transform(Fₕ, Fₗ)
    Fₗ = Fₗ |> transform
    Fₕ = Fₕ |> transform

    prefilter = prefilter(Fₕ, Fₗ)
    Fₗ = Fₗ |> prefilter
    Fₕ = Fₕ |> prefilter

    prescale = prescale(Fₕ, Fₗ)
    Fₗ = Fₗ |> prescale
    Fₕ = Fₕ |> prescale

    rotate = rotate(Fₕ, Fₗ)
    Fₗ = Fₗ |> rotate
    Fₕ = Fₕ |> rotate

    postfilter = postfilter(Fₕ, Fₗ)
    Fₗ = Fₗ |> postfilter
    Fₕ = Fₕ |> postfilter

    postscale = postscale(Fₕ, Fₗ)
    Fₗ = Fₗ |> postscale
    Fₕ = Fₕ |> postscale

    return F -> F |> transform |> prefilter |> prescale |> rotate |> postfilter |> postscale
end
export bPCA

function bPCA(; kwargs...)
    return (x, y) -> bPCA(x, y; kwargs...)
end


text:
module CSDP

using LinearAlgebra # For Diagonal
using SparseArrays # For SparseMatrixCSC

# Try to load the binary dependency
if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("CSDP not properly installed. Please run Pkg.build(\"CSDP\")")
end

export Blockmatrix

include("blockdiag.jl")
include("blockmat.h.jl")
include("blockmat.jl")
include("declarations.h.jl")
include("declarations.jl")
include("debug-mat.jl")
include("options.jl")
include("MOI_wrapper.jl")
include("MPB_wrapper.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Distributions, Random, StatsBase, DataFrames, Plots
	Random.seed!(1)
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.19"

# ╔═╡ 70179b62-1ea1-11eb-13ce-7fb4d54bad4b
function prn(lambda)
    k, p = 0, 1
    while p > MathConstants.e^(-lambda)
        k += 1
        p *= rand()
    end
    return k-1
end

# ╔═╡ 0cbc946e-2164-11eb-1b0c-7d4c7a607b1e
begin
	xGrid, lambda, N = 0:16, 5.5, 10^6

	pDist = Poisson(lambda)
	bPmf = pdf.(pDist,xGrid)
	data = counts([prn(lambda) for _ in 1:N],xGrid)/N
end

# ╔═╡ 0cbcbdf2-2164-11eb-193d-2baa8a6bce0b
begin
	plot( xGrid, data,
		line=:stem, marker=:circle,
		c=:blue, ms=10, msw=0, lw=4, label="MC estimate")
	plot!( xGrid, bPmf, line=:stem,
		marker=:xcross, c=:red, ms=6, msw=0, lw=2, label="PMF",
		ylims=(0,0.2), xlabel="x", ylabel="Probability of x events")
end

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.19"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═70179b62-1ea1-11eb-13ce-7fb4d54bad4b
# ╠═0cbc946e-2164-11eb-1b0c-7d4c7a607b1e
# ╠═0cbcbdf2-2164-11eb-193d-2baa8a6bce0b
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3


text:
using Memoize

@memoize function pattern(i::Int)
    vcat(zeros(Int,i), ones(Int,i), zeros(Int,i), fill(-1, i))
end

function pattern(position::Int, element::Int)
    p = pattern(position)
    return p[(element % length(p)) + 1]
end

function fft(input::String, phases::Int)
    m = [pattern(r, c) for r=1:length(input),c=1:length(input)]
    v = parse.(Int,split(input,""))

    for i=1:phases
        v = abs.(rem.(m * v, 10))
    end

    return v
end

examples = ["80871224585914546619083218645595", "19617804207202209144916044189917",
    "69317163492948606335995924319873"]

println("Day 16 Examples: ")
for example in examples
    println("$(example) becomes $(foldl(*, map(string, fft(example,100)[1:8])))")
end

input_file = open("input.txt")
input = readline(input_file)
close(input_file)

println("Day 16 Part 1: $(foldl(*, map(string, fft(input,100)[1:8])))")

# I never really understood the problem statement for part 2. From Reddit,
# it looks like I was not the only one.
# The real signal is the puzzle input repeated 10000 times.
signal = input ^ 10000

# The message offset is the first seven digits.
offset = parse(Int, signal[1:7])


text:
using TailwindCSS
using Test

"text-xl"  # Tailwind CLI should write this CSS into it's output

@testset "TailwindCSS.jl" begin
    input = joinpath(@__DIR__, "input.css")
    output = joinpath(@__DIR__, "output.css")
    config = Dict("content" => ["*.jl"])

    rm(output, force=true)

    TailwindCSS.minify(input, output, config)

    @test occursin("text-xl", read(output, String))
end


text:
import Base: ==, peek, last

mutable struct LinkedList
    value::Any
    next::Union{LinkedList,Nothing}
end

function (==)(x::LinkedList, y::LinkedList)
    while !isnothing(x.next) && !isnothing(y.next)
        if x.value != y.value
            return false
        end
        x = x.next
        y = y.next
    end
    return x.value == y.value && isnothing(x.next) && isnothing(y.next)
end

@assert LinkedList(1, LinkedList(2, nothing)) == LinkedList(1, LinkedList(2, nothing))
@assert LinkedList(1, LinkedList(2, nothing)) != LinkedList(1, LinkedList(3, nothing))
@assert LinkedList(1, LinkedList(2, nothing)) != LinkedList(2, LinkedList(2, nothing))

function peek(list::LinkedList, n::Integer)
    out = Array{Int}(undef, n)
    for i in 1:n
        out[i] = list.value
        list = list.next
    end
    return out
end

@assert all(peek(LinkedList(1, LinkedList(2, LinkedList(3, LinkedList(4, nothing)))), 2) .== [1, 2])

function insertafter!(list::LinkedList, index, items::Vector)
    head = list
    while head.value != index
        head = head.next
    end

    tail = head.next
    for x in items
        head.next = LinkedList(x, tail)
        head = head.next
    end

    return list
end

@assert insertafter!(LinkedList(1, LinkedList(2, LinkedList(3, nothing))), 2, [7, 8]) ==
    LinkedList(1, LinkedList(2, LinkedList(7, LinkedList(8, LinkedList(3, nothing)))))

function tolinkedlist(arr::Vector)
    list = LinkedList(arr[end], nothing)
    for x in reverse(arr)[2:end]
        list = LinkedList(x, list)
    end
    return list
end

@assert tolinkedlist([1, 2, 3]) == LinkedList(1, LinkedList(2, LinkedList(3, nothing)))

function last(list::LinkedList)
    head = list
    while !isnothing(head.next)
        head = head.next
    end
    return head
end

@assert last(LinkedList(1, LinkedList(2, LinkedList(3, nothing)))) == LinkedList(3, nothing)

max_stars_repo_name: twolodzko/advent-of-code

max_issues_repo_name: twolodzko/advent-of-code

max_forks_repo_name: twolodzko/advent-of-code


text:
fm = innerjoin( frshh, frspeople, on=[:data_year, :hid ] )

summarystats( fm[:tenure])

adddummies!( fm, hhv[:tentyp2], alias="tenure", use_value_as_label=true, missings_to_zero=true)
adddummies!( fm, hhv[:gvtregn], alias="region", use_value_as_label=false, missings_to_zero = true )

fm.deaf_blind=fm.registered_blind .| fm.registered_deaf .| fm.registered_partially_sighted
fm.yr = fm.data_year .- 2014
fm.any_dis = (
    fm.disability_vision .|
    fm.disability_hearing .|
    fm.disability_mobility .|
    fm.disability_dexterity .|
    fm.disability_learning .|
    fm.disability_memory .|
    fm.disability_other_difficulty .|
    fm.disability_mental_health .|
    fm.disability_stamina .|
    fm.disability_socially )

adddummies!(fm, hha[:heathad], alias="health_status" )
countmap(frsadult.HEATHAD)

print(countmap(frsadult.HEALTH1)) # long-standing illness
# Do you have any physical or mental health conditions or illnesses lasting or
# expected to last for 12 months or more?
print(countmap(frsadult.CONDIT)) # 'Whether condition limits day to day activities'
# Does your condition or illness/do any of your conditions or illnesses reduce your
# ability to carry-out day-to-day activities?

print(countmap(frsadult.DDATREP1)) # ever, but not now illness,disability limited activities
print(countmap(frsadult.DDAPROG1)) # ever, but not now: DDAPROG1 progressive illness 1 = adls down a little 2=lot

countmap(frsadult.LIMITL)
"For how long has your ability to carry-out day-to-day activities been reduced?"

# TODO fm = innerjoin( fm, frsadult, on=[:data_year, :hid ] )
make_enumerated_type( "Illness_Length", hha[:limitl ] )

println( make_enumerated_type( "Illness_Length", hha[:limitl ], true, true ))

println( make_enumerated_type( "ADLS_Inhibited", hha[:condit ], true, true ))

println( make_enumerated_type( "Long_Standing_Illness", hha[:health1 ], true, true ))

fm.adls_bad=fm.adls_are_reduced.==1
fm.adls_mid=fm.adls_are_reduced.==2
fm.rec_dla = ( fm.income_dlamobility.>0.0) .| ( fm.income_dlaself_care .>0.0 )
fm.rec_dla_care = ( fm.income_dlaself_care .>0.0 )
fm.rec_dla_mob = ( fm.income_dlamobility.>0.0 )
fm.rec_pip = ( fm.income_personal_independence_payment_mobility.>0.0) .| ( fm.income_personal_independence_payment_daily_living .>0.0 )
fm.rec_pip_care = ( fm.income_personal_independence_payment_daily_living .>0.0 )
fm.rec_pip_mob = ( fm.income_personal_independence_payment_mobility.>0.0)
fm.rec_esa = ( fm.income_employment_and_support_allowance.>0.0)
fm.rec_aa = ( fm.income_attendance_allowance.>0.0)
fm.rec_carers = ( fm.income_carers_allowance.>0.0)
fm_rec_aa = ( fm.income_attendance_allowance.>0.0)
fm.scotland = fm.region .== 299999999
fm.male = fm.sex .== 1

fm_working_age = fm[((fm.age .< 65).&(fm.from_child_record.!= 1)),:]
fm_pension_age = fm[(fm.age .>= 65),:]

dla_1=glm(@formula(rec_dla ~ male+age+age^2+has_long_standing_illness), fm, Binomial(), ProbitLink())
dla_2=glm(@formula(rec_dla ~ male+age+age^2+has_long_standing_illness+adls_are_reduced), fm, Binomial(), ProbitLink() )
dla_3=glm(@formula(rec_dla ~ male+age+age^2+has_long_standing_illness+adls_are_reduced+registered_blind), fm, Binomial(), ProbitLink() )
dla_4=glm(@formula(rec_dla ~ male+age+age^2+has_long_standing_illness+adls_are_reduced+registered_blind+registered_deaf), fm, Binomial(), ProbitLink() )
dla_5=glm(@formula(rec_dla ~ male+age+age^2+has_long_standing_illness+adls_bad+registered_blind+registered_deaf), fm, Binomial(), ProbitLink() )
dla_6=glm(
  @formula(
    rec_dla ~
      data_year+
      scotland+
      male+
      age+
      age^2+
      has_long_standing_illness+
      adls_bad+
      adls_mid+
      registered_blind+
      registered_partially_sighted+
      registered_deaf+
      disability_dexterity+
      disability_learning+
      disability_memory+
      disability_mental_health+
      disability_mobility+
      disability_socially+
      disability_stamina+
      disability_other_difficulty),
  fm,
  Binomial(),
  ProbitLink(),
  contrasts=Dict( :data_year=>DummyCoding()
  )
)

pred=StatsBase.predict(dla_1)

tab_dla = RegressionTables.regtable(dla_1,dla_2,dla_3,dla_4,dla_5,dla_6)

esa_1=glm(@formula(rec_esa ~ male+age+age^2+has_long_standing_illness), fm_working_age, Binomial(), ProbitLink())
esa_2=glm(@formula(rec_esa ~ male+age+age^2+has_long_standing_illness+adls_are_reduced), fm_working_age, Binomial(), ProbitLink() )
esa_3=glm(@formula(rec_esa ~ male+age+age^2+has_long_standing_illness+adls_are_reduced+registered_blind), fm_working_age, Binomial(), ProbitLink() )
esa_4=glm(@formula(rec_esa ~ male+age+age^2+has_long_standing_illness+adls_are_reduced+registered_blind+registered_deaf), fm_working_age, Binomial(), ProbitLink() )
esa_5=glm(@formula(rec_esa ~ male+age+age^2+has_long_standing_illness+adls_bad+registered_blind+registered_deaf), fm_working_age, Binomial(), ProbitLink() )
esa_6=glm(
  @formula(
    rec_esa ~
      data_year+
      scotland+
      male+
      age+
      age^2+
      has_long_standing_illness+
      adls_bad+
      adls_mid+
      registered_blind+
      registered_partially_sighted+
      registered_deaf+
      disability_dexterity+
      disability_learning+
      disability_memory+
      disability_mental_health+
      disability_mobility+
      disability_socially+
      disability_stamina+
      disability_other_difficulty),
  fm_working_age,
  Binomial(),
  ProbitLink(),
  contrasts=Dict( :data_year=>DummyCoding()
  )
)

pred=StatsBase.predict(esa_1)

tab_esa = RegressionTables.regtable(esa_1,esa_2,esa_3,esa_4,esa_5,esa_6)

pip_1=glm(@formula(rec_pip ~ male+age+age^2+has_long_standing_illness), fm_working_age, Binomial(), ProbitLink())
pip_2=glm(@formula(rec_pip ~ male+age+age^2+has_long_standing_illness+adls_are_reduced), fm_working_age, Binomial(), ProbitLink() )
pip_3=glm(@formula(rec_pip ~ male+age+age^2+has_long_standing_illness+adls_are_reduced+registered_blind), fm_working_age, Binomial(), ProbitLink() )
pip_4=glm(@formula(rec_pip ~ male+age+age^2+has_long_standing_illness+adls_are_reduced+registered_blind+registered_deaf), fm_working_age, Binomial(), ProbitLink() )
pip_5=glm(@formula(rec_pip ~ male+age+age^2+has_long_standing_illness+adls_bad+registered_blind+registered_deaf), fm_working_age, Binomial(), ProbitLink() )
pip_6=glm(
  @formula(
    rec_pip ~
      data_year+
      scotland+
      male+
      age+
      age^2+
      has_long_standing_illness+
      adls_bad+
      adls_mid+
      registered_blind+
      registered_partially_sighted+
      registered_deaf+
      disability_dexterity+
      disability_learning+
      disability_memory+
      disability_mental_health+
      disability_mobility+
      disability_socially+
      disability_stamina+
      disability_other_difficulty),
  fm_working_age,
  Binomial(),
  ProbitLink(),
  contrasts=Dict( :data_year=>DummyCoding()
  )
)

pred=StatsBase.predict(pip_1)

tab_pip = RegressionTables.regtable(pip_1,pip_2,pip_3,pip_4,pip_5,pip_6)

aa_6=glm(
  @formula(
    rec_aa ~
      data_year+
      scotland+
      male+
      age+
      age^2+
      has_long_standing_illness+
      adls_bad+
      adls_mid+
      registered_blind+
      registered_partially_sighted+
      registered_deaf+
      disability_dexterity+
      disability_learning+
      disability_memory+
      disability_mental_health+
      disability_mobility+
      disability_socially+
      disability_stamina+
      disability_other_difficulty),
  fm_pension_age,
  Binomial(),
  ProbitLink(),
  contrasts=Dict( :data_year=>DummyCoding()
  )
)

targets = Dict(
  :rec_aa=>fm_pension_age,
  :rec_pip=>fm_working_age,
  :rec_pip_mob=>fm_working_age,
  :rec_pip_care=>fm_working_age,
  :rec_esa=>fm_working_age,
  :rec_dla=>fm,
  :rec_dla_mob=>fm,
  :rec_dla_care=>fm )

main_out = []

for( target, dataset ) in targets

    push!(main_out, glm(
      (
        Term( target ) ~
          Term( :data_year )+
          Term( :scotland )+
          Term( :male )+
          Term( :age )+
          Term( :age )&Term(:age)+
          Term( :has_long_standing_illness )+
          Term( :adls_bad )+
          Term( :adls_mid )+
          Term( :registered_blind )+
          Term( :registered_partially_sighted )+
          Term( :registered_deaf )+
          Term( :disability_dexterity )+
          Term( :disability_learning )+
          Term( :disability_memory )+
          Term( :disability_mental_health )+
          Term( :disability_mobility )+
          Term( :disability_socially )+
          Term( :disability_stamina )+
          Term( :disability_other_difficulty )),
      dataset,
      Binomial(),
      ProbitLink(),
      contrasts=Dict( :data_year=>DummyCoding()
      )
    ))
end

RegressionTables.regtable( main_out ... )

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function _default(mesh::TOrSignal{M}, s::Style, data::Dict) where M <: GLNormalAttributeMesh
    @gen_defaults! data begin
        main = mesh
        boundingbox = const_lift(GLBoundingBox, mesh)
        shading = true
        color = nothing
        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "attribute_mesh.vert", "standard.frag",
            view = Dict("light_calc" => light_calc(shading))
        )
    end
end

function _default(mesh::TOrSignal{M}, s::Style, data::Dict) where M <: GLNormalMesh
    @gen_defaults! data begin
        shading = true
        main = mesh
        color = default(RGBA{Float32}, s)
        boundingbox = const_lift(GLBoundingBox, mesh)

        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "standard.vert", "standard.frag",
            view = Dict("light_calc" => light_calc(shading))
        )
    end
end
function _default(mesh::TOrSignal{M}, s::Style, data::Dict) where M <: GLNormalUVMesh
    @gen_defaults! data begin
        shading = true
        main = mesh
        color = default(RGBA{Float32}, s) => Texture
        boundingbox = const_lift(GLBoundingBox, mesh)

        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "uv_normal.vert", "standard.frag",
            view = Dict("light_calc" => light_calc(shading))
        )
    end
end
function _default(mesh::TOrSignal{M}, s::Style, data::Dict) where M <: GLNormalVertexcolorMesh
    @gen_defaults! data begin
        shading = true
        main = mesh
        boundingbox = const_lift(GLBoundingBox, mesh)
        color = nothing
        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "vertexcolor.vert", "standard.frag",
            view = Dict("light_calc" => light_calc(shading))
        )
    end
end

function _default(mesh::GLNormalColorMesh, s::Style, data::Dict)
    data[:color] = decompose(RGBA{Float32}, mesh)
    _default(GLNormalMesh(mesh), s, data)
end

function _default(main::TOrSignal{M}, ::style"grid", data::Dict) where M <: GLPlainMesh
    @gen_defaults! data begin
        primitive::GLPlainMesh = main
        color = default(RGBA, s, 1)
        bg_colorc = default(RGBA, s, 2)
        grid_thickness = Vec3f0(2)
        gridsteps = Vec3f0(5)
        boundingbox = const_lift(GLBoundingBox, mesh)
        shader = GLVisualizeShader("fragment_output.frag", "grid.vert", "grid.frag")
    end
end

function _default(mesh::TOrSignal{M}, s::Style, data::Dict) where M <: GLPlainMesh
    @gen_defaults! data begin
        primitive::GLPlainMesh = mesh
        color = default(RGBA, s, 1)
        boundingbox = const_lift(GLBoundingBox, mesh)
        shader = GLVisualizeShader("fragment_output.frag", "plain.vert", "plain.frag")
    end
end