text:
#using BenchmarkTools, Roots, Distributions
# Input from user

function input(prompt::String = "")::String
    print(prompt)
    return chomp(readline())
end

function fdesc(r,t,compos = 'c')
    if compos == 'c'
        czero = 1.*exp.(-r.*t)
    else
        czero = 1./((r.+1).^t)
    end
    return czero
end

function actualz(fd,t,compos='c')
    if compos == 'c'
        re = - log.(fd)./t
    else
        re = fd.^(-1./t).- 1
    end
    return re
end

function bootstrap(p, c, m,compos='c')
    n = length(c)
    cfMat = zeros(n,n)
    for i = 1:n, j = 1:i #fill the cfMat matrix
        cfMat[i,j] = cfMat[i,j] + c[i]
        if i == j #along diagonal, add face value (1)
            cfMat[i,j] = cfMat[i,j] + 1
        end
    end
    B = cfMat\p
    return B, actualz.(B,m, compos)
end

function nss(m,b0,b1,b2,b3, tau, tau2)
    f =  b0 + b1*exp.(-m/tau) + b2*(m/tau).*exp.(-m/tau) + b3*(m/tau2).*exp.(-m/tau2)
    s =  b0 + b1*(1 - exp.(-m/tau))./(m/tau) + b2*((1 - exp.(-m/tau)) ./(m/tau) - exp.(-m/tau)) + b3*((1 - exp.(-m/tau2))./(m/tau2) - exp.(-m/tau2))
    d = exp.(-s.*m)
    return s,f,d
end

function CouponPago(DataReferencia,Mes,Dia)
    if DataReferencia >= Date(Year(DataReferencia),Mes,Dia)
        return true
    else
        return false
    end
end

function DiasCorridos(DataReferencia,Mes,Dia)
    DataCouponPago = Date(Year(DataReferencia),Mes,Dia)
    if DataReferencia >= DataCouponPago
        DiasJurosCorridos = DataReferencia - DataCouponPago
    else
        DataCouponPago = Date(Year(DataReferencia)-Year(1),Mes,Dia)
        DiasJurosCorridos = DataReferencia - DataCouponPago
    end
    return Dates.value(DiasJurosCorridos), DataCouponPago
end

function JCorridos(dados,DataReferencia)
    NEmissões = size(dados)[1]
    a = []
    for i = 1:NEmissões
        Mes = Month(Date(dados[i,2]))
        Dia = Day(Date(dados[i,2]))
        DiasJurosCorridos = DiasCorridos(DataReferencia,Mes,Dia)[1]
        JurosCorridos = DiasJurosCorridos / 365 * dados[i,3] * 100
        push!(a,dados[i,1])
        push!(a,DiasCorridos(DataReferencia,Mes,Dia)[2])
        push!(a,DiasJurosCorridos)
        push!(a,dados[i,3])
        push!(a,JurosCorridos)
        push!(a,dados[i,4])
        push!(a,JurosCorridos + dados[i,4])
    end
    a = reshape(a,7,:)
    return JuliaDB.table(a[1,:], a[2,:], a[3,:], a[4,:], a[5,:], a[6,:],a[7,:], names=[:ISIN, :UCouponPago,:DCorridos, :TxCoupon, :JCorridos, :PLimpo, :PLiquid], pkey=(:ISIN))
end

function FluxosCaixa(dados,DataReferencia)
    NEmissões = size(dados)[1]
    b=[]
    for i = 1:NEmissões
        Mes = Month(Date(dados[i,2]))
        Dia = Day(Date(dados[i,2]))
        if CouponPago(DataReferencia,Mes,Dia) == true
            inicio = Dates.value(Year(DataReferencia)) + 1
        else
            inicio = Dates.value(Year(DataReferencia))
        end
        fim = Dates.value(Year(Date(dados[i,2])))
        for j = inicio:fim
            push!(b,dados[i,1])
            DataCF = tobday(:TARGET, Date(j, Mes, Dia))
            push!(b,DataCF)
            tenor = Dates.value(DataCF - DataReferencia)/365
            push!(b,tenor)
            if j == fim
                push!(b, (1 + dados[i,3])*100)
            else
                push!(b, dados[i,3]*100)
            end
        end
    end
    b = reshape(b,4,:)
    return JuliaDB.table(b[1,:], b[2,:], b[3,:], b[4,:],names=[:ISIN, :Data, :Tenores, :FluxosCaixa], pkey=(:Data, :ISIN))
end

function fdesconto(nss, tenores, janela_obs = 1:400)
    nobs = length(nss) # Número datas disponíveis para os parâmetros NSS
    NTenores = length(tenores)
    a = []
    nobs_plus = nobs + 1 # Se janela_obs começar em 1 então ordem obs = nobs
    for i = janela_obs
        obs = nobs_plus - i
        push!(a, nss[obs][1])
        b0 = nss[obs][2]
        b1 = nss[obs][3]
        b2 = nss[obs][4]
        b3 = nss[obs][5]
        tau = nss[obs][6]
        tau2 = nss[obs][7]
        for m in tenores
            spot =  b0 + b1*(1 - exp(-m/tau))/(m/tau) + b2*((1 - exp(-m/tau)) /(m/tau) - exp(-m/tau)) + b3*((1 - exp(-m/tau2))/(m/tau2) - exp(-m/tau2))
            desc = exp(-spot/100*m)
            push!(a, desc)
        end
    end
    fdes = reshape(a, NTenores + 1, :)
    return permutedims(fdes, [2,1]) # Transposta: fdesc' não funciona pois temos datas
end

function fspot(nss, tenores, janela_obs = 1:400)
    nobs = length(nss) # Número datas disponíveis para os parâmetros NSS
    NTenores = length(tenores)
    a = []
    nobs_plus = nobs + 1 # Se janela_obs começar em 1 então ordem obs = nobs
    for i = janela_obs
        obs = nobs_plus - i
        push!(a, nss[obs][1])
        b0 = nss[obs][2]
        b1 = nss[obs][3]
        b2 = nss[obs][4]
        b3 = nss[obs][5]
        tau = nss[obs][6]
        tau2 = nss[obs][7]
        for m in tenores
            spot =  b0 + b1*(1 - exp(-m/tau))/(m/tau) + b2*((1 - exp(-m/tau)) /(m/tau) - exp(-m/tau)) + b3*((1 - exp(-m/tau2))/(m/tau2) - exp(-m/tau2))
            push!(a, spot)
        end
    end
    funspot = reshape(a, NTenores + 1,:)
    return permutedims(funspot, [2,1]) # Transposta: funspot' não funciona pois temos datas
end

function fcaixa(tcoupon,freq, mat, capital, amort=false)
    numcf = convert(Int64,ceil(freq * mat))   #
    t1 = mod(mat, 1/freq)
    tVec = fill(t1,numcf)
    for i = 2:numcf
        tVec[i] = tVec[i-1] + 1/freq
    end
    if amort == false
        coupon = (tcoupon/freq) * capital
        cfVec = fill(coupon, numcf)
        cfVec[numcf] = coupon + capital
    else
        cfVec = zeros(numcf)
        for i = 1:numcf
            coupon = (tcoupon/freq) * (capital - (i-1)*capital/numcf)
            cfVec[i] = coupon + capital / numcf
        end
    end
    return cfVec, tVec
end

function valobrig(cf,r,t,compos='c')
    if compos == 'c'
        cdisc = cf.*exp.(-r.*t)
        P = sum(cdisc)
    else
        cdisc = cf./((r.+1).^t)
        P = sum(cdisc)
    end
    return P
 end

function ytm(p,cf,t)
    ytm = fzero(y->valobrig(cf,y,t)-p,0.0)
    return ytm
end

function duration(preço,cflow,m,ytm,compos='c')
    if compos == 'c'
       cdisc = m.*cflow .*exp.(- ytm .* m)
       dMod = dMac = sum(cdisc)/preço
    else
        cdisc = cflow.*m./((1+ytm).^(m.+1))
        D = sum(cdisc) # Duration, row vector
        dMod = D/preço # Modified duration
        dMac = D*(1+ytm)/preço # Macaulays duration
    end
    return dMod, dMac
end

function VaR(alfa,sigma, portfolio)
    var = -quantile(Normal(0,1),alfa) * sigma * portfolio
    return var
end

function ETL(alfa, sigma, portfolio)
    z = quantile(Normal(0,1),alfa)
    etl = pdfNormal(z)/alfa * sigma * portfolio
    return etl
end

function nss(m,svensson)
    datas = select(svensson,:Data)
    nlinhas = length(datas)
    b0 = select(svensson,:BETA0)
    b1 = select(svensson,:BETA1)
    b2 = select(svensson,:BETA2)
    b3 = select(svensson,:BETA3)
    tau = select(svensson,:TAU1)
    tau2 = select(svensson,:TAU2)
    s = Array{Float64}(nlinhas)
    f = Array{Float64}(nlinhas)
    d = Array{Float64}(nlinhas)
    for i = 1:nlinhas
        # Taxas spot
        s[i] =  b0[i] + b1[i]*(1 - exp(-m/tau[i]))/(m/tau[i]) + b2[i]*((1 - exp(-m/tau[i])) /(m/tau[i]) - exp(-m/tau[i])) + b3[i]*((1 - exp(-m/tau2[i]))/(m/tau2[i]) - exp(-m/tau2[i]))
        # Taxas forward
        f[i] =  b0[i] + b1[i]*exp(-m/tau[i]) + b2[i]*(m/tau[i])*exp(-m/tau[i]) + b3[i]*(m/tau2[i])*exp(-m/tau2[i])
        # Função de Desconto
        d[i] = exp(-s[i]*m/100)
    end
    tabela = JuliaDB.table(datas,s,f,d,names=[:Data, :Spot, :Forward, :Desconto])
    return tabela
end

function bucketing(vk, v1,sigma1, v2, sigma2, corr12)
    cov12 = sigma1*sigma2*corr12
    sigmak = sigma1 + (vk-v1)/(v2-v1)*(sigma2-sigma1)
    a = sigma1^2+sigma2^2-2*cov12
    b = 2*(-sigma2^2+cov12)
    c = sigma2^2-sigmak^2
    raiz2 = (-b - sqrt(b^2-4*a*c))/(2*a)
    raiz1 = (-b + sqrt(b^2-4*a*c))/(2*a)
    minraiz = min(raiz1,raiz2)
    println(a,"\n",b,"\n",c)
    if  minraiz < 0
        return max(raiz1,raiz2)
    else
        return minraiz
    end
end

max_stars_repo_name: iscalprog/Iscalc

max_issues_repo_name: iscalprog/Iscalc

max_forks_repo_name: iscalprog/Iscalc


text:
ROOT = pwd()
include(joinpath(ROOT, "src_PowSysMod", "PowSysMod_body.jl"))
include(joinpath("func_definitions.jl"))

function main()

    ########################################
    # Construction du problème type
    OPFpbs = load_OPFproblems(MatpowerInput, joinpath("data_Matpower", "matpower", "WB2.m"))
    problem = build_globalpb!(OPFpbs)

    println("WB2 problem built")
    # print(problem)

    ########################################
    # Setting relaxation order by constraint
    d_i = Dict("cstr1"=>1, "cstr2"=>1)
    k_i = get_ki(problem)
    check_di_ki!(d_i, k_i)

    ########################################
    # Construction du sparsity pattern
    sparsity_pattern = compute_sparsitypattern(problem, d_i, k_i)

    # Extension chordale et détection des cliques maximales
    compute_chordalextension!(sparsity_pattern)
    max_cliques = compute_maxcliques(sparsity_pattern)

    ########################################
    # Relaxation degree par clique and variables par constrainte
    varsbycstr = compute_varsbycstr(problem)
    cliquevarsbycstr = compute_varsbycstr(sparsity_pattern, max_cliques, varsbycstr)
    orderbyclique = compute_cliqueorders(sparsity_pattern, d_i, varsbycstr, max_cliques)

    ########################################
    # Calcul des matrices B_i et pose du probleme
    B_i = compute_Bibycstr(problem, max_cliques, cliquevarsbycstr, orderbyclique)
    SDP_SOS = build_SDP_SOS(problem, d_i, max_cliques, B_i, cliquevarsbycstr, orderbyclique)

    ########################################
    # Calcul d'une solution par un solveur
    m = make_JuMPproblem(SDP_SOS)

end

main()

ext: jl

max_stars_repo_name: Gimaju/GOC

max_issues_repo_name: Gimaju/GOC

max_forks_repo_name: Gimaju/GOC


text:
using Plasma
using Plots

TD = 0.3
Te = 0.1
D = species.D
e = species.e

D_D = Distribution(Maxwellian(TD, D.m), D)
D_e = Distribution(Maxwellian(Te, e.m), e)
G = Geometry()

plasma = ElectrostaticPlasma([D_D, D_e], G)

sol = Plasma.solve(plasma, dim=1, GPU=false)

Plasma.plot(sol)

## Two-stream instability

e = species.e

function TwoStream(vth2, vs1, vs2)

    function P(x,v)
        if !(v isa Array)
            v = [v]
        end

        if !(x isa Array)
            x = [x]
        end

        v = sqrt(sum(v .^2))
        x = sqrt(sum(x .^2))

        0.5/sqrt(vth2 * π) * exp(-(v-vs1)*(v-vs1)/vth2) + 0.5/sqrt(vth2 * π) * exp(-(v-vs2)*(v-vs2)/vth2) * (1+0.02*cos(3*π*x))
    end
end

D_e = Distribution(TwoStream(0.02, 1.6, -1.4), e)
G = Geometry()

plasma = ElectrostaticPlasma([D_e], G)

sol = Plasma.solve(plasma, dim=1, GPU=false, time_ub = 100.0, ub=2.0)

## 2D with custom P and species
using Plasma

Tα = 70000 # eV

α = Species(1.602176634e-19, 6.6446562e-27)

function HotCarrier(T)
    Kb = 8.617333262145e-5
    P(x,v) = exp(-v/(Kb*T))
end

Dα = Distribution(HotCarrier(Tα), α)
G = Geometry() # TODO define a custom geometry

plasma = ElectrostaticPlasma([Dα], G)

Plasma.solve(plasma, dim=2) # with GPU

Plasma.plot(sol)

## 3D CollisionlessPlasma
using Plasma

TD = 15000 # eV
TT = 15000 # eV
Te = 13000 # eV

e = species.e
T = species.T
D = species.D

De = Distribution(Maxwellian(Te, e.m), e)
DT = Distribution(Maxwellian(TT, T.m), T)
DD = Distribution(Maxwellian(TD, D.m), D)
G = Geometry()

plasma = CollisionlessPlasma([De,DT,DD], G)

Plasma.solve(plasma)

Plasma.plot(sol)


text:
# Biryani -- A conversion and validation toolbox
# By: Emmanuel Raviart <emmanuel@raviart.com>
#
# Copyright (C) 2015 Emmanuel Raviart
# https://github.com/eraviart/Biryani.jl
#
# This file is part of Biryani.
#
# 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.

function input_to_slug(convertible::Convertible; separator::Char = '-', transform::Function = lowercase)
  """Simplify a string to a slug (ie a string containing only 0-9, A-Z, a-z & separators)."""
  if convertible.error !== nothing || convertible.value === nothing
    return convertible
  end
  slug = slugify(convertible.value, separator = separator, transform = transform)
  return Convertible(isempty(slug)? nothing : slug , convertible.context)
end

input_to_slug(; separator::Char = '-', transform::Function = lowercase) = convertible::Convertible -> input_to_slug(
  convertible, separator = separator, transform = transform)


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module Images

import Base.Iterators.take
import Base: +, -, *
import Base: abs, atan, clamp, convert, copy, copy!, delete!,
             eltype, get, getindex, haskey, hypot,
             imag, length, map, map!, maximum,
             minimum, ndims, one, parent, permutedims, real, reinterpret,
             reshape, resize!,
             setindex!, show, similar, size, sqrt,
             strides, sum, write, zero

export float32, float64
export HomogeneousPoint

using Base: depwarn
using Base.Order: Ordering, ForwardOrdering, ReverseOrdering

using StaticArrays
using Base64: Base64EncodePipe

# CHECKME: use this or follow deprecation and substitute?
using SparseArrays: findnz

# "deprecated imports" are below

using Reexport
@reexport using ImageCore
if isdefined(ImageCore, :permuteddimsview)
    export permuteddimsview
end

using FileIO
export load, save
import .Colors: Fractional
import Graphics
import Graphics: width, height, Point
using StatsBase  # TODO: eliminate this dependency
using IndirectArrays, ImageCore.MappedArrays

# TODO: can we get rid of these definitions?
const NumberLike = Union{Number,AbstractGray}
const RealLike = Union{Real,AbstractGray}

const is_little_endian = ENDIAN_BOM == 0x04030201

@reexport using ImageTransformations
@reexport using ImageAxes
@reexport using ImageMetadata
@reexport using ImageFiltering
@reexport using ImageMorphology
@reexport using ImageDistances
@reexport using ImageContrastAdjustment
@reexport using ImageQualityIndexes

if isdefined(ImageQualityIndexes, :assess_psnr)
    # deprecated since ImageQualityIndexes v0.1.4
    Base.@deprecate_binding psnr assess_psnr
    Base.@deprecate_binding ssim assess_ssim
else
    const psnr = ImageQualityIndexes.psnr
    export psnr
end

# Both ImageMetadata v0.9.0 and ImageAxes v0.6.0 deprecate the symbol data and
# this causes a name conflict
if isdefined(ImageMetadata, :arraydata)
    Base.@deprecate_binding data arraydata
end
# ImageMetadata < v0.9.0 compatibility
if !hasmethod(arraydata, (ImageMeta, ) )
    ImageAxes.arraydata(img::ImageMeta) = ImageMetadata.data(img)
end

# Non-exported symbol bindings to ImageShow
import ImageShow
if isdefined(ImageShow, :play)
    @doc (@doc ImageShow.play)
    const play = ImageShow.play
else
    play(args...; kwargs...) = error("The `Images.play` function requires ImageShow at least 0.3.0.")
end
if isdefined(ImageShow, :explore)
    @doc (@doc ImageShow.explore)
    const explore = ImageShow.explore
else
    explore(args...; kwargs...) = error("The `Images.explore` function requires ImageShow at least 0.3.0.")
end
if isdefined(ImageShow, :gif)
    @doc (@doc ImageShow.gif)
    const gif = ImageShow.gif
else
    gif(args...; kwargs...) = error("The `Images.gif` function requires ImageShow at least 0.3.0.")
end

# While we are bridging the old API and the new API in ImageContrastAdjustment
# we need to import these functions because we make new definitions for them
# in deprecations.jl
import ImageContrastAdjustment: build_histogram, adjust_histogram, adjust_histogram!

using ImageMetadata: ImageMetaAxis
import ImageMorphology: dilate, erode
import ImageTransformations: restrict
using TiledIteration: EdgeIterator

include("compat.jl")
include("misc.jl")
include("labeledarrays.jl")
include("algorithms.jl")
include("exposure.jl")
include("deprecations.jl")
include("corner.jl")
include("edge.jl")

export
    # types
    BlobLoG,
    ColorizedArray,
    Percentile,

    # macros
    @test_approx_eq_sigma_eps,

    # core functions
    assert_timedim_last,
    colordim,
    coords_spatial,
    copyproperties,
    data,
    height,
    maxabsfinite,
    maxfinite,
    minfinite,
    nimages,
    pixelspacing,
    properties,
    restrict,
    sdims,
    size_spatial,
    shareproperties,
    spacedirections,
    spatialorder,
    spatialproperties,
    timedim,
    width,
    widthheight,

    # algorithms
    backdiffx,
    backdiffy,
    forwarddiffx,
    forwarddiffy,
    imcorner,
    imcorner_subpixel,
    corner2subpixel,
    harris,
    shi_tomasi,
    kitchen_rosenfeld,
    fastcorners,
    meancovs,
    gammacovs,
    imedge,  # TODO: deprecate?
    blob_LoG,
    findlocalmaxima,
    findlocalminima,
    imgaussiannoise,
    imlineardiffusion,
    imROF,
    otsu_threshold,
    yen_threshold,

    #Exposure
    complement,
    imhist,
    histeq,
    adjust_gamma,
    histmatch,
    clahe,
    imadjustintensity,
    imstretch,
    cliphist,

    magnitude,
    magnitude_phase,
    meanfinite,
    entropy,
    orientation,
    padarray,
    phase,
    thin_edges,
    thin_edges_subpix,
    thin_edges_nonmaxsup,
    thin_edges_nonmaxsup_subpix,
    canny,
    integral_image,
    boxdiff,
    bilinear_interpolation,
    gaussian_pyramid,

    # phantoms
    shepp_logan

_length(A::AbstractArray) = length(eachindex(A))
_length(A) = length(A)

"""
Constructors, conversions, and traits:

    - Construction: use constructors of specialized packages, e.g., `AxisArray`, `ImageMeta`, etc.
    - "Conversion": `colorview`, `channelview`, `rawview`, `normedview`, `permuteddimsview`
    - Traits: `pixelspacing`, `sdims`, `timeaxis`, `timedim`, `spacedirections`

Contrast/coloration:

    - `clamp01`, `clamp01nan`, `scaleminmax`, `colorsigned`, `scalesigned`

Algorithms:

    - Reductions: `maxfinite`, `maxabsfinite`, `minfinite`, `meanfinite`, `integral_image`, `boxdiff`, `gaussian_pyramid`
    - Resizing: `restrict`, `imresize` (not yet exported)
    - Filtering: `imfilter`, `imfilter!`, `mapwindow`, `imROF`, `padarray`
    - Filtering kernels: `Kernel.` or `KernelFactors.`, followed by `ando[345]`, `guassian2d`, `imaverage`, `imdog`, `imlaplacian`, `prewitt`, `sobel`
    - Exposure : `imhist`, `histeq`, `adjust_gamma`, `histmatch`, `imadjustintensity`, `imstretch`, `imcomplement`, `clahe`, `cliphist`
    - Gradients: `backdiffx`, `backdiffy`, `forwarddiffx`, `forwarddiffy`, `imgradients`
    - Edge detection: `imedge`, `imgradients`, `thin_edges`, `magnitude`, `phase`, `magnitudephase`, `orientation`, `canny`
    - Corner detection: `imcorner`,`imcorner_subpixel`, `harris`, `shi_tomasi`, `kitchen_rosenfeld`, `meancovs`, `gammacovs`, `fastcorners`
    - Blob detection: `blob_LoG`, `findlocalmaxima`, `findlocalminima`
    - Morphological operations: `dilate`, `erode`, `closing`, `opening`, `tophat`, `bothat`, `morphogradient`, `morpholaplace`, `feature_transform`, `distance_transform`, `convexhull`
    - Connected components: `label_components`, `component_boxes`, `component_lengths`, `component_indices`, `component_subscripts`, `component_centroids`
    - Interpolation: `bilinear_interpolation`

Test images and phantoms (see also TestImages.jl):

    - `shepp_logan`
"""
Images

end


text:
function avgfilter2d_setup(n, w = 3)
    xs = rand(SVector{w,UInt8}, n, n)
    ys = zero(xs)
    return (; ys, xs)
end

# Not sure if it's relevant but intel's benchmark convert `unsigned char` to
# `unsigned int` during the computation:
# https://github.com/neboat/cilkbench/blob/063cba43f2b6276c913d3ad87c7ad3257fd4f813/intel/AveragingFilter_01_07_15/src/AveragingFilter.cpp#L61
@inline u32(x::SVector{S}) where {S} = SVector{S,UInt32}(x)
@inline u8(x::SVector{S}) where {S} = x .% UInt8

#! format: off
@inline function avg3x3(xs0, i, j)
    @static if USE_ALIASSCOPE
        xs = Const(xs0)
    else
        xs = xs0
    end
    y = @inbounds (
        u32(xs[i-1,j-1]) .+ u32(xs[i,j-1]) .+ u32(xs[i+1,j-1]) .+
        u32(xs[i-1,j  ]) .+ u32(xs[i,j  ]) .+ u32(xs[i+1,j  ]) .+
        u32(xs[i-1,j+1]) .+ u32(xs[i,j+1]) .+ u32(xs[i+1,j+1])
    ) .÷ 9
    return u8(y)
end
#! format: on

function avgfilter2d_seq!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = size(xs)[1] - 2
    m = size(xs)[2] - 2
    @maybe_aliasscope begin
        for j0 in 0:m-1
            j = firstindex(xs, 2) + j0 + 1
            for i0 in 0:n-1
                i = firstindex(xs, 1) + i0 + 1
                @inbounds ys[i, j] = avg3x3(xs, i, j)
            end
        end
    end
    return ys
end

function avgfilter2d_threads!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = size(xs)[1] - 2
    m = size(xs)[2] - 2
    Threads.@threads for j0 in 0:m-1
        @maybe_aliasscope begin
            j = firstindex(xs, 2) + j0 + 1
            for i0 in 0:n-1
                i = firstindex(xs, 1) + i0 + 1
                @inbounds ys[i, j] = avg3x3(xs, i, j)
            end
        end
    end
    return ys
end

function avgfilter2d_tapir_dac!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = size(xs)[1] - 2
    m = size(xs)[2] - 2
    @maybe_aliasscope begin
        Tapir.@par dac for j0 in 0:m-1
            j = firstindex(xs, 2) + j0 + 1
            for i0 in 0:n-1
                i = firstindex(xs, 1) + i0 + 1
                @inbounds ys[i, j] = avg3x3(xs, i, j)
            end
            @grainsize 64
        end
    end
    return ys
end

function avgfilter2d_tapir_seq!(ys, xs)
    @assert axes(ys) == axes(xs)
    n = size(xs)[1] - 2
    m = size(xs)[2] - 2
    @maybe_aliasscope begin
        Tapir.@par seq for j0 in 0:m-1
            j = firstindex(xs, 2) + j0 + 1
            for i0 in 0:n-1
                i = firstindex(xs, 1) + i0 + 1
                @inbounds ys[i, j] = avg3x3(xs, i, j)
            end
            @grainsize 64
        end
    end
    return ys
end


text:
"""
    integerfeasibility(A,d,nullVecs=false)

Given a linear system `Ax=d`, return an integer vector `x` which satisfies the
system.  This is the *integer programming feasibility problem*.

If `nullVecs==true`, then as well as returning a solution `x`, also return a
matrix `xNull` of vectors in the null space of `A` which could be added to the
`x` vector to find a solution which satisfies a constraint such as `0 .≤ x
.≤ u`; see the paper below.

This is not a robust tool, just a demo.

"Solving A System Of Diophantine Equations With Bounds On The Variables" by
Karen Aardal, Cor Hurkens, and Arjen Lenstra in Integer Programming and
Combinatorial Optimization, 6th International IPCO Conference, vol 1412, pp
229-242, 1998. See http://softlib.rice.edu/pub/CRPC-TRs/reports/CRPC-TR98782.pdf

# Examples
```jldoctest
julia> A=[10 1 -9; 1 8 8]; xtrue=[0; 2; 9]; d=A*xtrue;

julia> integerfeasibility(A,d)
3-element Array{Int64,1}:
 0
 2
 9

julia> n=20;m=30; A = rand(-10:10,n,m); xtrue = rand(0:10,m); d=A*xtrue;

julia> sum(abs.(xtrue - integerfeasibility(A,d) ))
0

```
"""
function integerfeasibility(A::AbstractArray{Td,2},
                            d::AbstractArray{Td,1}, flag=false) where {Td<:Number}
    m,n=size(A)
    N1=10^3
    N2=10^4
    B = [[I           zeros(Td,n,1)];
         zeros(Td,1,n) N1;
         N2*A          -N2*d]

    Bhat,_ = lll(B)
    Bprime = Bhat[1:n+1,1:n-m+1]
    x_d = Bprime[1:n,end]

    if flag
        xNull=Bprime[1:n,1:n-m]
        return x_d,xNull
    else
        return x_d
    end
end

"""
    x = subsetsum(a,s)

For a vector of integers `a`, and an integer `s`, try to find a binary
vector `x` such that `x'*a=s`. We use the LLL algorithm to find the
solution. This is not a robust tool, just a demo.

This follows the technique described by Lagarias and Odlyzko  in
"Solving Low-Density Subset Sum Problems"  in Journal of ACM, Jan 1985.
Code based on http://web.eecs.umich.edu/~cpeikert/lic15/lec05.pdf
We can likely get better results using techniques described and referenced in
https://www-almasty.lip6.fr/~joux/pages/papers/ToolBox.pdf

It's odd that permuting the `a` vector in the second example given below
causes the alg to often not find a solution. The example doesn't have the
required density, so maybe we should be asking why it finds an answer at
all.

# Examples
```jldoctest
julia> a=[1,2,4,9,20,38]; s=30; x=subsetsum(a,s); s-x'*a
0.0

julia> a=[32771,65543,131101,262187,524387,1048759, # from Bremner p 117
          2097523,4195057,8390143,16780259,33560539,
          67121039,134242091,268484171,536968403];

julia> s=891221976; x=subsetsum(a,s); s-x'*a
0.0

julia> N=40;a=rand(1:2^BigInt(256),N);xtrue=rand(Bool,N); s=a'*xtrue;

julia> setprecision(BigFloat,300); x=subsetsum(a,s); s-x'*a
0.0

```
"""
function subsetsum(a::AbstractArray{Ti,1},ss::Ti,returnBinary=false) where {Ti<:Integer}
    # page numbers below refer to lecture note above

    n = length(a)

    if ss<sum(a)/2
        flag=true
        s = sum(a)-ss
    else
        flag = false
        s = ss
    end
    # b is the "B" parameter defined at bottom of page 2
    nt = Ti(n)
    b = ceil(sqrt(nt*2^nt))
    # BB is the 'B' matrix defined also at bottom of page 2
    BB = [Matrix{Ti}(I,n,n+1); -b*a' b*s]
    B,T = lll(BB)

    # Pick off the interesting parts of the B matrix
    Bp = B[1:n,1:n+1]
    # Check if any of the columns satisfy the requirement
    sums = abs.(a'*Bp)
    ixMatch = findall(sums.==s)
    # d = diag(Bp'*Bp)  # solution generally matches a small element of d

    binarySolution = false
    global xb
    for ix=1:length(ixMatch)
        # Pick off a column of B
        x = Bp[:,ixMatch[ix][2]]
        if !(any(x.<0)|any(x.>1))
            binarySolution= true
            xb = flag ? mod.(x .+ 1,2) : x
        end
        x = -x
        if !(any(x.<0)|any(x.>1))
            binarySolution= true
            xb = flag ? mod.(x .+ 1,2) : x
        end
    end

    if binarySolution
        return xb
    else
        x = mdsubsetsum(a,ss,.5,1)
        if ~ismissing(x)
            print("A solution was found via mdsubsetsum\n")
            return x
        end

        if maximum(a)<2^(n^2/2)
            density = n/maximum(log2.(a))
            @printf("The density (%4.2f) of the 'a' vector is not as low as required\n",
                    density)
            @printf("(%4.2f) for Lagarias-Odlyzko to work. \n",2/n)
        end

        if returnBinary && length(ixMatch)>=1
            print("A non-binary solution was found; check that it's correct\n")
            return Bp[:,ixMatch[1][2]]
        else
            print("A solution was not found\n")
        end
    end
    return missing.*Bp[:,1]
end

"""
    x = mdsubsetsum(a,sM,ratio=.5,Kpm=3)

For a vector of integers `a`, and an integer `sM`, try to find a binary
vector `x` such that `x'*a=s` using the technique from "Multidimensional
subset sum problem" [1][2]. A major goal of the technique is to solve
problems in there are about 50% ones in `x`; other ratios of ones to zeros
can be specified in `ratio`.  The thesis also suggests searching `Kpm=3`
values around the nominal k. This technique is related to that in
[`subsetsum`](@ref) in that both use the LLL algorithm.  This is not a
robust tool, just a demo.

[1] https://scholarworks.rit.edu/theses/64/
[2] https://pdfs.semanticscholar.org/21a7/c2f9ff29507f1153aefcca04d1cd308e45c0.pdf

# Examples
```jldoctest
julia> a=[1,2,4,9,20,38]; s=30; x=mdsubsetsum(a,s); s-x'*a
0.0

julia> a=[32771,65543,131101,262187,524387,1048759, # from Bremner p 117
          2097523,4195057,8390143,16780259,33560539,
          67121039,134242091,268484171,536968403];

julia> sM=891221976; x=mdsubsetsum(a,sM); sM-x'*a

julia> N=40;a=rand(1:2^BigInt(256),N);xtrue=rand(Bool,N); s=a'*xtrue;

julia> setprecision(BigFloat,300); x=mdsubsetsum(a,s); s-x'*a
```
"""
function mdsubsetsum(a::AbstractArray{Ti,1},sM::Ti,ratio=.5,Kpm=3) where {Ti<:Integer}

    n = Ti(length(a))

    # r is heuristic value for r from end of Section 11 of thesis above
    r = Ti(2^round(.65*log2(maximum(a))))
    c = 2^10 # a heuristic; look right below eq (10)

    p = a .÷ r # See start of Section 10
    s = a-p.*r
    k0 = Ti(round(sum(s)*ratio/r)) # Equation (9) with ratio=.5 and j=1.
    p0 = sM ÷ r
    m = sM - p0*r

    for k = k0-Kpm:k0+Kpm
        # Equation (14)
        Bt = [Matrix{Ti}(I,n,n)*2 c*s       c*p      zeros(Ti,n);
              ones(Ti,1,n)        c*(k*r+m) c*(p0-k) 1]
        BB = Bt'  # Notation in paper is row-based; my brain is column-based
        B,T = lll(BB)

        for nx = 1:n
            if abs(B[n+3,nx])==1 && B[n+2,nx]==0 && B[n+1,nx]==0 &&
                      all((B[1:n,nx] .==1) .| (B[1:n,nx] .==-1))
                xhat = abs.(B[1:n,nx] .- B[n+3,nx])/2
                return Ti.(xhat)
            end
        end
    end
    @warn "Solution not found"
    return missing
end

"""
    rationalapprox(x::AbstractArray{<:Real,1},M,Ti=BigInt,verbose=false)

For a vector of Reals `x`, and an integer `M`, find an integer q such that
`maximum(abs.(x*q-round.(x*q)))` is small; the vector `x` is approximated by
`round.(x*q)//q`.  The integer `q` is less than or equal to `M` and the
approximation satisfies `max(abs.(x*q-round.(x*q)))≤sqrt(5)*2^(n/4 -
5)*M^(-1/n)`; this equation comes from the paper below.  The LLL algorithm
reduction is used to find the solution. The approximation vector is returned.
This is also known as "simultaneous diophantine approximation"; see for
example the title of the Hanrot paper below.

This is not a robust tool, just a demo.

"LLL: A Tool for Effective Diophantine Approximation" by Guillaume Hanrot in
the book "The LLL Algorithm: Survey and Applications" edited by Phong
Q. Nguyen and Brigitte Vallée, Springer, Heidelberg, 2010.

See also Chapter 9 of M. R. Bremner, "Lattice Basis Reduction: An
Introduction to the LLL Algorithm and Its Applications" CRC Press, 2012.

# Examples
```jldoctest
julia> x = [0.3912641745333527; 0.5455179974014548; 0.1908698210882469];
julia> rationalapprox(x,1e4,Int64)
3-element Array{Rational{Int64},1}:
 43//110
  6//11
 21//110
```
"""
function rationalapprox(x::AbstractArray{<:Real,1},M,Ti=BigInt,verbose=false)
    # Follows the multivariate technique in Theorem 8 of the
    # Nguyen and Vallee paper above

    n = length(x)
    Q = (M/2^(1/4))^(1/n)
    C = Ti(ceil(Q^(n+1)))
    Cx = Ti.(round.(C*x))

    X = [one(Ti) zeros(Ti,1,n);
         [Cx     C*I] ]
    B,_ = lll(X)
    q = abs(B[1,1])
    maxErr = maximum(abs.(x-round.(x*q)/q))
    bound = sqrt(5)*2^(n/4 - 1)/Q /q
    if maxErr>bound
        warning("The max error in the approximation $(maxErr) is unexpectedly "*
                "greater than the bound $(bound).")
    end
    p = Ti.(round.(x*q));
    if verbose
        println("q        = $(q)\n"*
                "maxErr/q = $(maxErr)\n"*
                "bound/q  = $(bound)")
        display([x round.(x*q)/q abs.(x-round.(x*q)/q)])
    end
    return p//q
end


text:
# TODO: this should end up in BioSequences.jl?

"Extract the element stored in a packed bitarray referred to by bidx."
@inline function BioSequences.extract_encoded_element(bidx::BioSequences.BitIndex{N,W}, data::NTuple{n,W}) where {N,n,W}
    @inbounds chunk = data[BioSequences.index(bidx)]
    offchunk = chunk >> (BioSequences.bitwidth(W) - N - BioSequences.offset(bidx))
    return offchunk & BioSequences.bitmask(bidx)
end

"""
    _cliphead(by::Integer, head::UInt64, tail...)

A method used to mask the first `by` MSB's in `head`, before catting it with
tail to return a NTuple.

This is used internally to mask the first `by` bits in the first word of a
NTuple of UInt64's.

Notably it's used when constructing a Kmer from an existing NTuple of UInt64
"""
@inline function _cliphead(by::Integer, head::UInt64, tail...)
    return (head & (typemax(UInt64) >> by), tail...)
end

#=
rightshift_carry & leftshift_carry

These methods are micro-optimised (or should be!!!) for shifting the bits in
an NTuple of unsigned integers, carrying the bits "shifted off" one word
over to the next word. The carry can also be "seeded" so as other methods like
pushfirst and pushlast can be efficiently implemented without duplication of code
or less efficient implementations that first shift and then insert an element.
=#

@inline function rightshift_carry(x::NTuple{N,UInt64}, nbits::Integer, prevcarry = zero(UInt64)) where {N}
    return _rightshift_carry(nbits, prevcarry, x...)
end

@inline function _rightshift_carry(nbits::Integer, carry::UInt64, head::UInt64, tail...)
    return ((head >> nbits) | carry, _rightshift_carry(nbits, (head & ((one(UInt64) << nbits) - 1)) << (64 - nbits), tail...)...)
end

@inline _rightshift_carry(nbits::Integer, carry::UInt64) = ()

@inline function leftshift_carry(x::NTuple{N,UInt64}, nbits::Integer, prevcarry::UInt64 = zero(UInt64)) where {N}
    _, newbits = _leftshift_carry(nbits, prevcarry, x...)
    return newbits
end

@inline function _leftshift_carry(nbits::Integer, prevcarry::UInt64, head::UInt64, tail...)
    carry, newtail = _leftshift_carry(nbits, prevcarry, tail...)
    return head >> (64 - nbits), ((head << nbits) | carry, newtail...)
end

@inline _leftshift_carry(nbits::Integer, prevcarry::UInt64) = prevcarry, ()

@inline function _reverse(bpe::BioSequences.BitsPerSymbol{N}, head::UInt64, tail...) where {N}
    return (_reverse(bpe, tail...)..., BioSequences.reversebits(head, bpe))
end

@inline _reverse(::BioSequences.BitsPerSymbol{N}) where {N} = ()

#=
@inline function _reverse(f::F, bpe::BioSequences.BitsPerSymbol{N}, head::UInt64, tail...) where {N,F<:Function}
    return (_reverse(f, bpe, tail...)..., f(reversebits(head, bpe)))
end

@inline _reverse(f::F, ::BioSequences.BitsPerSymbol{N}) where {N,F<:Function} = ()
=#


text:
# Transform
struct Translate <: Group
    vector::AbstractArray{Real, 1}
    elements::AbstractArray{Op, 1}
end

function toString(t::Translate)
    output = "translate([$(t.vector[1]), $(t.vector[2]), $(t.vector[3])]) {"
    for c in t.elements
        output = output * toString(c)
    end
    output * "}"
end


text:
# Don't move this line from the top, please. {{X}} {{Y}} {{Z}}

@info "Running plugin tests"

struct FileTest <: PT.FilePlugin
    a::String
    b::Bool
end

PT.gitignore(::FileTest) = ["a", "aa", "aaa"]
PT.source(::FileTest) = @__FILE__
PT.destination(::FileTest) = "foo.txt"
PT.badges(::FileTest) = PT.Badge("{{X}}", "{{Y}}", "{{Z}}")
PT.view(::FileTest, ::Template, ::AbstractString) = Dict("X" => 0, "Y" => 2)
PT.user_view(::FileTest, ::Template, ::AbstractString) = Dict("X" => 1, "Z" => 3)

@testset "Plugins" begin
    @testset "FilePlugin" begin
        p = FileTest("foo", true)
        t = tpl(; plugins=[p])

        # The X from user_view should override the X from view.
        s = PT.render_plugin(p, t, "")
        @test occursin("1 2 3", first(split(s, "\n")))

        with_pkg(t) do pkg
            pkg_dir = joinpath(t.dir, pkg)
            badge = string(PT.Badge("1", "2", "3"))
            @test occursin("a\naa\naaa", read(joinpath(pkg_dir, ".gitignore"), String))
            @test occursin(badge, read(joinpath(pkg_dir, "README.md"), String))
            @test read(joinpath(pkg_dir, "foo.txt"), String) == s
        end
    end

    @testset "Tests Project.toml warning on Julia < 1.2" begin
        p = Tests(; project=true)
        @test_logs (:warn, r"The project option is set") tpl(; julia=v"1", plugins=[p])
        @test_logs (:warn, r"The project option is set") tpl(; julia=v"1.1", plugins=[p])
        @test_logs tpl(; julia=v"1.2", plugins=[p])
        @test_logs tpl(; julia=v"1.3", plugins=[p])
    end

    @testset "Equality" begin
        a = FileTest("foo", true)
        b = FileTest("foo", true)
        @test a == b
        c = FileTest("foo", false)
        @test a != c
    end
end


text:
export ImmutablePolynomial

"""
    ImmutablePolynomial{T<:Number}(coeffs::AbstractVector{T}, var=:x)

Construct an immutable (static) polynomial from its coefficients `a`,
lowest order first, optionally in terms of the given variable `x`
where `x` can be a character, symbol, or string.

If ``p = a_n x^n + \\ldots + a_2 x^2 + a_1 x + a_0``, we construct
this through `ImmutablePolynomial((a_0, a_1, ..., a_n))` (assuming
`a_n ≠ 0`). As well, a vector or number can be used for construction.

The usual arithmetic operators are overloaded to work with polynomials
as well as with combinations of polynomials and scalars. However,
operations involving two polynomials of different variables causes an
error, though for `+` and `*` operations, constant polynomials are
treated as having no variable.

As the coefficient size is a compile-time constant, immutable
polynomials can take advantage of faster polynomial evaluation
provided by `evalpoly` from Julia 1.4.

    # Examples

```jldoctest
julia> using  Polynomials

julia> ImmutablePolynomial((1, 0, 3, 4))
ImmutablePolynomial(1 + 3*x^2 + 4*x^3)

julia> ImmutablePolynomial((1, 2, 3), :s)
ImmutablePolynomial(1 + 2*s + 3*s^2)

julia> one(ImmutablePolynomial)
ImmutablePolynomial(1.0)
```

!!! note
    This was modeled after https://github.com/tkoolen/StaticUnivariatePolynomials.jl by @tkoolen.

"""
struct ImmutablePolynomial{T <: Number,  N} <: StandardBasisPolynomial{T}
    coeffs::NTuple{N, T}
    var::Symbol
    function ImmutablePolynomial{T,N}(coeffs::NTuple{N,T}, var::SymbolLike=:x) where {T <: Number,N}
        N == 0 && return new{T,0}(coeffs, Symbol(var))
        iszero(coeffs[end]) &&  throw(ArgumentError("Leading  term must  be  non-zero"))
        new{T,N}(coeffs, Symbol(var))
    end
    function ImmutablePolynomial{T,N}(coeffs::AbstractVector{S}, var::SymbolLike=:x) where {T <: Number, N, S}
        M = findlast(!iszero, coeffs)
        M == N || throw(ArgumentError("Leading  term must  be  non-zero; length of coefficients must be N"))
        cs = NTuple{N,T}(T(c) for c in coeffs[1:N])
        new{T,N}(cs, var)
    end
    function ImmutablePolynomial{T}(coeffs::NTuple{M,S}, var::SymbolLike=:x) where {T, S<: Number, M}
        N = findlast(!iszero, coeffs)
        if N == nothing
            return zero(ImmutablePolynomial{T}, var)
        else
            cs = NTuple{N,T}(coeffs[i] for  i in  1:N)
        end
        new{T,N}(cs, var)
    end
    function ImmutablePolynomial{T}(coeffs::Vector{M}, var::SymbolLike=:x) where {T, M}
        N = findlast(!iszero, coeffs)
        if N == nothing
            return zero(ImmutablePolynomial{T}, var)
        else
            cs = NTuple{N,T}(coeffs[i] for  i in  1:N)
        end
        new{T,N}(cs, var)
    end

end

@register ImmutablePolynomial

function ImmutablePolynomial{T,N}(coeffs::Tuple, var)  where {T,N}
    ImmutablePolynomial{T,N}(NTuple{N,T}(c for c in coeffs), var)
end

function ImmutablePolynomial(coeffs::NTuple{M,T}, var::SymbolLike=:x) where {M, T}
    N = findlast(!iszero, coeffs)
    if N == nothing
        return zero(ImmutablePolynomial, var) # no  type  passed in NTuple{M,T}
    else
        cs = NTuple{N,T}(T(c) for  c in  coeffs[1:N])
        end
    ImmutablePolynomial{T,N}(cs, var)
end

# Convenience; pass tuple to Polynomial
# Not documented, not sure this is a good idea as P(...)::P is not true...
Polynomial(coeffs::NTuple{N,T}, var::SymbolLike = :x) where{N,T} =
    ImmutablePolynomial(coeffs, var)
function Polynomial{T}(coeffs::NTuple{N,S}, var::SymbolLike = :x) where{N,T,S}
    ImmutablePolynomial{N,T}(T.(coeffs), var)
end

##
## ----
##
# overrides from common.jl due to coeffs possibly being padded, coeffs being non mutable, ...

Base.collect(p::P) where {P <: ImmutablePolynomial} = [pᵢ for pᵢ ∈ p]

Base.copy(p::P) where {P <: ImmutablePolynomial} = P(coeffs(p), p.var)

# catch q == 0 case
LinearAlgebra.norm(q::ImmutablePolynomial{T}, p::Real = 2) where {T} = degree(q) == -1 ? zero(T) : norm(coeffs(q), p)

#  zero, one, variable
function Base.zero(P::Type{<:ImmutablePolynomial},var::SymbolLike=:x)
    R = eltype(P)
    ImmutablePolynomial{R,0}(NTuple{0,R}(),var)
end

function  Base.one(P::Type{<:ImmutablePolynomial},var::SymbolLike=:x)
    R = eltype(P)
    ImmutablePolynomial{R,1}(NTuple{1,R}(1),var)
end
function variable(P::Type{<:ImmutablePolynomial},var::SymbolLike=:x)
    R  = eltype(P)
    ImmutablePolynomial{R,2}(NTuple{2,R}((0,1)),var)
end

# degree, isconstant
degree(p::ImmutablePolynomial{T,N}) where {T,N} = N - 1
isconstant(p::ImmutablePolynomial{T,N}) where {T,N}  = N <= 1

function Base.getindex(p::ImmutablePolynomial{T,N}, idx::Int) where {T <: Number,N}
    (idx <  0 || idx > N-1) && return zero(T)
    return p.coeffs[idx + 1]
end

Base.setindex!(p::ImmutablePolynomial, val::Number,  idx::Int) = throw(ArgumentError("Immutable polynomials are immutable"))

for op in [:isequal, :(==)]
    @eval function Base.$op(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
        check_same_variable(p1,p2) || return false
        p1s, p2s = coeffs(p1), coeffs(p2)
        (N == M  && $op(p1s,p2s)) &&  return  true
        n1 = findlast(!iszero, p1s) # now trim out zeros
        n2 = findlast(!iszero, p2s)
        (n1 == nothing && n2 == nothing) && return true
        (n1 == nothing || n2  == nothing) && return false
        $op(p1s[1:n1],p2s[1:n2]) &&  return true
        false
    end
end

##
## Change ≈ to handle tuples for coefficients *and* get handling of Inf correct
## right now
## Inf ≈ Inf # true
## [Inf] ≈ [Inf] # true
## P([Inf]) ≈ P([Inf]) # false
## P([Inf]) ≈ Inf # false
## This fixes the last two cases for P=ImmutablePolynomial, and could replace
## isapprox in common.jl
##
## check Inf values (matching Vector ≈ Vector)
## @test P([Inf]) ≈ Inf
## @test !(P([Inf]) ≈ P([-Inf])) # default compares zero(P) to zero(P)
## @test !(P([1,2,Inf]) ≈ P([1,3,Inf])) # default uses truncate with  Inf which clobbers numbers
function Base.isapprox(p1::ImmutablePolynomial{T},
    p2::ImmutablePolynomial{S};
    rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,S}

    p1, p2 = promote(p1, p2)
    check_same_variable(p1, p2)  || error("p1 and p2 must have same var")

    # copy over from abstractarray.jl
    Δ  = norm(p1-p2)
    if isfinite(Δ)
        return Δ <= max(atol, rtol*max(norm(p1), norm(p2)))
    else
        for i in 0:max(degree(p1), degree(p2))
            isapprox(p1[i], p2[i]; rtol=rtol, atol=atol) || return false
        end
        return true
    end
end

function Base.isapprox(p1::P,
                       n::S;
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,S, P<:ImmutablePolynomial{T}}
    return isapprox(p1, ⟒(P){T}(n,p1.var))
end

function Base.chop(p::ImmutablePolynomial{T,N};
              rtol::Real = Base.rtoldefault(real(T)),
              atol::Real = 0)  where {T,N}
    cs = coeffs(p)
    for i in N:-1:1
        if !isapprox(cs[i], zero(T), rtol=rtol, atol=atol)
            return ImmutablePolynomial{T,i}(cs[1:i], p.var)
        end
    end
    zero(ImmutablePolynomial{T})
end

function Base.truncate(p::ImmutablePolynomial{T,N};
                       rtol::Real = Base.rtoldefault(real(T)),
                       atol::Real = 0)  where {T,N}
    q = chop(p, rtol=rtol, atol=atol)
    iszero(q) && return  q
    cs = coeffs(q)
    thresh = maximum(abs,cs) * rtol + atol
    ImmutablePolynomial(map(c->abs(c) <= thresh ? zero(T) : c, coeffs(q)), p.var)
end

# no in-place chop! and truncate!
chop!(p::ImmutablePolynomial; kwargs...) =  chop(p; kwargs...)
truncate!(p::ImmutablePolynomial; kwargs...) =  truncate(p; kwargs...)

LinearAlgebra.conj(p::P) where {P <: ImmutablePolynomial} = P(conj([aᵢ for aᵢ in coeffs(p)]))

##
## --------------------
##

(p::ImmutablePolynomial{T,N})(x::S) where {T,N,S} = evalpoly(x, coeffs(p))

function Base.:+(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}

    R = promote_type(S,T)
    iszero(N) && return ImmutablePolynomial{R}(coeffs(p2), p2.var)
    iszero(M) && return ImmutablePolynomial{R}(coeffs(p1), p1.var)

    isconstant(p1) && p1.var != p2.var && return p2 + p1[0]*one(ImmutablePolynomial{T},p2.var)
    isconstant(p2) && p1.var != p2.var && return p1 + p2[0]*one(ImmutablePolynomial{S},p1.var)

    p1.var != p2.var && error("Polynomials must have same variable")

    if  N == M
        cs = coeffs(p1) .+ coeffs(p2)
        return ImmutablePolynomial(cs, p1.var) # N  unknown, as leading terms can cancel
    elseif N < M
        ⊕(p2,p1)
    else
        ⊕(p1,p2)
    end

end

# assume N > M.
# As N ≠ M, we are assured of size of output (max(N,M)), so we generate the function
@generated function ⊕(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
    R = promote_type(T,S)

    exprs = Any[nothing for i = 1:N]
    for i in  1:M
        exprs[i] = :(p1.coeffs[$i] + p2.coeffs[$i])
    end
    for i in M+1:N
        exprs[i] =:(p1.coeffs[$i])
    end

    return quote
        Base.@_inline_meta
        ImmutablePolynomial{$(R),$(N)}(tuple($(exprs...)), p1.var)
    end

end

function Base.:*(p1::ImmutablePolynomial, p2::ImmutablePolynomial)
    isconstant(p1) && return p2 * p1[0]
    isconstant(p2) && return p1 * p2[0]
    p1.var != p2.var && error("Polynomials must have same variable")
    p1 ⊗ p2
end

## Static size of product makes generated functions  a good choice
## from https://github.com/tkoolen/StaticUnivariatePolynomials.jl/blob/master/src/monomial_basis.jl
@generated function ⊗(p1::ImmutablePolynomial{T,N}, p2::ImmutablePolynomial{S,M}) where {T,N,S,M}
    P = M + N - 1
    R = promote_type(T,S)
    exprs = Any[nothing for i = 1 : P]
    for i in 1 : N
        for j in 1 : M
            k = i + j - 1
            if exprs[k] === nothing
                exprs[k] = :(p1.coeffs[$i] * p2.coeffs[$j])
            else
                exprs[k] = :(muladd(p1.coeffs[$i], p2.coeffs[$j], $(exprs[k])))
            end
        end
    end

    return quote
        Base.@_inline_meta
        ImmutablePolynomial{$(R),$(max(0,P))}(tuple($(exprs...)), p1.var)
    end

end

# scalar ops
function Base.:+(p::ImmutablePolynomial{T,N}, c::S) where {T, N, S<:Number}
    R = promote_type(T,S)
    iszero(c) && return ImmutablePolynomial{R,N}(R.(p.coeffs), p.var)
    N == 0 && return ImmutablePolynomial{R,1}((R(c),), p.var)
    N == 1 && return ImmutablePolynomial(R[p[0]+c])

    return p + ImmutablePolynomial{R,1}(R[c],p.var)

    cs = NTuple{N,R}(iszero(i) ? p[i]+c : p[i] for i in 0:N-1)
    return ImmutablePolynomial{R,N}(cs, p.var)
end

function Base.:*(p::ImmutablePolynomial{T,N}, c::S) where {T, N, S <: Number}
    R = promote_type(T,S)
    iszero(c)  && return zero(ImmutablePolynomial{R})
    ImmutablePolynomial{R,N}(p.coeffs .* c, p.var)
end

function Base.:/(p::ImmutablePolynomial{T,N}, c::S) where {T,N,S <: Number}
    R = eltype(one(T)/one(S))
    isinf(c)  && return zero(ImmutablePolynomial{R})
    ImmutablePolynomial{R,N}(p.coeffs ./ c, p.var)
end

Base.:-(p::ImmutablePolynomial{T,N}) where {T,N} = ImmutablePolynomial(NTuple{N,T}(-pi for pi in p.coeffs), p.var)

Base.to_power_type(p::ImmutablePolynomial{T,N}) where {T,N} = p


text:
macro implement_array_methods(t)
    quote
        Base.Array(a::$t) = DiskArrays._Array(a)
        Base.copyto!(dest::$t, source::AbstractArray) = DiskArrays._copyto!(dest, source)
        Base.copyto!(dest::AbstractArray, source::$t) = DiskArrays._copyto!(dest, source)
        Base.copyto!(dest::$t, source::$t) = DiskArrays._copyto!(dest, source)
        function Base.copyto!(dest::$t, Rdest::CartesianIndices, src::AbstractArray, Rsrc::CartesianIndices)
            DiskArrays._copyto!(dest, Rdest, src, Rsrc)
        end
        function Base.copyto!(dest::AbstractArray, Rdest::CartesianIndices, src::$t, Rsrc::CartesianIndices)
            DiskArrays._copyto!(dest, Rdest, src, Rsrc)
        end
        function Base.copyto!(dest::$t, Rdest::CartesianIndices, src::$t, Rsrc::CartesianIndices)
            DiskArrays._copyto!(dest, Rdest, src, Rsrc)
        end
        Base.reverse(a::$t, dims=:) = DiskArrays._reverse(a, dims)
        # Here we extend the unexported `_replace` method, but we replicate
        # much less Base functionality by extending it rather than `replace`.
        function Base._replace!(new::Base.Callable, res::AbstractArray, A::$t, count::Int)
            DiskArrays._replace!(new, res, A, count)
        end
        function Base._replace!(new::Base.Callable, res::$t, A::AbstractArray, count::Int)
            DiskArrays._replace!(new, res, A, count)
        end
        function Base._replace!(new::Base.Callable, res::$t, A::$t, count::Int)
            DiskArrays._replace!(new, res, A, count)
        end
    end
end

# Use broadcast to copy to a new Array
function _Array(a::AbstractArray{T,N}) where {T,N}
    dest = Array{T,N}(undef, size(a))
    dest .= a
    return dest
end

# Use broadcast to copy
_copyto!(dest::AbstractArray{<:Any,N}, source::AbstractArray{<:Any,N}) where N = dest .= source
function _copyto!(dest::AbstractArray, source::AbstractArray)
    # TODO make this more specific so we are reshaping the Non-DiskArray more often.
    reshape(dest, size(source)) .= source
    return dest
end
_copyto!(dest, Rdest, src, Rsrc) = view(dest, Rdest) .= view(src, Rsrc)

# Use a view for lazy reverse
_reverse(a, dims::Colon) = _reverse(a, ntuple(identity, ndims(a)))
_reverse(a, dims::Int) = _reverse(a, (dims,))
function _reverse(A, dims::Tuple)
    rev_axes = map(ntuple(identity, ndims(A)), axes(A)) do d, a
        ax = StepRange(a)
        d in dims ? reverse(ax) : ax
    end
    return view(A, rev_axes...)
end

# Use broadcast instead of a loop.
# The `count` argument is disallowed as broadcast is not sequential.
function _replace!(new, res::AbstractArray, A::AbstractArray, count::Int)
    count < length(res) && throw(ArgumentError("`replace` on DiskArrays objects cannot use a count value"))
    broadcast!(new, res, A)
end


text:
# Benchmark the parsing of a PDB file given as an argument

using MIToS.PDB

pdb_filepath = ARGS[1]

# Run to JIT compile
read(pdb_filepath, PDBFile)

elapsed = @elapsed struc = read(pdb_filepath, PDBFile)

println(elapsed)

max_stars_repo_name: wojdyr/pdb-benchmarks

max_issues_repo_name: wojdyr/pdb-benchmarks

max_forks_repo_name: wojdyr/pdb-benchmarks


text:
"""
```
Kalman Smoother(Model, Output, Output)
```

Takes as input the output from the Kalman Filter, i.e. the predicted
and filtered of the means and covariances of the states and measurements.

The function returns the type Output which smoothed means and covariances
for the states and measurements, respectively.
"""

function kalmansmoother(model::Model, predicted::Predicted, filtered::Filtered)
    T, p = size(model.y)
    m = length(model.x)
    ks = Smoother(T, p, m)
    ks.x[T, :]    = filtered.x[T, :]
    ks.P[T, :, :] = filtered.P[T, :, :]
    ks.μ[T, :]    = update_μ(ks.x[T, :], model.H)
    ks.Σ[T, :, :] = update_Σ(ks.P[T, :, :], model.H, model.R)
    for t in T-1:-1:1
        J             = filtered.P[t, :, :] * model.A' * inv(predicted.P[t + 1, :, :])
        ks.x[t, :]    = filtered.x[t, :] + J * (ks.x[t + 1, :] - predicted.x[t + 1, :])
        ks.P[t, :, :] = filtered.P[t, :, :] + J * (ks.P[t + 1, :, :] - predicted.P[t + 1, :, :]) * J'
        ks.μ[t, :]    = update_μ(ks.x[t, :], model.H)
        ks.Σ[t, :, :] = update_Σ(ks.P[t, :, :], model.H, model.R)
    end
    return ks
end


text:
"""
This file shows the code for S-S model.
The S-S model is based on Sobel and Schneider 2009 and 2013 papers, please refer to the papers for more details.
"""

## ----------------------------------------------
# model configuration
const integration_year_number = 10          # integrate how many years
const resolution_factor = 1                 # resolution factor. If less than 1, use reduced resolution to save time.
const dt=60*floor(1/resolution_factor+0.5)  # satisfies CFL requirement if wind speed is not insanely large
                                            # unit: s
Dir="D:/SS_model_output/YearCycle_ITCZDisplaced_MomentumInputCycle/LowRes_Kv100/" # path to save output files
const ITCZ_displaced = true                 # If true, annual mean ITCZ is located at 6 deg N.
                                            # ITCZ is controlled by radiative forcing
# ----------------------------------------------

## ----------------------------------------------
# Define constants
# For the meaning of constants, refer to Sobel and Schneider 2009
const β=2e-11
const kᵥ=7786*100
const ϵᵤ=1e-8
const Δz=60
const Δy=50
const T₀=300
const δ=4e3
const g=9.81
const H=16e3
const τ=37*24*3600
const y₁=9439e3
# ----------------------------------------------

# ----------------------------------------------
# define dimensions
if iseven(Int(floor(801*resolution_factor)))
    const y_Grids_for_V = Int(floor(801*resolution_factor)+1)
else
    const y_Grids_for_V = Int(floor(801*resolution_factor))
end
const y_Grids_for_UandT = y_Grids_for_V-1
const dy = 15751e3*2/y_Grids_for_UandT

const y_for_V=range(-15751e3, stop=15751e3,length=y_Grids_for_V)
const y_for_UandT=range(-15751e3+0.5*dy, step=dy ,length=y_Grids_for_UandT)

time = [0]
timestamp=0  # which day in the simulation
if ITCZ_displaced==ture
    solar_insolation_degree=23.5*sin(2*pi*(timestamp%365-79)/365)+6
else
    solar_insolation_degree=23.5*sin(2*pi*(timestamp%365-79)/365)
end
# ----------------------------------------------

# ----------------------------------------------
# Define variables u, v and θ. (time, y)
u=fill(0.0,1,y_Grids_for_UandT)
v=fill(0.0,1,y_Grids_for_V)
# ----------------------------------------------

# ----------------------------------------------
# initialize theta_e and u.
# Define the radiative-convective equilibrium (RCE) temperature θₑ
θₑ=zeros(y_Grids_for_UandT)
θ00=330
u_initial=2*g*H*Δy/(1.6*T₀*β*y₁^2)
S=zeros(y_Grids_for_UandT)
u[1,:].=u_initial
for i=1:y_Grids_for_UandT
    if abs(y_for_UandT[i]-111e3*solar_insolation_degree) <= y₁
        θₑ[i]=θ00-Δy*((y_for_UandT[i]-111e3*solar_insolation_degree)/y₁).^2
        # u[:,i].=u_initial
    else
        θₑ[i]=θ00-Δy
    end
end
θ=reshape(θₑ,1,y_Grids_for_UandT)

function get_θₑ(timestamp)
    θₑ=zeros(y_Grids_for_UandT)
    θ00=330
    if ITCZ_displaced==ture
        solar_insolation_degree=23.5*sin(2*pi*(timestamp%365-79)/365)+6
    else
        solar_insolation_degree=23.5*sin(2*pi*(timestamp%365-79)/365)
    end
    for i=1:y_Grids_for_UandT
        if abs(y_for_UandT[i]-111e3*solar_insolation_degree) <= y₁
            θₑ[i]=θ00-Δy*((y_for_UandT[i]-111e3*solar_insolation_degree)/y₁).^2
        else
            θₑ[i]=θ00-Δy
        end
    end
    return θₑ
end
# ----------------------------------------------

# ----------------------------------------------
# Define necessary Functions
# Heaviside function
heaviside(x::Float64) = ifelse(x < 0, zero(x), ifelse(x > 0, one(x), oftype(x,0.5)))

# Get potential temperature
# According to the paper, the factor (pₛ/pₜ)^(R/cₚ) is 1.6
function get_θ(T)
    return 1.6*T
end

function get_T(θ)
    return θ/1.6
end

# The eddy momentum flux divergence function is based on the momentum budget in MERRA-2 reanalysis data.
# Please refer to Zhang and Lutsko 2022 for details.
function get_S(timestamp,i)
    day_of_year=timestamp%365
    if day_of_year<60 || day_of_year>285
        S=-3e-6*exp(-(y_for_UandT[i]/(3*111e3))^2)+3e-6*exp(-((y_for_UandT[i]-15*111e3)/(3*111e3))^2)+0.2e-6*exp(-((y_for_UandT[i]+15*111e3)/(3*111e3))^2)
    elseif day_of_year>155 && day_of_year<260
        S=-2.8e-6*exp(-(y_for_UandT[i]/(3*111e3))^2)+1.6e-6*exp(-((y_for_UandT[i]+15*111e3)/(3*111e3))^2)+1.6e-6*exp(-((y_for_UandT[i]-15*111e3)/(3*111e3))^2)
    else
        S=-1.6e-6*exp(-(y_for_UandT[i]/(3*111e3))^2)+0.6e-6*exp(-((y_for_UandT[i]+15*111e3)/(3*111e3))^2)+0.6e-6*exp(-((y_for_UandT[i]-15*111e3)/(3*111e3))^2)
    end
    return S*0.8
end

# First PDE about u:
function get_du(u,v,θ,timestamp,dt)
    adv=fill(0.0,y_Grids_for_V)
    adv_for_UandT=fill(0.0,y_Grids_for_UandT)
    u_grad=fill(0.0,y_Grids_for_V)
    v_grad=fill(0.0,y_Grids_for_UandT)
    Hvy=fill(0.0,y_Grids_for_UandT) # Heaviside of partial v / partial y
    beta=fill(0.0,y_Grids_for_V)
    beta_for_UandT=fill(0.0,y_Grids_for_UandT)
    # damping=fill(0.0,y_Grids_for_UandT)
    # Damping for u is not required.
    F=ϵᵤ*u
    for i=1:y_Grids_for_UandT
        v_grad[i]=(v[i+1]-v[i])/dy
        if i!=y_Grids_for_UandT
            u_grad[i+1]=(u[i+1]-u[i])/dy
        end
        adv=v.*u_grad
        beta[i+1]=β*y_for_V[i+1]*v[i+1]
        Hvy[i]=heaviside(v_grad[i])
    end
    for i=1:y_Grids_for_UandT
        adv_for_UandT[i]=(adv[i]+adv[i+1])/2
        beta_for_UandT[i]=(beta[i]+beta[i+1])/2
        if abs(y_for_UandT[i]-111e3*solar_insolation_degree) <= y₁
            S[i]=get_S(timestamp,i)
        end

    end
    return (beta_for_UandT-adv_for_UandT-Hvy.*v_grad.*u-F-S)*dt
end

# Second PDE about v
function get_dv(u,v,θ,dt)
    T=get_T(θ)
    beta=fill(0.0,y_Grids_for_V)
    T_grad=fill(0.0,y_Grids_for_V)
    damping=fill(0.0,y_Grids_for_V)
    for i=1:y_Grids_for_UandT-2
        T_grad[i+1]=(T[i+1]-T[i])/dy
        beta[i+1]=β*y_for_V[i+1]*(u[i]+u[i+1])*0.5
        damping[i+1]=(v[i+2]+v[i]-2*v[i+1])/(dy^2)
    end
    return (-beta-g*H*T_grad/T₀+kᵥ*damping)*0.5*dt
end

# Third PDE about θ
function get_dθ(u,v,θ,timestamp,dt)
    v_grad=fill(0.0,y_Grids_for_UandT)
    if ITCZ_displaced==ture
        solar_insolation_degree=23.5*sin(2*pi*(timestamp%365-79)/365)+6
    else
        solar_insolation_degree=23.5*sin(2*pi*(timestamp%365-79)/365)
    end
    for i=1:y_Grids_for_UandT
        if abs(y_for_UandT[i]-111e3*solar_insolation_degree) <= y₁
            v_grad[i]=(v[i+1]-v[i])/dy
        else
            v_grad[i]=0
        end
    end
    θₑ=get_θₑ(timestamp)
    return ((θₑ.-θ)/τ-δ*Δz*v_grad/H)*dt  # need to prescribe θₑ
end
# ----------------------------------------------

# ----------------------------------------------
# Define function for integration
# Function of Eular scheme. Only used for the first calculation.
function Euler(u,v,θ,get_du,get_dv,get_dθ,timestamp,dt)
    return u+get_du(u,v,θ,timestamp,dt), v+get_dv(u,v,θ,dt), θ+get_dθ(u,v,θ,timestamp,dt)
end

# leapfrog function.
# a_n+1 = a_n-1 + F(a_n) * 2 * dt
function leapfrog(u,u_before,v,v_before,θ,θ_before,get_du,get_dv,get_dθ,timestamp,dt)
    return u_before+get_du(u,v,θ,timestamp,2*dt), v_before+get_dv(u,v,θ,2*dt), θ_before+get_dθ(u,v,θ,timestamp,2*dt)
end
# ----------------------------------------------

# ----------------------------------------------
# core code for the integration
timesteps=Int(86400*365*integration_year_number/dt)

u_thisstep=u[end,:]
v_thisstep=v[end,:]
θ_thisstep=θ[end,:]
u_before, v_before, θ_before=Euler(u_thisstep,v_thisstep,θ_thisstep,get_du,get_dv,get_dθ,timestamp,-dt)
for i = 1:timesteps
    u_thisstep[1]=0
    u_thisstep[end]=0
    v_thisstep[1]=0
    v_thisstep[end]=0
    global timestamp
    u_after,v_after,θ_after=leapfrog(u_thisstep,u_before,v_thisstep,v_before,θ_thisstep,θ_before,get_du,get_dv,get_dθ,timestamp,dt)
    global u_before=u_thisstep+0.04*(u_after+u_before-2*u_thisstep)
    global v_before=v_thisstep+0.04*(v_after+v_before-2*v_thisstep)
    global θ_before=θ_thisstep+0.04*(θ_after+θ_before-2*θ_thisstep)
    global u_thisstep=u_after
    global v_thisstep=v_after
    global θ_thisstep=θ_after
    time_end=time[end]
    if isempty(u)
        u=reshape(u_thisstep,1,y_Grids_for_UandT)
        v=reshape(v_thisstep,1,y_Grids_for_V)
        θ=reshape(θ_thisstep,1,y_Grids_for_UandT)
    else
        append!(time,[time_end+dt])
        global u=vcat(u,reshape(u_thisstep,1,y_Grids_for_UandT))
        global v=vcat(v,reshape(v_thisstep,1,y_Grids_for_V))
        global θ=vcat(θ,reshape(θ_thisstep,1,y_Grids_for_UandT))
    end

    if mod((i+1)*dt,86400)==0
        stamp=lpad(Int((i+1)*dt/(86400)),5,"0")
        println("Day $stamp finished.")
        timestamp+=1
        if i>= (integration_year_number-1)*timesteps/integration_year_number
            using NCDatasets
            ds_u = NCDataset("$Dir/u_output_Day$stamp.nc","c")
            ds_v = NCDataset("$Dir/v_output_Day$stamp.nc","c")
            ds_θ = NCDataset("$Dir/pot.temp_Day$stamp.nc","c")
            # Define the dimension
            defDim(ds_u,"time",Inf)
            defDim(ds_u,"y",y_Grids_for_UandT)
            defDim(ds_v,"time",Inf)
            defDim(ds_v,"y",y_Grids_for_V)
            defDim(ds_θ,"time",Inf)
            defDim(ds_θ,"y",y_Grids_for_UandT)

            # Define a global attribute, the long name of the file.
            ds_u.attrib["title"] = "Zonal wind speed"
            ds_v.attrib["title"] = "Meridional wind speed"
            ds_θ.attrib["title"] = "Potential temperature"

            # Define the variables velocity
            variable_u = defVar(ds_u,"U",Float64,("y","time"))
            variable_uy= defVar(ds_u,"y",Float64,("y",))
            variable_ut= defVar(ds_u,"time",Float64,("time",))
            variable_v = defVar(ds_v,"V",Float64,("y","time"))
            variable_vy= defVar(ds_v,"Y",Float64,("y",))
            variable_vt= defVar(ds_v,"time",Float64,("time",))
            variable_θ = defVar(ds_θ,"θ",Float64,("y","time"))
            variable_θy= defVar(ds_θ,"Y",Float64,("y",))
            variable_θt= defVar(ds_θ,"time",Float64,("time",))

            # write a the complete data set
            variable_u[:,:] = permutedims(u, [2,1])
            variable_ut[:]=time
            variable_uy[:]=y_for_UandT
            variable_v[:,:] = permutedims(v, [2,1])
            variable_vt[:]=time
            variable_vy[:]=y_for_V
            variable_θ[:,:] = permutedims(θ, [2,1])
            variable_θt[:]=time
            variable_θy[:]=y_for_UandT

            # write attributes
            variable_u.attrib["units"] = "m/s"
            variable_v.attrib["units"] = "m/s"
            variable_θ.attrib["units"] = "K"
            variable_uy.attrib["units"]= "m"
            variable_vy.attrib["units"]= "m"
            variable_θy.attrib["units"]= "m"
            variable_ut.attrib["units"]= "seconds since 0000-01-01 00:00:00.0"
            variable_vt.attrib["units"]= "seconds since 0000-01-01 00:00:00.0"
            variable_θt.attrib["units"]= "seconds since 0000-01-01 00:00:00.0"
            variable_ut.attrib["calendar"]= "noleap"
            variable_vt.attrib["calendar"]= "noleap"
            variable_θt.attrib["calendar"]= "noleap"

            close(ds_u)
            close(ds_v)
            close(ds_θ)
        end

        global time=[time_end+2*dt]
        global u=[]
        global v=[]
        global θ=[]
    end
end
# ----------------------------------------------

timeend=time[end]/(86400*365)
println("We have already run $timeend years.")

max_stars_repo_name: zpcllyj/SobelSchneiderModel

max_issues_repo_name: zpcllyj/SobelSchneiderModel

max_forks_repo_name: zpcllyj/SobelSchneiderModel


text:
@theme TangoTheme Dict(
    :name => "Tango",
    :description => "A theme inspired by the Tango Icon Theme Guidelines.",
    :comments => "Based on Tango theme from Pygments.",
    :style => S"fg: f8f8f8",
    :tokens => Dict(
        COMMENT                => S"italic; fg: 8f5902",
        COMMENT_MULTILINE      => S"italic; fg: 8f5902",
        COMMENT_PREPROC        => S"italic; fg: 8f5902",
        COMMENT_SINGLE         => S"italic; fg: 8f5902",
        COMMENT_SPECIAL        => S"italic; fg: 8f5902",

        ERROR                  => S"fg: a40000; bg: ef2929",

        GENERIC                => S"fg: 000000",
        GENERIC_DELETED        => S"fg: a40000",
        GENERIC_EMPH           => S"italic; fg: 000000",
        GENERIC_ERROR          => S"fg: ef2929",
        GENERIC_HEADING        => S"bold; fg: 000080",
        GENERIC_INSERTED       => S"fg: 00a000",
        GENERIC_OUTPUT         => S"italic; fg: 000000",
        GENERIC_PROMPT         => S"fg: 8f5902",
        GENERIC_STRONG         => S"bold; fg: 000000",
        GENERIC_SUBHEADING     => S"bold; fg: 800080",
        GENERIC_TRACEBACK      => S"bold; fg: a40000",

        KEYWORD                => S"bold; fg: 204a87",
        KEYWORD_CONSTANT       => S"bold; fg: 204a87",
        KEYWORD_DECLARATION    => S"bold; fg: 204a87",
        KEYWORD_NAMESPACE      => S"bold; fg: 204a87",
        KEYWORD_PSEUDO         => S"bold; fg: 204a87",
        KEYWORD_RESERVED       => S"bold; fg: 204a87",
        KEYWORD_TYPE           => S"bold; fg: 204a87",

        LITERAL                => S"fg: 000000",
        LITERAL_DATE           => S"fg: 000000",

        NAME                   => S"fg: 000000",
        NAME_ATTRIBUTE         => S"fg: c4a000",
        NAME_BUILTIN           => S"fg: 204a87",
        NAME_BUILTIN_PSEUDO    => S"fg: 3465a4",
        NAME_CLASS             => S"fg: 000000",
        NAME_CONSTANT          => S"fg: 000000",
        NAME_DECORATOR         => S"bold; fg: 5c35cc",
        NAME_ENTITY            => S"fg: ce5c00",
        NAME_EXCEPTION         => S"bold; fg: cc0000",
        NAME_FUNCTION          => S"fg: 000000",
        NAME_LABEL             => S"fg: f57900",
        NAME_NAMESPACE         => S"fg: 000000",
        NAME_OTHER             => S"fg: 000000",
        NAME_PROPERTY          => S"fg: 000000",
        NAME_TAG               => S"bold; fg: 204a87",
        NAME_VARIABLE          => S"fg: 000000",
        NAME_VARIABLE_CLASS    => S"fg: 000000",
        NAME_VARIABLE_GLOBAL   => S"fg: 000000",
        NAME_VARIABLE_INSTANCE => S"fg: 000000",

        NUMBER                 => S"bold; fg: 0000cf",
        NUMBER_FLOAT           => S"bold; fg: 0000cf",
        NUMBER_HEX             => S"bold; fg: 0000cf",
        NUMBER_INTEGER         => S"bold; fg: 0000cf",
        NUMBER_INTEGER_LONG    => S"bold; fg: 0000cf",
        NUMBER_OCT             => S"bold; fg: 0000cf",

        OPERATOR               => S"bold; fg: ce5c00",
        OPERATOR_WORD          => S"bold; fg: 204a87",

        OTHER                  => S"fg: 000000",

        PUNCTUATION            => S"bold; fg: 000000",

        STRING                 => S"fg: 4e9a06",
        STRING_BACKTICK        => S"fg: 4e9a06",
        STRING_CHAR            => S"fg: 4e9a06",
        STRING_DOC             => S"italic; fg: 8f5902",
        STRING_DOUBLE          => S"fg: 4e9a06",
        STRING_ESCAPE          => S"fg: 4e9a06",
        STRING_HEREDOC         => S"fg: 4e9a06",
        STRING_INTERPOL        => S"fg: 4e9a06",
        STRING_OTHER           => S"fg: 4e9a06",
        STRING_REGEX           => S"fg: 4e9a06",
        STRING_SINGLE          => S"fg: 4e9a06",
        STRING_SYMBOL          => S"fg: 4e9a06",

        TEXT                   => S"",

        WHITESPACE             => S"underline; fg: f8f8f8",
    )
)


text:
# Plot text labels at arbitrary coordinate points
function PlotText(text::String, x::Real, y::Real, vertical::Bool=false,
                  pixel_offset::ImVec2=ImVec2(0, 0))
    return LibCImPlot.PlotText(text, x, y, vertical, pixel_offset)
end

function PlotImage(user_texture_id, bounds_min::ImPlotPoint=(0, 0),
                   bounds_max::ImPlotPoint=(0, 0), uv0::ImVec2=(0, 0), uv1::ImVec2=(1, 1),
                   tint_col::ImVec4=(1, 1, 1, 1); label_id::String="")
    return LibCImPlot.PlotImage(label_id, user_texture_id, bounds_min, bounds_max, uv0, uv1,
                                tint_col)
end

function PlotImage(label_id::String, user_texture_id, bounds_min::ImPlotPoint=(0, 0),
                   bounds_max::ImPlotPoint=(0, 0), uv0::ImVec2=(0, 0), uv1::ImVec2=(1, 1),
                   tint_col::ImVec4=(1, 1, 1, 1))
    return LibCImPlot.PlotImage(label_id, user_texture_id, bounds_min, bounds_max, uv0, uv1,
                                tint_col)
end


text:
import PyCall, Pkg

pip = joinpath(split(PyCall.PYTHONHOME, ":")[end], "bin/pip")

if !isfile(pip)
    println("`pip` is not available in the current PyCall.jl")
    println("Configuring PyCall.jl to use Conda.jl")
    ENV["PYTHON"] = ""
    Pkg.build("PyCall")
end


text:
"""
    BagNode{T <: Union{AbstractMillNode, Missing}, B <: AbstractBags, C} <: AbstractBagNode

Data node that represents a multi-instance learning problem.
Contains instances stored in a subtree of type `T`,
bag indices of type `B` and optional metadata of type `C`.

See also: [`WeightedBagNode`](@ref), [`AbstractBagNode`](@ref),
    [`AbstractMillNode`](@ref), [`BagModel`](@ref).
"""
struct BagNode{T <: Maybe{AbstractMillNode}, B <: AbstractBags, C} <: AbstractBagNode
    data::T
    bags::B
    metadata::C

    function BagNode(d::T, b::B, m::C=nothing) where {T <: Maybe{AbstractMillNode}, B <: AbstractBags, C}
        ismissing(d) && any(length.(b) .> 0) && error("BagNode with `missing` in data cannot have a non-empty bag")
        new{T, B, C}(d, b, m)
    end
end

"""
    BagNode(d::Union{AbstractMillNode, Missing}, b::AbstractBags, m=nothing)
    BagNode(d::Union{AbstractMillNode, Missing}, b::AbstractVector, m=nothing)

Construct a new [`BagNode`](@ref) with data `d`, bags `b`, and metadata `m`.
If `b` is an `AbstractVector`, [`Mill.bags`](@ref) is applied first.

# Examples
```jldoctest
julia> BagNode(ArrayNode(maybehotbatch([1, missing, 2], 1:2)), AlignedBags([1:1, 2:3]))
BagNode with 2 obs
  └── ArrayNode(2×3 MaybeHotMatrix with Union{Missing, Bool} elements)

julia> BagNode(ArrayNode(randn(2, 5)), [1, 2, 2, 1, 1])
BagNode with 2 obs
  └── ArrayNode(2×5 Array with Float64 elements)
```

See also: [`WeightedBagNode`](@ref), [`AbstractBagNode`](@ref),
    [`AbstractMillNode`](@ref), [`BagModel`](@ref).
"""
BagNode(d::Maybe{AbstractMillNode}, b::AbstractVector, m=nothing) = BagNode(d, bags(b), m)

Flux.@functor BagNode

mapdata(f, x::BagNode) = BagNode(mapdata(f, x.data), x.bags, x.metadata)

dropmeta(x::BagNode) = BagNode(x.data, x.bags)

function Base.getindex(x::BagNode, i::VecOrRange{<:Int})
    nb, ii = remapbags(x.bags, i)
    emptyismissing() && isempty(ii) && return(BagNode(missing, nb, nothing))
    BagNode(subset(x.data,ii), nb, subset(x.metadata, i))
end

function reduce(::typeof(catobs), as::Vector{<:BagNode})
    d = filter(!ismissing, data.(as))
    md = filter(!isnothing, metadata.(as))
    bags = reduce(vcat, [n.bags for n in as])
    BagNode(reduce(catobs, d), bags, reduce(catobs, md))
end

removeinstances(a::BagNode, mask) = BagNode(subset(a.data, findall(mask)), adjustbags(a.bags, mask), a.metadata)

Base.hash(e::BagNode, h::UInt) = hash((e.data, e.bags, e.metadata), h)
(e1::BagNode == e2::BagNode) = e1.data == e2.data && e1.bags == e2.bags && e1.metadata == e2.metadata
Base.isequal(e1::BagNode, e2::BagNode) = isequal(e1.data, e2.data) && isequal(e1.bags, e2.bags) && isequal(e1.metadata, e2.metadata)


text:
# # interface for skipping null entries

# function skipna_init(f, op, X::DataValueArray,
#                        ifirst::Int, ilast::Int)
#     # Get first non-null element
#     ifirst = Base.findnext(x -> x == false, X.isna, ifirst)
#     @inbounds v1 = X.values[ifirst]

#     # Get next non-null element
#     ifirst = Base.findnext(x -> x == false, X.isna, ifirst + 1)
#     @inbounds v2 = X.values[ifirst]

#     # Reduce first two elements
#     return op(f(v1), f(v2)), ifirst
# end

# # sequential map-reduce
# function mapreduce_seq_impl_skipna(f, op, X::DataValueArray,
#                                      ifirst::Int, ilast::Int)
#     # initialize first reduction
#     v, i = skipna_init(f, op, X, ifirst, ilast)

#     while i < ilast
#         i += 1
#         @inbounds isna = X.isna[i]
#         isna && continue
#         @inbounds entry = X.values[i]
#         v = op(v, f(entry))
#     end
#     return DataValue(v)
# end

# # pairwise map-reduce
# function mapreduce_pairwise_impl_skipna{T}(f, op, X::DataValueArray{T},
#                                              ifirst::Int, ilast::Int,
#                                             #  n_notnull::Int, blksize::Int)
#                                             blksize::Int)
#     if ifirst + blksize > ilast
#         # fall back to Base implementation if no nulls in block
#         # if any(isna, slice(X, ifirst:ilast))
#             return mapreduce_seq_impl_skipna(f, op, X, ifirst, ilast)
#         # else
#             # DataValue(Base.mapreduce_seq_impl(f, op, X.values, ifirst, ilast))
#         # end
#     else
#         imid = (ifirst + ilast) >>> 1
#         # n_notnull1 = imid - ifirst + 1 - countnz(X.isna[ifirst:imid])
#         # n_notnull2 = ilast - imid - countnz(X.isna[imid+1:ilast])
#         v1 = mapreduce_pairwise_impl_skipna(f, op, X, ifirst, imid,
#                                               blksize)
#         v2 = mapreduce_pairwise_impl_skipna(f, op, X, imid+1, ilast,
#                                               blksize)
#         return op(v1, v2)
#     end
# end

# # from comment: https://github.com/JuliaLang/julia/pull/16217#issuecomment-223768129
# sum_pairwise_blocksize(T) = Base.pairwise_blocksize(T, +)

# mapreduce_impl_skipna{T}(f, op, X::DataValueArray{T}) =
#     mapreduce_seq_impl_skipna(f, op, X, 1, length(X.values))
# mapreduce_impl_skipna(f, op::typeof(+), X::DataValueArray) =
#     mapreduce_pairwise_impl_skipna(f, op, X, 1, length(X.values),
#                                    max(128, sum_pairwise_blocksize(f)))

# # general mapreduce interface

# function _mapreduce_skipna{T}(f, op, X::DataValueArray{T}, missingdata::Bool)
#     n = length(X)
#     !missingdata && return DataValue(Base.mapreduce_impl(f, op, X.values, 1, n))

#     nnull = countnz(X.isna)
#     nnull == n && return Base.mr_empty(f, op, T)
#     nnull == n - 1 && return Base.r_promote(op, f(X.values[findnext(x -> x == false), X, 1]))
#     # nnull == 0 && return Base.mapreduce_impl(f, op, X, 1, n)

#     return mapreduce_impl_skipna(f, op, X)
# end

# function Base._mapreduce(f, op, X::DataValueArray, missingdata)
#     missingdata && return Base._mapreduce(f, op, X)
#     DataValue(Base._mapreduce(f, op, X.values))
# end

# # to fix ambiguity warnings
# function Base.mapreduce(f, op::Union{typeof(&), typeof(|)},
#                         X::DataValueArray, skipna::Bool = false)
#     missingdata = any(isna, X)
#     if skipna
#         return _mapreduce_skipna(f, op, X, missingdata)
#     else
#         return Base._mapreduce(f, op, X, missingdata)
#     end
# end

# const specialized_binary = identity

"""
    mapreduce(f, op::Function, X::DataValueArray; [skipna::Bool=false])

Map a function `f` over the elements of `X` and reduce the result under the
operation `op`. One can set the behavior of this method to skip the null entries
of `X` by setting the keyword argument `skipna` equal to true. If `skipna`
behavior is enabled, `f` will be automatically lifted over the elements of `X`.
Note that, in general, mapreducing over a `DataValueArray` will return a
`DataValue` object regardless of whether `skipna` is set to `true` or `false`.
"""
function Base.mapreduce(f, op::Function, X::T; skipna::Bool = false) where {N,S<:DataValue,T<:AbstractArray{S,N}}
    if skipna
        return DataValue(mapreduce(f, op, dropna(X)))
    else
        return Base._mapreduce(f, op, IndexStyle(X), X)
    end
end

"""
    reduce(f, op::Function, X::DataValueArray; [skipna::Bool=false])

Reduce `X`under the operation `op`. One can set the behavior of this method to
skip the null entries of `X` by setting the keyword argument `skipna` equal
to true. If `skipna` behavior is enabled, `f` will be automatically lifted
over the elements of `X`. Note that, in general, mapreducing over a
`DataValueArray` will return a `DataValue` object regardless of whether `skipna`
is set to `true` or `false`.
"""
function Base.reduce(op, X::T; skipna::Bool = false) where {N,S<:DataValue,T<:AbstractArray{S,N}}
    return mapreduce(identity, op, X; skipna = skipna)
end

# standard reductions

for (fn, op) in ((:(Base.sum), +),
                 (:(Base.prod), *),
                 (:(Base.minimum), min),
                 (:(Base.maximum), max))
    @eval begin
        # supertype(typeof(@functorize(abs))) returns Func{1} on Julia 0.4,
        # and Function on 0.5
        function $fn(f::Union{Function,Type},X::T; skipna::Bool=false) where {N,S<:DataValue,T<:AbstractArray{S,N}}
            return mapreduce(f, $op, X; skipna = skipna)
        end

        function $fn(X::T; skipna::Bool=false) where {N,S<:DataValue,T<:AbstractArray{S,N}}
            return mapreduce(identity, $op, X; skipna = skipna)
        end
    end
end

for op in (Base.min, Base.max)
    @eval begin
        function Base._mapreduce(::typeof(identity), ::$(typeof(op)),
                                    X::DataValueArray{T}, missingdata) where {T}
            missingdata && return DataValue{T}()
            DataValue(Base._mapreduce(identity, $op, X.values))
        end
    end
end

# function Base.mapreduce_impl{T}(f, op::typeof(Base.scalarmin), X::DataValueArray{T},
#                                 first::Int, last::Int)
#     i = first
#     v = f(X[i])
#     i += 1
#     while i <= last
#         @inbounds x = f(X[i])
#         if isna(x) | isna(v)
#             return DataValue{eltype(x)}()
#         elseif x.value < v.value
#             v = x
#         end
#         i += 1
#     end
#     return v
# end

# function Base.mapreduce_impl{T}(f, op::typeof(Base.scalarmax), X::DataValueArray{T},
#                                 first::Int, last::Int)
#     i = first
#     v = f(X[i])
#     i += 1
#     while i <= last
#         @inbounds x = f(X[i])
#         if isna(x) | isna(v)
#             return DataValue{eltype(x)}()
#         elseif x.value > v.value
#             v = x
#         end
#         i += 1
#     end
#     return v
# end

function Base.extrema(X::T; skipna::Bool = false) where {N,T2,T<:DataValueArray{T2,N}}
    length(X) > 0 || throw(ArgumentError("collection must be non-empty"))
    vmin = DataValue{T2}()
    vmax = DataValue{T2}()
    @inbounds for i in 1:length(X)
        x = X.values[i]
        missing = X.isna[i]
        if skipna && missing
            continue
        elseif missing
            return (DataValue{T2}(), DataValue{T2}())
        elseif isna(vmax) # Equivalent to isna(vmin)
            vmax = vmin = DataValue(x)
        elseif x > vmax.value
            vmax = DataValue(x)
        elseif x < vmin.value
            vmin = DataValue(x)
        end
    end
    return (vmin, vmax)
end

function Base.extrema(X::T; skipna::Bool = false) where {N,T2,S<:DataValue{T2},T<:AbstractArray{S,N}}
    length(X) > 0 || throw(ArgumentError("collection must be non-empty"))
    vmin = DataValue{T2}()
    vmax = DataValue{T2}()
    @inbounds for i in 1:length(X)
        missing = isna(X[i])

        if skipna && missing
            continue
        elseif missing
            return (DataValue{T2}(), DataValue{T2}())
        else
            x = get(X[i])
            if isna(vmax) # Equivalent to isna(vmin)
                vmax = vmin = DataValue(x)
            elseif x > vmax.value
                vmax = DataValue(x)
            elseif x < vmin.value
                vmin = DataValue(x)
            end
        end
    end
    return (vmin, vmax)
end


text:
"""
PlasmaProperties, a simple plasma calculator module in Julia

The module contains basic functions to compute common plasma properties, such as characteristic frequencies and lengths.

It builds on Unitful, which provides a simple unit system for physical quantities.
"""
module PlasmaProperties

using Unitful
import Unitful: Length, Time, Mass, Energy, Temperature, BField, EField
using PhysicalConstants
ε0 = PhysicalConstants.CODATA2018.ε_0 |> u"F/m"
μ0 = PhysicalConstants.CODATA2018.μ_0 |> u"N/A^2"
me = PhysicalConstants.CODATA2018.m_e |> u"kg"
qe = PhysicalConstants.CODATA2018.e |> u"C"
kB = PhysicalConstants.CODATA2018.k_B |> u"J/K"
mu = PhysicalConstants.CODATA2018.m_u |> u"kg"

export QuasineutralPlasma # Plasma object type and constructors
export atomic_masses # Atomic masses
export f_ce, f_pe, r_Le, r_Lsi, λ_De, c_s, c_the, ν_ei, ν_ie, ν_ee, χ # Calculator functions
export ε0, μ0, me, qe, kB # Physical constants, reexported from PhysicalConstants

include("auxiliary.jl")
include("types.jl")
include("functions.jl")

end # module


text:
function aicSSM(m)
	# Akaike Information Criterion for State Space Models
	nx = size(m.A,1)
	ny = size

	P = ny*nx + nx*nx + nx*(nx+1)/2 + ny*(ny+1)/2 + nx*nu + ny*nu
	aic[k] = -2*loglik[end] + 2*P

end


text:
## Calibrate camera
    # Initial guess at reference camera

    # Richard Shenton email:
    # Camera height: 2.065m
    #
    # Camera offset from centre line: 0.6m
    # Pitch: 19.16 degrees down
    # Yaw: 7.28 degrees to the left
    # Roll: Assumed to be zero

    # The camera frame rate is 25fps for this data. Regarding train details, there are some details here https://en.wikipedia.org/wiki/British_Rail_Class_357 .
    #
    # It looks to me that the carriage with the cab (seems to be called DMOS) is 20.75(L) x 2.8 (W) x 3.78(H) and weighs 40.7 t.

using TrainView

using CSV, DataFrames
using Statistics, LinearAlgebra

file = "../data/rail_detection_results_20210127.csv";
df = CSV.read(file, DataFrame);

Lu1 = df." Left_Rail_Bottom_X";
Lv1 = df." Left_Rail_Bottom_Y";
Lu2 = df." Left_Rail_Top_X";
Lv2 = df." Left_Rail_Top_Y";

Ru1 = df." Right_Rail_Bottom_X";
Rv1 = df." Right_Rail_Bottom_Y";
Ru2 = df." Right_Rail_Top_X";
Rv2 = df." Right_Rail_Top_Y";

# max_v = max(maximum(Lv1),maximum(Rv1))
max_v = 720
vs = [Lv1, Lv2, Rv1, Rv2];
for v in vs
    v[1:end] = v .- (max_v / 2.0)
end

# max_u = max(maximum(Lu2),maximum(Ru2))
max_u = 1280
us = [Lu1, Lu2, Ru1, Ru2];
for v in us
    v[1:end] = v .- (max_u / 2.0)
end

Lu = LinRange(mean(Lu1),mean(Lu2),20);
Lv = LinRange(mean(Lv1),mean(Lv2),20);
Ru = LinRange(mean(Ru1),mean(Ru2),20);
Rv = LinRange(mean(Rv1),mean(Rv2),20);

left_uvs =  [[Lu[i],Lv[i]] for i in eachindex(Lu)];
right_uvs = [[Ru[i],Rv[i]] for i in eachindex(Ru)];

    trackprop = TrackProperties(track_gauge = 1.435 + 0.065)

    cameraposition_reference = [0.0,0.6,-2.165];
    ψθφ_ref = [0.0,-19.1,-7.28] .* (pi/180.0)

    camera_reference = VideoCamera(cameraposition_reference;
        focal_length = 5.8e-3,
        # focal_length = 5.8e-3,
        pixelspermeter = 1 / 5.5e-6,
        ψθφ = ψθφ_ref
    )

    uLs = left_track_image_u(camera_reference, trackprop, Lv);
    uRs = right_track_image_u(camera_reference, trackprop, Rv);

    distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera_reference, trackprop;
        choose_distortions = [:Y,:Z,:θ,:φ],
        # choose_distortions = [:Y,:Z,:f,:θ,:φ,:ψ,:α,:β],
        iterations = 6)

    camera = VideoCamera(camera_reference,distortion)

    uLs2 = left_track_image_u(camera, trackprop, Lv);
    uRs2 = right_track_image_u(camera, trackprop, Rv);

    using Plots
    gr(yflip = true)

    plot(Lu,Lv, lab = "")
    plot!(Ru,Rv, lab = "")

    plot!(uLs2,Lv,linestyle=:dash, lab = "")
    plot!(uRs2,Rv,linestyle=:dash, lab = "")

## Analise journey
    choose_distortions = [:Y,:Z,:θ,:ψ];
    # choose_distortions = [:Y,:Z,:θ,:φ,:ψ,:α,:β]

    distortions = map(eachindex(Lu1)) do j
        Lu = LinRange(Lu1[j],Lu2[j],20);
        Lv = LinRange(Lv1[j],Lv2[j],20);
        Ru = LinRange(Ru1[j],Ru2[j],20);
        Rv = LinRange(Rv1[j],Rv2[j],20);

        left_uvs =  [[Lu[i],Lv[i]] for i in eachindex(Lu)];
        right_uvs = [[Ru[i],Rv[i]] for i in eachindex(Ru)];

        distortion = rail_uvs_to_distortion(left_uvs, right_uvs, camera, trackprop;
            choose_distortions = choose_distortions,
            iterations = 1
        );
        return [distortion[d] for d in choose_distortions]
    end

    data = transpose(hcat(distortions...))

    dict = Dict(choose_distortions[j] => data[:,j] for j in eachindex(choose_distortions))

    df = DataFrame(dict);

    moving_average(vs,n) = [sum(@view vs[i:(i+n-1)])/n for i in 1:(length(vs)-(n-1))]

    ## Plotting using the @df macro specifying colum names as symbol:
    using StatsPlots
    n = 15;
    n = 25;
    fps = 25;
    time = (1:size(data,1)) ./ fps
    df.time = time;

    @df df plot(moving_average(:time,n), [moving_average(:Y,n),moving_average(:Z,n)], labels = ["Horizontal" "Vertical"], xlab = "seconds")

    @df df plot(:time,[:Y,:Z], labels = ["Horizontal" "Vertical"], xlab = "seconds")

    @df df plot(moving_average(:time,n), [moving_average(:θ,n) * 180/π, moving_average(:ψ,n) * 180/π], labels = ["Pitch" "Roll"], xlab = "seconds", ylab = "Degree")


text:
module sigma_estimate

function sigma_est(point_array::AbstractMatrix{<:Real},
    center::AbstractMatrix{<:Real},
    fruit_array::AbstractMatrix{<:Real})
    n = size(point_array)[1]
    #println(n)
    #println(point_array[1,2])
    sum_x = 0
    sum_y = 0
    for i in 1:n
        if ((point_array[i, 1]-center[1])^2+(point_array[i, 2]-center[2])^2) < 45
            sum_x = sum_x + (point_array[i, 1]-center[1])^2*fruit_array[i]
            sum_y = sum_y + (point_array[i, 2]-center[2])^2*fruit_array[i]
        end
    end
    t = sum(fruit_array)
    esigma_x = (sum_x+1)/(t+1)
    esigma_y = (sum_y+1)/(t+1)
    return esigma_x, esigma_y
end

#test code
#point_array = [4 3; 6 5; 3 3]
#center = [5 4]
#fruit_array = [2 1 3]
#sigma_x, sigma_y = sigma_estimate(point_array, center, fruit_array)
#println(sigma_x, sigma_y)

end


text:
struct Jupyter <: Format end

extensionformat(::Val{:ipynb}) = Jupyter()

function parse(io::IO, format::Jupyter)
    return parse(JSON3.read(io), format)
end

function parse(obj::JSON3.Object, format::Jupyter)
    cs = XTree[]
    lang = get(get(get(obj, :metadata, Dict()), :language_info, Dict()), :name, "")
    for cell in obj[:cells]
        cs = vcat(cs, children(parsejupytercell(cell, lang)))
    end
    return XNode(
        :body,
        cs
    )
end

function parsejupytercell(cell, lang = "")
    type = cell[:cell_type]
    if  type == "markdown"
        return parsejupytercellmd(cell)
    elseif type == "code"
        return parsejupytercellcode(cell, lang)
    else
        error("Unsupported cell type $type.")
    end
end

function parsejupytercellmd(cell)
    return parse(join(cell[:source], "\n"), Markdown())
end

function parsejupytercellcode(cell, lang)
    code = join(cell[:source])
    xcode = XNode(:pre, Dict(:lang => lang), [XNode(:code, [XLeaf(code)])])
    cs = XTree[xcode]

    outputs = cell[:outputs]
    stream = ""
    for output in outputs
        if output[:output_type] == "stream"
            stream *= join(output[:text])
        elseif output[:output_type] == "execute_result"
            # Add concatenated outputs first
            push!(cs, viewcodeoutput(stream))
            stream = ""

            # The handle result
            reprs = Dict(MIME(k) => join(v) for (k, v) in output[:data])
            push!(cs, viewcoderesult(PreRendered(reprs)))
        end
    end
    if stream != ""
        push!(cs, viewcodeoutput(stream))
    end

    return XNode(
        :div,
        Dict(:class => "cellcontainer"),
        cs,
    )
end

struct PreRendered
    reprs::Dict
end

Base.show(io::IO, prerendered::PreRendered) = print(io,
    "Prerendered() with $(length(prerendered.reprs)) reprs")

function render!(io, x::XLeaf{PreRendered}, ::HTML)
    reprs = x[].reprs
    for mime in HTML_MIMES
        if mime in keys(reprs)
            print(io, adapthtmlstr(mime, reprs[mime]))
            return
        end
    end
    error("Could not find mime for $(x[])!")
end

Base.showable(mime::MIME, x::PreRendered) = mime in keys(x.reprs)


text:
import YAML, MPI, OrderedCollections

function default_config_filename()
    return "riopa_default.yaml"
end

const Config = OrderedCollections.LittleDict{Symbol,Any}

function read_config(filename::AbstractString = default_config_filename())
    return YAML.load_file(filename, dicttype = Config)
end

read_config(::Nothing) = read_config()

function default_config()
    D = Config
    config = D(
        :datasets => [
            D(
                :type => "output",
                :name => "data 1",
                :basename => "data_one",
                :io_backend => "HDF5",
                :nsteps => 10,
                :compute_seconds => 1.0,
                :data_streams => [
                    D(
                        :name => "Level_0",
                        :evolution => "none",
                        :nprocs_ratio => 0.5,
                        :proc_payloads => [
                            D(:size_range => [1000, 1200], :ratio => 0.5),
                            D(:size_range => [2000, 2400], :ratio => 0.5),
                        ],
                    ),
                    D(
                        :name => "Level_1",
                        :proc_payloads => [
                            D(:size_range => [2000, 2500], :ratio => "1/4"),
                            D(:size_range => [4000, 4800], :ratio => "3/4"),
                        ],
                    ),
                ],
            ),
            D(
                :type => "output",
                :name => "data 2",
                :io_backend => "HDF5",
                :basename => "data_two",
                :nsteps => 3,
                :compute_seconds => 3.0,
                :data_streams => [
                    D(
                        :name => "Level_0",
                        :proc_payloads => [
                            D(:size_range => [1000, 1200], :ratio => 0.1),
                            D(:size_range => [2000, 2400], :ratio => 0.9),
                        ],
                    ),
                    D(
                        :name => "Level_1",
                        :proc_payloads => [
                            D(:size_range => [2000, 2500], :ratio => "1/8"),
                            D(:size_range => [4000, 4800], :ratio => "7/8"),
                        ],
                    ),
                ],
            ),
        ],
    )
    return config
end

function generate_config(filename::AbstractString = default_config_filename())
    if MPI.Comm_rank(MPI.COMM_WORLD) == 0
        YAML.write_file(filename, default_config())
        println("Generated config file: ", filename)
    end
    nothing
end

generate_config(::Nothing) = generate_config()


text:
using Test
using FileIO, Random, Pkg
using MakieGallery
using Makie, AbstractPlotting
using Statistics

database = MakieGallery.load_database()

# which one are the slowest and kicked those out!
slow_examples = [
    "Animated time series",
    "Animation",
    "Lots of Heatmaps",
    "Chess Game",
    "Line changing colour",
    "Line changing colour with Observables",
    "Colormap collection",
    "Record Video",
    "Animated surface and wireframe",
    "Moire",
    "Line GIF",
    "Electrostatic repulsion",
    "pong",
    "pulsing marker",
    "Travelling wave",
    "Axis theming",
    "Legend",
    "Color Legend",
    "DifferentialEquations path animation",
    "Interactive Differential Equation",
    "Spacecraft from a galaxy far, far away",
    "Type recipe for molecule simulation",
    "WorldClim visualization",
    "Image on Geometry (Moon)",
    "Image on Geometry (Earth)",
]
# DIffeq errors with stackoverflow
# The others look fine on the CI, but the measured difference is too high -.-
# Maybe related to the axis changes, will investigate later
filter!(database) do entry
    !("diffeq" in entry.tags) &&
    !(entry.unique_name in (:analysis, :colormap_collection, :lots_of_heatmaps))
 end

# Download is broken on CI
if get(ENV, "CI", "false") == "true"
    printstyled("CI detected\n"; bold = true, color = :yellow)
    println("Filtering out examples which download")
    filter!(entry-> !("download" in entry.tags), database)
end

tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")

isdir(tested_diff_path) && rm(tested_diff_path, force = true, recursive = true)
mkpath(tested_diff_path)

isdir(test_record_path) && rm(test_record_path, force = true, recursive = true)
mkpath(test_record_path)

examples = MakieGallery.record_examples(test_record_path)

@test length(examples) == length(database)

# MakieGallery.generate_preview(test_record_path, joinpath(homedir(), "Desktop", "index.html"))
# MakieGallery.generate_thumbnails(test_record_path)
# MakieGallery.gallery_from_recordings(test_record_path, joinpath(test_record_path, "index.html"))

printstyled("Running ", color = :green, bold = true)
println("visual regression tests")

MakieGallery.run_comparison(test_record_path, tested_diff_path)


text:
include(joinpath(Sys.BINDIR, "..", "share", "julia", "stdlib", "v$(VERSION.major).$(VERSION.minor)", "REPL", "src", "latex_symbols.jl"));

φ = open("latex_symbols.py", "w")

println(φ, "from __future__ import unicode_literals\n\n")
println(φ, "latex_symbols = [")
for (ω, (α, β)) in enumerate(latex_symbols)
    print(φ, "  (\"", escape_string(α), "\", u\"",  β, "\")")
    ω < length(latex_symbols) && print(φ, ",")
    println(φ, "")
end
println(φ, "]")

close(φ)

max_stars_repo_name: albertas-jn/Julia-sublime

max_issues_repo_name: albertas-jn/Julia-sublime

max_forks_repo_name: albertas-jn/Julia-sublime


text:
using Test
using IonBase

@testset "command/install" begin
    include("command/install.jl")
end

@testset "command/activate" begin
    include("command/activate.jl")
end

@testset "command/create" begin
    include("command/create.jl")
end

@testset "command/release" begin
    include("command/release.jl")
end

@testset "command/search" begin
    include("command/search.jl")
end

@testset "command/pkg" begin
    include("command/pkg.jl")
end


text:
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.MPIStateArrays
using ClimateMachine.Atmos
using ClimateMachine.ConfigTypes
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers: ManyColumnLU
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.Mesh.Grids
using ClimateMachine.TemperatureProfiles
using ClimateMachine.TurbulenceClosures
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Geometry: LocalGeometry

using LinearAlgebra
using StaticArrays
using Test

using CLIMAParameters
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function init_to_ref_state!(problem, bl, state, aux, localgeo, t)
    FT = eltype(state)
    state.ρ = aux.ref_state.ρ
    state.ρu = SVector{3, FT}(0, 0, 0)
    state.energy.ρe = aux.ref_state.ρe
end

function config_balanced(
    FT,
    poly_order,
    temp_profile,
    numflux,
    (config_type, config_fun, config_args),
)
    ref_state = HydrostaticState(temp_profile; subtract_off = false)

    physics = AtmosPhysics{FT}(
        param_set;
        ref_state = ref_state,
        turbulence = ConstantDynamicViscosity(FT(0)),
        hyperdiffusion = NoHyperDiffusion(),
        moisture = DryModel(),
    )
    model = AtmosModel{FT}(
        config_type,
        physics;
        source = (Gravity(),),
        init_state_prognostic = init_to_ref_state!,
    )

    config = config_fun(
        "balanced state",
        poly_order,
        config_args...,
        param_set,
        nothing;
        model = model,
        numerical_flux_first_order = numflux,
    )

    return config
end

function main()
    FT = Float64
    poly_order = 4

    timestart = FT(0)
    timeend = FT(100)
    domain_height = FT(50e3)

    LES_params = let
        LES_resolution = ntuple(_ -> domain_height / 3poly_order, 3)
        LES_domain = ntuple(_ -> domain_height, 3)
        (LES_resolution, LES_domain...)
    end

    GCM_params = let
        GCM_resolution = (3, 3)
        (GCM_resolution, domain_height)
    end

    GCM = (AtmosGCMConfigType, ClimateMachine.AtmosGCMConfiguration, GCM_params)
    LES = (AtmosLESConfigType, ClimateMachine.AtmosLESConfiguration, LES_params)

    imex_solver_type = ClimateMachine.IMEXSolverType(
        splitting_type = HEVISplitting(),
        implicit_model = AtmosAcousticGravityLinearModel,
        implicit_solver = ManyColumnLU,
        solver_method = ARK2GiraldoKellyConstantinescu,
    )
    explicit_solver_type = ClimateMachine.ExplicitSolverType(
        solver_method = LSRK54CarpenterKennedy,
    )

    @testset for config in (LES, GCM)
        @testset for ode_solver_type in (explicit_solver_type, imex_solver_type)
            @testset for numflux in (
                CentralNumericalFluxFirstOrder(),
                RoeNumericalFlux(),
                HLLCNumericalFlux(),
            )
                @testset for temp_profile in (
                    IsothermalProfile(param_set, FT),
                    DecayingTemperatureProfile{FT}(param_set),
                )
                    driver_config = config_balanced(
                        FT,
                        poly_order,
                        temp_profile,
                        numflux,
                        config,
                    )

                    solver_config = ClimateMachine.SolverConfiguration(
                        timestart,
                        timeend,
                        driver_config,
                        Courant_number = FT(0.1),
                        init_on_cpu = true,
                        ode_solver_type = ode_solver_type,
                        CFL_direction = EveryDirection(),
                        diffdir = HorizontalDirection(),
                    )

                    Qinit = similar(solver_config.Q)
                    Qinit .= solver_config.Q

                    ClimateMachine.invoke!(solver_config)

                    error =
                        euclidean_distance(solver_config.Q, Qinit) / norm(Qinit)
                    @test error <= 100 * eps(FT)
                end
            end
        end
    end
end

main()


text:
using AOPE
using Documenter

DocMeta.setdocmeta!(AOPE, :DocTestSetup, :(using AOPE); recursive=true)

makedocs(;
    modules=[AOPE],
    authors="António Saragga Seabra",
    repo="https://github.com/ASaragga/AOPE.jl/blob/{commit}{path}#{line}",
    sitename="AOPE.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ASaragga.github.io/AOPE.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/ASaragga/AOPE.jl",
)


text:
#=
Set
=#

using SimpleDataStructures

ss = SimpleSet{Int}()
push!(ss, 1)
1 in ss
2 in ss
push!(ss, 2)
2 in ss

delete!(ss, 1)
delete!(ss, 2)

SimpleSet([1,2,3,1,2,3])