text:
# This file contains functions used for the production of the plots in:
# Figure 6 of the main manuscript, as well as Supplementary Figures 4, 5, 6, and 7.
# For plotting we relied on the PGFPlotsX julia package, that essentially creates
# pgfplots code that is then compliled into the final figure.

# The parameters used for the specific plots in our manuscript are explained within the text.
# For simplicity, we have added a Jupyter Notebook that recreates the exact plots seen in the manuscript

using PGFPlotsX, Colors, LaTeXStrings

# Formatting instructions for the production of the final vector graphic using PGFPlotsX
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{cmbright}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage[T1]{fontenc}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\DeclareFontFamily{U}{futm}{}
    \DeclareFontShape{U}{futm}{m}{n}{
        <-> s * [.92] fourier-bb
    }{}
    \DeclareSymbolFont{Ufutm}{U}{futm}{m}{n}
    \DeclareSymbolFontAlphabet{\mathbb}{Ufutm}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\definecolor{colorBm}{RGB}{178,223,138}")
push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\definecolor{colorBp}{RGB}{31,120,180}")

function plotBmPop(xs, ys, error, title)
    ggplot2_plot_theme = @pgf {
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }
    mycolor = raw"colorBm"

    @pgf SemiLogXAxis(
        {
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys[1] .+ error)),
            title = title,
            xlabel = "Dissemination Time", #"Cell Population",
            ylabel = L"\large$\mathbb{E}[B_{md}]$",
        },
        Plot(
            {"name path=f", color = mycolor, ggplot2_plot_theme...},
            Coordinates(xs, ys[1]),
        ),
        Plot(
            {"name path=g", color = mycolor, opacity = 0.35, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .+ error),
        ),
        Plot(
            {"name path=q", color = mycolor, opacity = 0.35, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .- error),
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.35},
            raw"fill between [of=f and g]",
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.35},
            raw"fill between [of=f and q]",
        ),
    )
end

function plotBp(xs, ys, error, title)

    ggplot2_plot_theme = @pgf {
        #         mark="square",
        #         mark_size = 1,
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }

    mycolor = raw"colorBp"

    @pgf SemiLogXAxis(
        {
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys[1] .+ error)),
            title = title,
            xlabel = "Dissemination Time", #"Cell Population",
            ylabel = L"\large$\mathbb{E}[B_{p}]$",
        },
        # [
        #     PlotInc(
        #         {ggplot2_plot_theme..., color = colors[i]},
        #         Coordinates(xs, _y; yerror = error),
        #     ) for (i, _y) in enumerate(ys)
        # ]...,
        Plot(
            {"name path=f", color = mycolor, ggplot2_plot_theme...},
            Coordinates(xs, ys[1]),
        ),
        Plot(
            {"name path=g", color = mycolor, opacity = 0.35, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .+ error),
        ),
        Plot(
            {"name path=q", color = mycolor, opacity = 0.35, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .- error),
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.35},
            raw"fill between [of=f and g]",
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.35},
            raw"fill between [of=f and q]",
        ),
    )
end

function BmMutIndexPlotPrepare(detectionFrequency, MAX_MUT, mutationProbability, ρ)
    k_max = MAX_MUT
    y = [Bm(detectionFrequency, mutationProbability, ρ, k_seed) for k_seed = 1:k_max]
    ys = [y]
    xs = 1:k_max
    error = sqrt.(BmVar.(detectionFrequency, mutationProbability, ρ, 1:k_max))
    options = Dict("u" => mutationProbability, "rho" => ρ, "f" => detectionFrequency)
    return xs, ys, error, options
end

function plotBmVar(xs, ys, options::Dict{String,Float64})

    ggplot2_plot_theme =
        @pgf {mark = "square", mark_size = 0, mark_options = "solid", line_width = "1.64pt"}

    n = length(ys)

    # Evenly spread out colors
     colors = range(HSL(colorant"red"), stop=HSL(colorant"green"), length=n)
    @pgf SemiLogXAxis(
        {
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = minimum(minimum(ys)),
            ymax = 20,
            title = L"Neutral Growth $(u = $"*
                   "$(options["u"]),"*
                    L"$f = $"*
                    "$(options["f"]))",
            xlabel = "Dissemination Time",
            ylabel = L"$\sqrt{\mathrm{Var}(B_{md}^{neutral})}$",
            legend_pos = "north west",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
        Legend([L"$\rho = 0.72$", L"$\rho = 0.81$", L"$\rho = 0.90$", L"$\rho = 0.99$"]),
    )
end

# This is the old one with the error bars
function plotBmPop_(xs, ys, error, title)
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        xmajorgrids,
        x_grid_style = "white",
        ymajorgrids,
        y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme = @pgf {
        #         mark="square",
        #         mark_size = 1,
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }

    n = length(ys)

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    @pgf SemiLogXAxis(
        #@pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys)),
            title = title,
            xlabel = "Cell Population",
            ylabel = L"$B_m$",
        },
        [
            PlotInc(
                {ggplot2_plot_theme..., color = colors[i]},
                Coordinates(xs, _y; yerror = error),
            ) for (i, _y) in enumerate(ys)
        ]...,
    )
end

function plotBmPop(xs, ys, error, title)
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        # xmajorgrids,
        # x_grid_style = "white",
        # ymajorgrids,
        # y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme = @pgf {
        #         mark="square",
        #         mark_size = 1,
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }

    n = length(ys)
    n = 3
    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    mycolor = colors[2]
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{cmbright}")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage[T1]{fontenc}")

    @pgf SemiLogXAxis(
        #@pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys[1] .+ error)),
            title = title,
            xlabel = "Seeding Time", #"Cell Population",
            ylabel = L"$B_{md}$",
        },
        # [
        #     PlotInc(
        #         {ggplot2_plot_theme..., color = colors[i]},
        #         Coordinates(xs, _y; yerror = error),
        #     ) for (i, _y) in enumerate(ys)
        # ]...,
        Plot(
            {"name path=f", color = mycolor, ggplot2_plot_theme...},
            Coordinates(xs, ys[1]),
        ),
        Plot(
            {"name path=g", color = mycolor, opacity = 0.25, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .+ error),
        ),
        Plot(
            {"name path=q", color = mycolor, opacity = 0.25, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .- error),
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.25},
            raw"fill between [of=f and g]",
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.25},
            raw"fill between [of=f and q]",
        ),
    )
end

function plotBm(xs, ys, error, title)
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        xmajorgrids,
        x_grid_style = "white",
        ymajorgrids,
        y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme = @pgf {
        #         mark="square",
        #         mark_size = 1,
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }

    n = length(ys)

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    @pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys)),
            title = title,
            xlabel = "Time",
            ylabel = L"$B_{md}$",
        },
        [
            PlotInc(
                {ggplot2_plot_theme..., color = colors[i]},
                Coordinates(xs, _y; yerror = error),
            ) for (i, _y) in enumerate(ys)
        ]...,
    )
end

function plotBmNoVar(xs, ys, options::Dict{String,Float64})
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        xmajorgrids,
        x_grid_style = "white",
        ymajorgrids,
        y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme = @pgf {mark = "square", mark_size = 1, line_width = "1.64pt"}

    n = length(ys)

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    @pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys)),
            title = options["title"],
            xlabel = "Time",
            ylabel = L"$B_m$",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
    )
end

function plotPdetFreq(xs, ys, options::Dict{String,Float64})
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        xmajorgrids,
        x_grid_style = "white",
        ymajorgrids,
        y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme =
        @pgf {mark = "square", mark_size = 1, mark_options = "solid", line_width = "1.64pt"}

    n = length(ys)

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    @pgf SemiLogXAxis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys)),
            title = L"Neutral Growth $(u = $"*
                   "$(options["u"]),"*
                    L"$\rho = $"*
                    "$(options["rho"]))",
            xlabel = "Frequency",
            ylabel = L"$\phi(c_j,f,\rho)$",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
    )
end

function plotQuick(xs, ys, title)
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        xmajorgrids,
        x_grid_style = "white",
        ymajorgrids,
        y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme =
        @pgf {mark = "square", mark_size = 0.5, mark_options = "solid", line_width = "1pt"}

    n = length(ys)

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    @pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = minimum(minimum(ys)),
            ymax = 1,
            title = title,
            xlabel = "Seeding time",
            ylabel = L"$\Pr[\mathrm{sampling\ from\ type}_i](t)$",
            legend_pos = "north west",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
        Legend(["type-0", "type-1", "type-2"]),
    )
end

function plotDecays(xs, ys, title)
    ggplot2_axis_theme = @pgf {
        # tick_align = "outside",
        # tick_pos = "left",
        # xmajorgrids,
        # x_grid_style = "white",
        # ymajorgrids,
        # y_grid_style = "white",
        # axis_line_style = "white",
        # "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme =
        @pgf {mark = "square", mark_size = 0.5, mark_options = "solid", line_width = "1pt"}

    n = length(ys)

    # Evenly spread out colors
    colors = range(HSL(colorant"red"), stop=HSL(colorant"green"), length=n)
    @pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = minimum(minimum(ys)),
            ymax = 1,
            title = title,
            xlabel = L"$j$",
            ylabel = L"$d_{c_j}$",
            # legend_pos = " east",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
        Legend(["Exponential Decay", "Logarithmic Decay", "Polynomial Decay"]),
    )
end

function plotBmPMFs(xs, ys, title)
    ggplot2_axis_theme = @pgf {
        # tick_align = "outside",
        # tick_pos = "left",
        # xmajorgrids,
        # x_grid_style = "white",
        # ymajorgrids,
        # y_grid_style = "white",
        # axis_line_style = "white",
        # "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme =
        @pgf {mark = "square", mark_size = 0.5, mark_options = "solid", line_width = "1pt"}

    n = length(ys)

    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{cmbright}")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage[T1]{fontenc}")

    # Evenly spread out colors
    colors = range(HSL(colorant"red"), stop=HSL(colorant"green"), length=n)
    @pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = minimum(minimum(ys)),
            ymax = 1,
            title = title,
            xlabel = L"$j$",
            ylabel = L"$B_{md}$",
            legend_pos = "north east",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
        Legend(["Exponential Decay", "Logarithmic Decay", "Polynomial Decay"]),
    )
end

function BmPopPlotPrepare(detectionFrequency, detectionPopulation, mutationProbability, ρ)
    k_max = kmaxEstimate(detectionPopulation, mutationProbability, ρ)
    y = [Bm(detectionFrequency, mutationProbability, ρ, k_seed) for k_seed = 1:k_max]
    ys = [y]
    xs = [popEstimate(s, mutationProbability, ρ) for s = 1:k_max]
    error = sqrt.(BmVar.(detectionFrequency, mutationProbability, ρ, 1:k_max))
    options = Dict("u" => mutationProbability, "rho" => ρ, "f" => detectionFrequency)
    return xs, ys, error, options
end
function BmMutIndexPlotPrepare(detectionFrequency, MAX_MUT, mutationProbability, ρ)
    k_max = MAX_MUT
    y = [Bm(detectionFrequency, mutationProbability, ρ, k_seed) for k_seed = 1:k_max]
    ys = [y]
    xs = 1:k_max
    error = sqrt.(BmVar.(detectionFrequency, mutationProbability, ρ, 1:k_max))
    options = Dict("u" => mutationProbability, "rho" => ρ, "f" => detectionFrequency)
    return xs, ys, error, options
end
function BpMutIndexPlotPrepare(MAX_MUT, Bm_s, McMs_sum)
    k_max = MAX_MUT
#     y = [max(McMs_sum - k_seed + Bm_s[k_seed], 0) for k_seed = 1:k_max]
    y = [McMs_sum - k_seed + Bm_s[k_seed] for k_seed = 1:k_max]
    ys = [y]
    xs = 1:k_max
    return xs, ys
end

function plotBpUpperLower(xs, ys, title)
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        xmajorgrids,
        x_grid_style = "white",
        ymajorgrids,
        y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme = @pgf {
        #         mark="square",
        #         mark_size = 1,
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }

    n = length(ys)

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    @pgf SemiLogXAxis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            #              ymin = minimum(ys),   ymax = 1.05*maximum(maximum(ys)),
            title = title,
            xlabel = "Seeding Time",
            ylabel = L"$B_p$",
        },
        [
            PlotInc({ggplot2_plot_theme..., color = colors[i]}, Coordinates(xs, _y))
            for (i, _y) in enumerate(ys)
        ]...,
    )
end

function plotBp(xs, ys, error, title)
    ggplot2_axis_theme = @pgf {
        tick_align = "outside",
        tick_pos = "left",
        # xmajorgrids,
        # x_grid_style = "white",
        # ymajorgrids,
        # y_grid_style = "white",
        axis_line_style = "white",
        "axis_background/.style" = {fill = "white!89.803921568627459!black"},
    }

    ggplot2_plot_theme = @pgf {
        #         mark="square",
        #         mark_size = 1,
        "no marks",
        "error bars/y dir=both",
        "error bars/y explicit",
        line_width = "1.64pt",
    }

    n = length(ys)
    n = 3

    # Evenly spread out colors
    colors = [LCHuv(65, 100, h) for h in range(15; stop = 360 + 15, length = n + 1)][1:n]
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepgfplotslibrary{fillbetween}")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage{cmbright}")
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usepackage[T1]{fontenc}")
    mycolor = colors[3]
    @pgf SemiLogXAxis(
        #@pgf Axis(
        {
            ggplot2_axis_theme...,
            xmin = minimum(xs),
            xmax = 1.05 * maximum(xs),
            ymin = 0,
            ymax = 1.05 * maximum(maximum(ys[1] .+ error)),
            title = title,
            xlabel = "Seeding Time", #"Cell Population",
            ylabel = L"$B_{p}$",
        },
        # [
        #     PlotInc(
        #         {ggplot2_plot_theme..., color = colors[i]},
        #         Coordinates(xs, _y; yerror = error),
        #     ) for (i, _y) in enumerate(ys)
        # ]...,
        Plot(
            {"name path=f", color = mycolor, ggplot2_plot_theme...},
            Coordinates(xs, ys[1]),
        ),
        Plot(
            {"name path=g", color = mycolor, opacity = 0.25, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .+ error),
        ),
        Plot(
            {"name path=q", color = mycolor, opacity = 0.25, ggplot2_plot_theme...},
            Coordinates(xs, ys[1] .- error),
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.25},
            raw"fill between [of=f and g]",
        ),
        Plot(
            {thick, color = mycolor, fill = mycolor, opacity = 0.25},
            raw"fill between [of=f and q]",
        ),
    )
end

max_stars_repo_name: SunPathLab/Comet

max_issues_repo_name: SunPathLab/Comet

max_forks_repo_name: SunPathLab/Comet


text:
module SMCExamples

using SequentialMonteCarlo
using RNGPool
using StaticArrays

include("particles.jl")
include("lgModel.jl")
include("mvlgModel.jl")
include("smcSampler.jl")
include("finiteFK.jl")
include("nettoModel.jl")
include("lorenz96Model.jl")

end # module

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# L is true if the view itself supports fast linear indexing
struct SubQuasiArray{T,N,P,I,L} <: AbstractQuasiArray{T,N}
    parent::P
    indices::I
    offset1::Int       # for linear indexing and pointer, only valid when L==true
    stride1::Int       # used only for linear indexing
    function SubQuasiArray{T,N,P,I,L}(parent, indices, offset1, stride1) where {T,N,P,I,L}
        @_inline_meta
        check_parent_index_match(parent, indices)
        new(parent, indices, offset1, stride1)
    end
end
# Compute the linear indexability of the indices, and combine it with the linear indexing of the parent
function SubQuasiArray(parent::AbstractQuasiArray, indices::Tuple)
    @_inline_meta
    SubQuasiArray(IndexStyle(viewindexing(indices), IndexStyle(parent)), parent, ensure_indexable(indices), index_dimsum(indices...))
end
function SubQuasiArray(::IndexCartesian, parent::P, indices::I, ::NTuple{N,Any}) where {P,I,N}
    @_inline_meta
    SubQuasiArray{eltype(P), N, P, I, false}(parent, indices, 0, 0)
end
function SubQuasiArray(::IndexLinear, parent::P, indices::I, ::NTuple{N,Any}) where {P,I,N}
    @_inline_meta
    # Compute the stride and offset
    stride1 = compute_stride1(parent, indices)
    SubQuasiArray{eltype(P), N, P, I, true}(parent, indices, compute_offset1(parent, stride1, indices), stride1)
end

check_parent_index_match(parent::AbstractQuasiArray{T,N}, ::NTuple{N, Bool}) where {T,N} = nothing
viewindexing(I::Tuple{AbstractQuasiArray, Vararg{Any}}) = IndexCartesian()

# Simple utilities
size(V::SubQuasiArray) = (@_inline_meta; map(n->cardinality(n), axes(V)))

similar(V::SubQuasiArray, T::Type, dims::Dims) = similar(V.parent, T, dims)

sizeof(V::SubQuasiArray) = length(V) * sizeof(eltype(V))

parent(V::SubQuasiArray) = V.parent
parentindices(V::SubQuasiArray) = V.indices
parentindices(a::AbstractQuasiArray) = map(OneTo, size(a))

## Aliasing detection
dataids(A::SubQuasiArray) = (dataids(A.parent)..., _splatmap(dataids, A.indices)...)
unaliascopy(A::SubQuasiArray) = typeof(A)(unaliascopy(A.parent), map(unaliascopy, A.indices), A.offset1, A.stride1)

# Transform indices to be "dense"
_trimmedindex(i::AbstractQuasiArray) = oftype(i, reshape(eachindex(IndexLinear(), i), axes(i)))

## SubArray creation
# We always assume that the dimensionality of the parent matches the number of
# indices that end up getting passed to it, so we store the parent as a
# ReshapedArray view if necessary. The trouble is that arrays of `CartesianIndex`
# can make the number of effective indices not equal to length(I).
_maybe_reshape_parent(A::AbstractQuasiArray, ::NTuple{1, Bool}) = reshape(A, Val(1))
_maybe_reshape_parent(A::AbstractQuasiArray{<:Any,1}, ::NTuple{1, Bool}) = reshape(A, Val(1))
_maybe_reshape_parent(A::AbstractQuasiArray{<:Any,N}, ::NTuple{N, Bool}) where {N} = A
_maybe_reshape_parent(A::AbstractQuasiArray, ::NTuple{N, Bool}) where {N} = reshape(A, Val(N))

function view(A::AbstractQuasiArray, I::Vararg{Any,N}) where {N}
    @_inline_meta
    J = map(i->unalias(A,i), to_indices(A, I))
    @boundscheck checkbounds(A, J...)
    unsafe_view(_maybe_reshape_parent(A, index_ndims(J...)), J...)
end

const ViewIndex{IND} = Union{IND,AbstractArray{<:IND}}
const QViewIndex{IND} = Union{IND,AbstractArray{<:IND},AbstractQuasiArray{<:IND}}

# This computes the linear indexing compatibility for a given tuple of indices
quasi_viewindexing(::Tuple{}, I::Tuple{}) = IndexLinear()
# Leading scalar indices simply increase the stride
quasi_viewindexing(axs::Tuple{AbstractQuasiVector{IND}, Vararg{Any}}, I::Tuple{IND, Vararg{Any}}) where IND =
    (@_inline_meta; quasi_viewindexing(tail(axs), tail(I)))
quasi_viewindexing(axs::Tuple{AbstractVector{IND}, Vararg{Any}}, I::Tuple{IND, Vararg{Any}}) where IND =
    (@_inline_meta; quasi_viewindexing(tail(axs), tail(I)))

# Slices may begin a section which may be followed by any number of Slices
# quasi_viewindexing(axs, I::Tuple{Slice, Slice, Vararg{Any}}) = (@_inline_meta; quasi_viewindexing(tail(I)))
# # A UnitRange can follow Slices, but only if all other indices are scalar
# quasi_viewindexing(I::Tuple{Slice, AbstractUnitRange, Vararg{ScalarIndex}}) = IndexLinear()
# quasi_viewindexing(I::Tuple{Slice, Slice, Vararg{ScalarIndex}}) = IndexLinear() # disambiguate
# # In general, ranges are only fast if all other indices are scalar
# quasi_viewindexing(I::Tuple{AbstractRange, Vararg{ScalarIndex}}) = IndexLinear()
# # All other index combinations are slow
# quasi_viewindexing(I::Tuple{Vararg{Any}}) = IndexCartesian()
# # Of course, all other array types are slow
quasi_viewindexing(axs::Tuple{AbstractQuasiVector{IND}, Vararg{Any}}, I::Tuple{AbstractArray{IND}, Vararg{Any}}) where IND = IndexCartesian()
quasi_viewindexing(axs::Tuple{AbstractVector{IND}, Vararg{Any}}, I::Tuple{AbstractArray{IND}, Vararg{Any}}) where IND = IndexCartesian()

# combined dimensionality of all indices
# rather than returning N, it returns an NTuple{N,Bool} so the result is inferrable
@inline quasi_index_dimsum(axs::Tuple{AbstractQuasiVector{IND},Vararg{Any}}, inds::Tuple{IND,Vararg{Any}}) where IND =
    (quasi_index_dimsum(tail(axs), tail(inds))...,)
@inline quasi_index_dimsum(axs::Tuple{AbstractVector{IND},Vararg{Any}}, inds::Tuple{IND,Vararg{Any}}) where IND =
    (quasi_index_dimsum(tail(axs), tail(inds))...,)
@inline quasi_index_dimsum(axs, inds::Tuple{Colon,Vararg{Any}}) = (true, quasi_index_dimsum(tail(axs), tail(inds))...)
@inline quasi_index_dimsum(axs::Tuple{AbstractQuasiVector{IND},Vararg{Any}}, inds::Tuple{AbstractArray{IND,N},Vararg{Any}}) where {N,IND} =
    (ntuple(x->true, Val(N))..., quasi_index_dimsum(tail(axs), tail(inds))...)
@inline quasi_index_dimsum(axs::Tuple{AbstractVector{IND},Vararg{Any}}, inds::Tuple{AbstractArray{IND,N},Vararg{Any}}) where {N,IND} =
    (ntuple(x->true, Val(N))..., quasi_index_dimsum(tail(axs), tail(inds))...)
quasi_index_dimsum(::Tuple{}, ::Tuple{}) = ()

function SubArray(parent::AbstractQuasiArray, indices::Tuple)
    @_inline_meta
    SubArray(IndexStyle(quasi_viewindexing(axes(parent), indices), IndexStyle(parent)), parent, ensure_indexable(indices), quasi_index_dimsum(axes(parent), indices))
end

unsafe_view(A::AbstractQuasiArray, I...) = _unsafe_view(indextype(A), A, Base.to_indices(A, I))
unsafe_view(V::SubQuasiArray, I...) = (@_inline_meta; _maybe_reindex(V, Base.to_indices(V, I)))
_unsafe_view(::Type{INDS}, A::AbstractQuasiArray, I::Tuple) where {INDS} = _unsafe_view_type(INDS, I)(A, I)
# If Array-based indexing, use SubArray, otherwise SubQuasiArray
_unsafe_view_type(::Type{Tuple{}}, ::Tuple) = SubArray
_unsafe_view_type(inds::Type{<:Tuple{IND,Vararg{Any}}}, I::Tuple{ViewIndex{IND},Vararg{Any}}) where IND = _unsafe_view_type(Base.tuple_type_tail(inds), tail(I))
_unsafe_view_type(inds::Type{<:Tuple{IND,Vararg{Any}}}, I::Tuple{AbstractQuasiArray{IND},Vararg{Any}}) where IND = SubQuasiArray
# _unsafe_view_type(inds::Type{<:Tuple{IND,Vararg{Any}}}, I::Tuple{Any,Vararg{Any}}) where IND = _unsafe_view_type(inds, tuple(convert(IND, I[1]), tail(I)...))

# When we take the view of a view, it's often possible to "reindex" the parent
# view's indices such that we can "pop" the parent view and keep just one layer
# of indirection. But we can't always do this because arrays of `CartesianIndex`
# might span multiple parent indices, making the reindex calculation very hard.
# So we use _maybe_reindex to figure out if there are any arrays of
# `CartesianIndex`, and if so, we punt and keep two layers of indirection.
_unsafe_view(::Type{IND}, V::SubQuasiArray, I::Vararg{ViewIndex{IND},N}) where {N,IND} =
    (@_inline_meta; _maybe_reindex(V, I))
_unsafe_view(::Type{IND}, V::SubQuasiArray, I::Vararg{QViewIndex{IND},N}) where {N,IND} =
    (@_inline_meta; _maybe_reindex(V, I))

_maybe_reindex(V, I) = (@_inline_meta; _maybe_reindex(V, I, I))
# _maybe_reindex(V, I, ::Tuple{AbstractArray{<:AbstractCartesianIndex}, Vararg{Any}}) =
#     (@_inline_meta; SubArray(V, I))
# # But allow arrays of CartesianIndex{1}; they behave just like arrays of Ints
# _maybe_reindex(V, I, A::Tuple{AbstractArray{<:AbstractCartesianIndex{1}}, Vararg{Any}}) =
#     (@_inline_meta; _maybe_reindex(V, I, tail(A)))
_maybe_reindex(V, I, A::Tuple{Any, Vararg{Any}}) = (@_inline_meta; _maybe_reindex(V, I, tail(A)))

_subarray(::Type{IND}, A::AbstractArray, idxs) where IND = SubArray(A, idxs)
_subarray(::Type{IND}, A::AbstractQuasiArray, idxs) where IND = _unsafe_view_type(IND, idxs)(A, idxs)

function _maybe_reindex(V, I, ::Tuple{})
    @_inline_meta
    @inbounds idxs = to_indices(V.parent, quasi_reindex(axes(parent(V)), V.indices, I))
    _subarray(indextype(V.parent), V.parent, idxs)
end

Base._maybe_reindex(V::SubArray{<:Any,N,<:AbstractQuasiArray}, I, ::Tuple{}) where N = _maybe_reindex(V, I, ())

## Re-indexing is the heart of a view, transforming A[i, j][x, y] to A[i[x], j[y]]
#
# Recursively look through the heads of the parent- and sub-indices, considering
# the following cases:
# * Parent index is array  -> re-index that with one or more sub-indices (one per dimension)
# * Parent index is Colon  -> just use the sub-index as provided
# * Parent index is scalar -> that dimension was dropped, so skip the sub-index and use the index as is

AbstractZeroDimQuasiArray{T} = AbstractQuasiArray{T, 0}

# Re-index into parent vectors with one subindex
reindex(idxs::Tuple{AbstractQuasiVector, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1]], reindex(tail(idxs), tail(subidxs))...))

# Parent matrices are re-indexed with two sub-indices
reindex(idxs::Tuple{AbstractQuasiMatrix, Vararg{Any}}, subidxs::Tuple{Any, Any, Vararg{Any}}) =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1], subidxs[2]], reindex(tail(idxs), tail(tail(subidxs)))...))

# In general, we index N-dimensional parent arrays with N indices
@generated function reindex(idxs::Tuple{AbstractQuasiArray{T,N}, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) where {T,N}
    if length(subidxs.parameters) >= N
        subs = [:(subidxs[$d]) for d in 1:N]
        tail = [:(subidxs[$d]) for d in N+1:length(subidxs.parameters)]
        :(@_propagate_inbounds_meta; (idxs[1][$(subs...)], reindex(tail(idxs), ($(tail...),))...))
    else
        :(throw(ArgumentError("cannot re-index SubArray with fewer indices than dimensions\nThis should not occur; please submit a bug report.")))
    end
end

# indices are taken from the range/vector
# Since bounds-checking is performance-critical and uses
# indices, it's worth optimizing these implementations thoroughly
axes(S::SubArray{T,N,<:AbstractQuasiArray}) where {T,N} =
    _quasi_indices_sub(axes(parent(S)), S.indices)
_quasi_indices_sub(axs::Tuple{AbstractQuasiVector{IND},Vararg{Any}}, inds::Tuple{IND,Vararg{Any}}) where IND =
    (@_inline_meta; _quasi_indices_sub(tail(axs), tail(inds)))
_quasi_indices_sub(axs::Tuple{AbstractVector{IND},Vararg{Any}}, inds::Tuple{IND,Vararg{Any}}) where IND =
    (@_inline_meta; _quasi_indices_sub(tail(axs), tail(inds)))
_quasi_indices_sub(::Tuple{}, ::Tuple{}) = ()
function _quasi_indices_sub(axs::Tuple{AbstractQuasiVector{IND},Vararg{Any}}, inds::Tuple{AbstractArray{IND},Vararg{Any}}) where IND
    @_inline_meta
    (unsafe_indices(inds[1])..., _quasi_indices_sub(tail(axs), tail(inds))...)
end
function _quasi_indices_sub(axs::Tuple{AbstractVector{IND},Vararg{Any}}, inds::Tuple{AbstractArray{IND},Vararg{Any}}) where IND
    @_inline_meta
    (unsafe_indices(inds[1])..., _quasi_indices_sub(tail(axs), tail(inds))...)
end

quasi_reindex(axs::Tuple{AbstractQuasiOrVector{IND}, Vararg{Any}}, idxs::Tuple{IND, Vararg{Any}}, subidxs::Tuple{Vararg{Any}}) where IND =
    (@_propagate_inbounds_meta; (idxs[1], quasi_reindex(tail(axs), tail(idxs), subidxs)...))
quasi_reindex(axs::Tuple{AbstractQuasiOrVector{IND}, Vararg{Any}}, idxs::Tuple{AbstractQuasiOrVector{IND}, Vararg{Any}}, subidxs::Tuple{Any, Vararg{Any}}) where IND =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1]], quasi_reindex(tail(axs), tail(idxs), tail(subidxs))...))

# Parent matrices are re-indexed with two sub-indices
quasi_reindex(axs::Tuple{AbstractQuasiOrVector{IND}, Vararg{Any}},idxs::Tuple{AbstractQuasiOrMatrix{IND}, Vararg{Any}}, subidxs::Tuple{Any, Any, Vararg{Any}}) where IND =
    (@_propagate_inbounds_meta; (idxs[1][subidxs[1], subidxs[2]], quasi_reindex(tail(axs), tail(idxs), tail(tail(subidxs)))...))

quasi_reindex(::Tuple{}, ::Tuple{}, ::Tuple{}) = ()

function getindex(V::SubArray{T,N,<:AbstractQuasiArray}, I::Vararg{Int,N}) where {T,N}
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[quasi_reindex(axes(parent(V)), V.indices, I)...]
    r
end

# In general, we simply re-index the parent indices by the provided ones
SlowSubQuasiArray{T,N,P,I} = SubQuasiArray{T,N,P,I,false}
function _getindex(::Type{IND}, V::SlowSubQuasiArray, I::IND) where IND
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds r = V.parent[quasi_reindex(axes(parent(V)), V.indices, I)...]
    r
end

function _setindex!(::Type{IND}, V::SlowSubQuasiArray, x, I::IND) where IND
    @_inline_meta
    @boundscheck checkbounds(V, I...)
    @inbounds V.parent[quasi_reindex(axes(parent(V)), V.indices, I)...] = x
    V
end

IndexStyle(::Type{<:SubQuasiArray}) = IndexCartesian()

# Strides are the distance in memory between adjacent elements in a given dimension
# which we determine from the strides of the parent
strides(V::SubQuasiArray) = substrides(V.parent, V.indices)
substrides(parent, strds, I::Tuple{Any, Vararg{Any}}) = throw(ArgumentError("strides is invalid for SubQuasiArrays with indices of type $(typeof(I[1]))"))

stride(V::SubQuasiArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]

compute_stride1(parent::AbstractQuasiArray, I::NTuple{N,Any}) where {N} =
    1

elsize(::Type{<:SubQuasiArray{<:Any,<:Any,P}}) where {P} = elsize(P)

function first_index(V::SubQuasiArray)
    P, I = parent(V), V.indices
    s1 = compute_stride1(P, I)
    s1 + compute_offset1(P, s1, I)
end

# Computing the first index simply steps through the indices, accumulating the
# sum of index each multiplied by the parent's stride.
# The running sum is `f`; the cumulative stride product is `s`.
# If the parent is a vector, then we offset the parent's own indices with parameters of I
compute_offset1(parent::AbstractQuasiVector, stride1::Integer, I::Tuple{AbstractRange}) = 0
# If the result is one-dimensional and it's a Colon, then linear
# indexing uses the indices along the given dimension. Otherwise
# linear indexing always starts with 1.

function compute_linindex(parent, I::NTuple{N,Any}) where N
    @_inline_meta
    IP = fill_to_length(axes(parent), OneTo(1), Val(N))
    compute_linindex(1, 1, IP, I)
end
function compute_linindex(f, s, IP::Tuple, I::Tuple{ScalarIndex, Vararg{Any}})
    @_inline_meta
    Δi = I[1]-first(IP[1])
    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
end
function compute_linindex(f, s, IP::Tuple, I::Tuple{Any, Vararg{Any}})
    @_inline_meta
    Δi = first(I[1])-first(IP[1])
    compute_linindex(f + Δi*s, s*unsafe_length(IP[1]), tail(IP), tail(I))
end
compute_linindex(f, s, IP::Tuple, I::Tuple{}) = f

find_extended_dims(dim, ::ScalarIndex, I...) = (@_inline_meta; find_extended_dims(dim + 1, I...))
find_extended_dims(dim, i1, I...) = (@_inline_meta; (dim, find_extended_dims(dim + 1, I...)...))
find_extended_dims(dim) = ()
find_extended_inds(::ScalarIndex, I...) = (@_inline_meta; find_extended_inds(I...))
find_extended_inds(i1, I...) = (@_inline_meta; (i1, find_extended_inds(I...)...))
find_extended_inds() = ()

unsafe_convert(::Type{Ptr{T}}, V::SubQuasiArray{T,N,P,<:Tuple{Vararg{RangeIndex}}}) where {T,N,P} =
    unsafe_convert(Ptr{T}, V.parent) + (first_index(V)-1)*sizeof(T)

pointer(V::SubQuasiArray, i::Int) = _pointer(V, i)
_pointer(V::SubQuasiArray{<:Any,1}, i::Int) = pointer(V, (i,))
_pointer(V::SubQuasiArray, i::Int) = pointer(V, Base._ind2sub(axes(V), i))

# indices are taken from the range/vector
# Since bounds-checking is performance-critical and uses
# indices, it's worth optimizing these implementations thoroughly
axes(S::SubQuasiArray) = (@_inline_meta; _indices_sub(S, S.indices...))
_indices_sub(S::SubQuasiArray) = ()
_indices_sub(S::SubQuasiArray, ::Number, I...) = (@_inline_meta; _indices_sub(S, I...))
function _indices_sub(S::SubQuasiArray, i1::Union{AbstractQuasiArray,AbstractArray}, I...)
    @_inline_meta
    (unsafe_indices(i1)..., _indices_sub(S, I...)...)
end

@propagate_inbounds maybeview(A::AbstractQuasiArray, args...) = view(A, args...)
@propagate_inbounds maybeview(A::AbstractQuasiArray, args::Number...) = getindex(A, args...)
@propagate_inbounds maybeview(A::AbstractQuasiArray) = getindex(A)

##
# MemoryLayout
##

@inline MemoryLayout(A::Type{<:SubQuasiArray{T,N,P,I}}) where {T,N,P,I} =
    sublayout(MemoryLayout(P), I)

@inline sub_materialize(_, V::AbstractQuasiArray, _) = QuasiArray(V)
@inline sub_materialize(V::SubQuasiArray) = sub_materialize(MemoryLayout(typeof(V)), V)
@inline sub_materialize(A::AbstractQuasiArray) = A # already materialized

###
# show
###

function summary(io::IO, v::SubQuasiArray)
    print(io, "view(")
    summary(io, parent(v))
    Base.showindices(io, v.indices...)
    print(io, ')')
    print(io, " with eltype ", eltype(v))
end


text:
## Общие функции
# Поменять местами строки
function swaprows(A::Array{Int64,2},row::Int64,row2::Int64)
   for j = 1:size(A)[end]
      r = A[row,j]
      A[row,j] = A[row2,j]
      A[row2,j] = r
   end
   return A
end

#Поменять местами столбцы
function swapcols(A::Array{Int64,2},col::Int64,col2::Int64)
   for i = 1:size(A)[1]
      c = A[i,col]
      A[i,col] = A[i,col2]
      A[i,col2] = c
   end
   return A
end

#Поиск элемента с указанным или отличным значением, возвращает координаты этого элемента
function search(sr::Int64,A::Array{Int64,2},startr::Int64,startc::Int64,same::Bool)
   for i = startr+1:size(A)[1], j = startc+1:size(A)[end]
      if (same == true)
         if (A[i,j] == sr)
            return (i,j)
         end
      else
         if (A[i,j] != sr)
            return (i,j)
         end
      end
   end
   return 0
end

## Транспонирование
function trans(A::Array{Int64,2})
   S = Array{Int64,2}(undef,size(A)[1],size(A)[end])
   for i=1:size(A)[end]
      for j=1:size(A)[1]
         S[j,i] = A[i,j]
      end
   end
   return S
end

## Умножение
function mul(A::Array{Int64,2},B::Array{Int64,2})
   @assert size(A)[end] == size(B)[1] "The number of columns of the first matrix should be equal to the number of rows of the second"
   S = Array{Int64,2}(undef,size(A)[1],size(B)[end])
   for i = 1:size(A)[1], j = 1:size(B)[end]
      S[i,j] = 0
      for k = 1:size(A)[end]
         S[i,j] += a[i,k] * b[k,j]
      end
   end
   return S
end

## Определитель
# Функция получения матрицы без i-ой строки и j-ого столбца
function dethelp(A::Array{Int64,2},str::Int64,stlb::Int64)
   S = Array{Int64,2}(undef,size(A)[1]-1,size(A)[end]-1)
   di=0
   for i = 1:size(A)[1]-1
      if (i == str)
         di = 1
      end
      dj = 0
      for j = 1:size(A)[1] - 1
         if (j == stlb)
            dj = 1
         end
          S[i,j] = A[i + di ,j + dj];
      end
   end
   return S
end

# Рекурсивная функция определителя
function det(A::Array{Int64,2})
   @assert size(A)[end] == size(A)[1] "gag"
   d = 0
   k = 1
   if (size(A)[end] == 1)
      d=A[1,1]
      return d
   end
   if (size(A)[end] == 2)
      d = A[1,1] * A[2,2] - A[2,1] * A[1,2]
      return d
   end
   if (size(A)[end] > 2)
      for i = 1:size(A)[1]
         S = dethelp(A,i,1)
         d = d + k * A[i,1] * det(S)
         k = -k
      end
   end
   return d
end

## След
function tr(A::Array{Int64,2})
   @assert size(A)[end] == size(A)[1] "Kv"
   t = 0
   for i = 1:size(A)[1]
      t += A[i,i]
   end
   return t
end

##Ранг
function rang(A::Array{Int64,2})
   iter = min(size(A)[end],size(A)[1])
   for i = 1:iter
      if (A[i,i] == 0)
         check = search(0,A,i,i,false)
         if (!(typeof((1,1))==typeof(check)))
            break
         end
         if (i != check[1])
            swaprows(A,i,check[1])
         end
         if (i != check[2])
            swapcols(A,i,check[2])
         end
      end
      tmp=A[i,i]
      for x = i:size(A)[end]
         A[i,x] = div(A[i,x],tmp)
      end
      for y = (i+1):size(A)[1]
         tmp = A[y,i]
         for x = i:size(A)[end]
            A[y,x] -= (A[i,x] * tmp)
         end
      end
      for x = (i+1):size(A)[end]
         tmp = A[i,x]
         for y = i:size(A)[1]
            A[y,x] -= (A[y,i] * tmp)
         end
      end
   end
   cnt = 0
   for i = 1:iter
      if (A[i,i] == 0)
         break
      else
         cnt += 1
      end
   end
   if (cnt == 0)
      cnt += 1
   end
   return cnt
end

##Сложение матриц
function plus(A::Array{Int64,2},B::Array{Int64,2})
   @assert (size(A)[end] == size(B)[end]) || (size(A)[1] == size(B)[1]) "Kv"
   S = Array{Int64,2}(undef,size(A)[1],size(A)[end])
   for i = 1:size(A)[1], j = 1:size(A)[end]
      S[i,j] = A[i,j] + B[i,j]
   end
   return S
end

##Разность матриц
function minus(A::Array{Int64,2},B::Array{Int64,2})
   @assert (size(A)[end] == size(B)[end]) || (size(A)[1] == size(B)[1]) "Kv"
   S = Array{Int64,2}(undef,size(A)[1],size(A)[end])
   for i = 1:size(A)[1], j = 1:size(A)[end]
      S[i,j] = A[i,j] - B[i,j]
   end
   return S
end

##Обращение
function obr(A::Array{Int64,2})
   @assert size(A)[end] == size(A)[1] "Kv"
   @assert det(A) != 0 "det = 0"
   B = A
   S = zeros(size(A)[1],size(A)[1])
   for i = 1:size(A)[1]
      S[i,i] = 1
   end
   for k = 1:size(A)[1]
      tmp = A[k,k]
      for j = 1:size(A)[1]
         A[k,j] = A[k,j]/tmp
         S[k,j] = S[k,j]/tmp
      end
      for i = k + 1:size(A)[1]
         tmp = A[i,k]
         for j = 1:size(A)[1]
            A[i,j] -= A[k,j] * tmp;
            S[i,j] -= S[k,j] * tmp;
         end
      end
   end
   k = size(A)[1] - 1
   while (k > 1)
      while (i >= 1)
         tmp = A[i,k]
         for j=1:size(A)[1]
            A[i,j] -= A[k,j] * tmp;
            S[i,j] -= S[k,j] * tmp;
            println("A[",i,j,"]=",A[i,j])
            println("S[",i,j,"]=",S[i,j])
         end
         i -= 1
      end
      k -= 1
   end
   return S
end
b = [2 4; 1 4]
print(obr(b))


text:
# # This is a title

# This is a description
# # [This is another title](@id custom_id)

# This isn't a description

x = 1+1


text:
using DelayEmbeddings, Random
export ShuffleDimensions

"""
    ShuffleDimensions()
Multidimensional surrogates of input *datasets* (`DelayEmbeddings.Dataset`, which are
also multidimensional) that have shuffled dimensions in each point.

These surrogates destroy the state space structure of the dataset and are thus
suited to distinguish deterministic datasets from high dimensional noise.
"""
struct ShuffleDimensions <: Surrogate end

function surrogenerator(x, sd::ShuffleDimensions)
    @assert x isa Dataset "input `x` must be `DelayEmbeddings.Dataset` for `ShuffleDimensions`"
    return SurrogateGenerator(sd, x, nothing)
end

function (sg::SurrogateGenerator{<:ShuffleDimensions})()
    data = copy(sg.x.data)
    for i in 1:length(data)
        @inbounds data[i] = shuffle(data[i])
    end
    return Dataset(data)
end


text:
using Dates

function _create_metadata_impl(
    item_size::Int32,
    item_name::AbstractString,
    fields::AbstractVector{Field},
    time_field_offsets::AbstractVector{Int32},
    content_description::AbstractString,
    name_values::AbstractVector{NameValue},
    epoch_utc::DateTime,
    tick_duration::Dates.FixedPeriod
)::TeaFileMetadata
    sections = [
        ItemSection(item_size, item_name, fields),
        ContentDescriptionSection(content_description),
        NameValueSection(name_values),
        TimeSection(tea_epoch(epoch_utc), ticks_per_day(tick_duration), time_field_offsets)
    ]

    item_start = minimum_item_start(sections)
    metadata = TeaFileMetadata(item_start, 0, sections)
    return metadata
end

"""
Create metadata given an item type, and other optional information.

# Arguments
- `item_type::Type`: A type for each item that will be stored. Must be a bits type.

# Keywords
- `time_fields::AbstractVector{Symbol}`: Fields in `item_type` that represent times. The
    first such field will be the primary time field, and must be non-decreasing in the
    file.
- `content_description::AbstractString`: A description of the file.
- `name_values::AbstractVector{NameValue}`: Additional name-value pairs to include in the
    metadata.
- `epoch_utc::DateTime`: The epoch, in UTC. Defaults to midnight 1970-01-01.
- `tick_duration::Dates.FixedPeriod`: The length of one tick of the clock.
"""
function create_metadata(
    item_type::Type;
    time_fields::AbstractVector{Symbol}=Symbol[],
    content_description::AbstractString="",
    name_values::AbstractVector{NameValue}=NameValue[],
    epoch_utc::DateTime=DateTime(1970, 1, 1),
    tick_duration::Dates.FixedPeriod=Millisecond(1)
)::TeaFileMetadata
    if !isbitstype(item_type)
        throw(ArgumentError("$item_type is not a bits type."))
    end

    fields = [
        Field(
            field_type_id(fieldtype(item_type, i)),
            fieldoffset(item_type, i),
            string(name)
        )
        for (i, name) in enumerate(fieldnames(item_type))
    ]

    time_field_offsets = Int32[
        fieldoffset(item_type, i)
        for time_field in time_fields
        for (i, name) in enumerate(fieldnames(item_type)) if name == time_field
    ]

    return _create_metadata_impl(
        Int32(sizeof(item_type)), string(item_type.name.name), fields, time_field_offsets,
        content_description, name_values, epoch_utc, tick_duration
    )
end

"""
Create metadata that derives whether or not a field represents a time from whether it is
a `DateTime`. Such a field will always be encoded as milliseconds since the Julia epoch,
as an Int64.
"""
function create_metadata_julia_time(
    item_type::Type;
    item_type_name::Union{Nothing, AbstractString}=nothing,
    content_description::AbstractString="",
    name_values::AbstractVector{NameValue}=NameValue[]
)::TeaFileMetadata
    if !isbitstype(item_type)
        throw(ArgumentError("$item_type is not a bits type."))
    end

    # The epoch and tick duration for Julia `DateTime` objects.
    epoch_utc = DateTime(0) - Millisecond(Dates.DATETIMEEPOCH)
    tick_duration = Millisecond(1)

    # Any `DateTime` field should be registered as an Int64, and also as a time field.
    fields = Field[]
    time_field_offsets = Int32[]

    for i in 1:fieldcount(item_type)
        offset = fieldoffset(item_type, i)
        T = fieldtype(item_type, i)
        if T == DateTime
            T = Int64
            push!(time_field_offsets, offset)
        end
        push!(fields, Field(field_type_id(T), offset, string(fieldname(item_type, i))))
    end

    # Use the name of the type specified if an item type name is not given.
    item_type_name = if isnothing(item_type_name)
        string(item_type.name.name)
    else
        item_type_name
    end

    return _create_metadata_impl(
        Int32(sizeof(item_type)), item_type_name, fields, time_field_offsets,
        content_description, name_values, epoch_utc, tick_duration
    )
end

# TODO Make this NOT fail if the field offsets don't match; we can test that separately
# depending on how we are going to use it.
"""
Create a `NamedTuple` type that is compatible with the specified metadata.

This may not be possible, due to field offsets, in which case we throw.
"""
function create_item_namedtuple(metadata::TeaFileMetadata)::Type{<:NamedTuple}
    time_section = get_section(metadata, TimeSection)
    # Create DateTime fields if we know they are binary compatible, otherwise don't do
    # anything special with time fields.
    create_datetime_fields = is_julia_time_compatible(time_section)

    item_section = get_section(metadata, ItemSection)

    field_names = Symbol[]
    field_types = Type[]

    for field in item_section.fields
        push!(field_names, Symbol(field.name))
        field_type_ = field_type(field)
        if create_datetime_fields && field_type_ == Int64
            if field.offset in time_section.time_field_offsets
                field_type_ = DateTime
            end
        end
        push!(field_types, field_type_)
    end

    candidate = NamedTuple{tuple(field_names...), Tuple{field_types...}}

    if !is_item_compatible(candidate, metadata)
        throw(ArgumentError(
            "Unable to create compatible NamedTuple; probably due to field offsets"
        ))
    end

    return candidate
end


text:
a=1

print(a)

# where are we installing the package
pyimport("sys").executable

#----------- Windows and Ubuntu----------------------
Conda.add("pytorch"; channel="pytorch")
Conda.add("torchvision"; channel="pytorch")
Conda.add("spyder")
Conda.add("scipy")
Conda.add("pandas")

pip.main(["install", "Quantlib-Python"])


text:
using Iris
using Test

@testset "Laplacians" begin
    @testset "1D" begin
        @testset "Cartesian" begin
            lat = Lattices.Cartesian(.1)
            lap = @test_nowarn Laplacian{Iris.Common.Laplacians.Symmetric}(lat,ones(ComplexF64,10),ones(ComplexF64,11))
            for i ∈ 1:10 @test lap.l0[i,i] ≈ -200 end
            for i ∈ 1:9 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 2:10 @test lap.l0[i,i-1] ≈ 100 end
        end
    end
    @testset "2D" begin
        @testset "Cartesian" begin
            lat = Cartesian(.1,.2)
            lap = @test_nowarn Laplacian{Iris.Common.Laplacians.Symmetric}(lat,ones(ComplexF64,10),ones(ComplexF64,11),ones(ComplexF64,11),ones(ComplexF64,12))
            for i ∈ 1:10*11 @test lap.l0[i,i] ≈ -250 end
            for i ∈ 1:9 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 1:9 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 11:19 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 11:19 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 21:29 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 21:29 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 31:39 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 31:39 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 41:49 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 41:49 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 51:59 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 51:59 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 61:69 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 61:69 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 71:79 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 71:79 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 81:89 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 81:89 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 91:99 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 91:99 @test lap.l0[i,i+10] ≈ 25 end
            for i ∈ 101:109 @test lap.l0[i,i+1] ≈ 100 end
            for i ∈ 2:10 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 12:20 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 12:20 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 22:30 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 22:30 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 32:40 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 32:40 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 42:50 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 42:50 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 52:60 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 52:60 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 62:70 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 62:70 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 72:80 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 72:80 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 82:90 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 82:90 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 92:100 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 92:100 @test lap.l0[i,i-10] ≈ 25 end
            for i ∈ 102:110 @test lap.l0[i,i-1] ≈ 100 end
            for i ∈ 102:110 @test lap.l0[i,i-10] ≈ 25 end
            lap.l0[1,1] = 1im
            @test_nowarn conj!(lap)
            @test lap.l0[1,1] ≈ -1im
        end
    end
end


text:
module CSTParser
global debug = true

using Tokenize
import Base: length, first, last, getindex, setindex!
import Tokenize.Tokens
import Tokenize.Tokens: RawToken, AbstractToken, iskeyword, isliteral, isoperator, untokenize
import Tokenize.Lexers: Lexer, peekchar, iswhitespace

export ParseState, parse_expression

include("lexer.jl")
include("spec.jl")
include("utils.jl")
include("recovery.jl")
include("components/internals.jl")
include("components/keywords.jl")
include("components/lists.jl")
include("components/operators.jl")
include("components/strings.jl")
include("conversion.jl")
include("display.jl")
include("interface.jl")
include("iterate.jl")

"""
    parse_expression(ps)

Parses an expression until `closer(ps) == true`. Expects to enter the
`ParseState` the token before the the beginning of the expression and ends
on the last token.

Acceptable starting tokens are:
+ A keyword
+ An opening parentheses or brace.
+ An operator.
+ An instance (e.g. identifier, number, etc.)
+ An `@`.

"""
function parse_expression(ps::ParseState, esc_on_error = false)
    if kindof(ps.nt) === Tokens.ENDMARKER
        ret = mErrorToken(ps, UnexpectedToken)
    elseif (esc_on_error && ps.nt.kind == Tokens.ERROR)
        ret = EXPR(:errortoken, 0, 0)
    elseif kindof(ps.nt) ∈ term_c && !(kindof(ps.nt) === Tokens.END && ps.closer.square)
        ret = mErrorToken(ps, EXPR(next(ps)), UnexpectedToken)
    else
        next(ps)
        if iskeyword(kindof(ps.t)) && kindof(ps.t) != Tokens.DO
            ret = parse_kw(ps)
        elseif kindof(ps.t) === Tokens.LPAREN
            ret = parse_paren(ps)
        elseif kindof(ps.t) === Tokens.LSQUARE
            ret = @default ps parse_array(ps)
        elseif kindof(ps.t) === Tokens.LBRACE
            ret = @default ps @closebrace ps parse_braces(ps)
        elseif isinstance(ps.t) || isoperator(ps.t)
            if both_symbol_and_op(ps.t)
                ret = EXPR(:IDENTIFIER, ps)
            else
                ret = INSTANCE(ps)
            end
            if is_colon(ret) && !(iscomma(ps.nt) || kindof(ps.ws) == SemiColonWS)
                ret = parse_unary(ps, ret)
            elseif isoperator(ret) && assign_prec(valof(ret)) && !isunaryop(ret)
                ret = mErrorToken(ps, ret, UnexpectedAssignmentOp)
            end
        elseif kindof(ps.t) === Tokens.AT_SIGN
            ret = parse_macrocall(ps)
        else
            ret = mErrorToken(ps, INSTANCE(ps), UnexpectedToken)
        end
        ret = parse_compound_recur(ps, ret)
    end
    return ret
end

function parse_compound_recur(ps, ret)
    !closer(ps) ? parse_compound_recur(ps, parse_compound(ps, ret)) : ret
end

"""
    parse_compound(ps::ParseState, ret::EXPR)

Attempts to parse a compound expression given the preceding expression `ret`.
"""
function parse_compound(ps::ParseState, ret::EXPR)
    if kindof(ps.nt) === Tokens.FOR
        ret = parse_generator(ps, ret)
    elseif kindof(ps.nt) === Tokens.DO
        ret = @default ps @closer ps :block parse_do(ps, ret)
    elseif isajuxtaposition(ps, ret)
        if disallowednumberjuxt(ret)
            ret = mErrorToken(ps, ret, CannotJuxtapose)
        end
        ret = parse_operator(ps, ret, EXPR(:OPERATOR, 0, 0, "*"))
    elseif issuffixableliteral(ps, ret)
        arg = EXPR(:IDENTIFIER, next(ps))
        push!(ret, EXPR(:STRING, arg.fullspan, arg.span, val(ps.t, ps)))
    elseif (isidentifier(ret) || is_getfield(ret)) && isemptyws(ps.ws) && isprefixableliteral(ps.nt)
        ret = parse_prefixed_string_cmd(ps, ret)
    elseif kindof(ps.nt) === Tokens.LPAREN
        no_ws = !isemptyws(ps.ws)
        ret = @closeparen ps parse_call(ps, ret)
        if no_ws && !isunarycall(ret)
            ret = mErrorToken(ps, ret, UnexpectedWhiteSpace)
        end
    elseif kindof(ps.nt) === Tokens.LBRACE
        if isemptyws(ps.ws)
            ret = @default ps @nocloser ps :inwhere @closebrace ps parse_curly(ps, ret)
        else
            ret = mErrorToken(ps, (@default ps @nocloser ps :inwhere @closebrace ps parse_curly(ps, ret)), UnexpectedWhiteSpace)
        end
    elseif kindof(ps.nt) === Tokens.LSQUARE && isemptyws(ps.ws) && !isoperator(ret)
        ret = @default ps @nocloser ps :block parse_ref(ps, ret)
    elseif iscomma(ps.nt)
        ret = parse_tuple(ps, ret)
    elseif isunaryop(ret) && kindof(ps.nt) != Tokens.EQ
        ret = parse_unary(ps, ret)
    elseif isoperator(ps.nt)
        op = EXPR(:OPERATOR, next(ps))
        ret = parse_operator(ps, ret, op)
    elseif is_prime(ret.head)
        # prime operator followed by an identifier has an implicit multiplication
        nextarg = @precedence ps TimesOp parse_expression(ps)
        ret = EXPR(:call, EXPR[EXPR(:OPERATOR, 0, 0, "*"), ret, nextarg], nothing)
# ###############################################################################
# Everything below here is an error
# ###############################################################################
    else
        ps.errored = true
        if kindof(ps.nt) in (Tokens.RPAREN, Tokens.RSQUARE, Tokens.RBRACE)
            nextarg = mErrorToken(ps, EXPR(next(ps)), Unknown)
        else
            nextarg = try
                parse_expression(ps)
            catch err
                if err isa StackOverflowError
                    throw(error(string(ps, "\nsize: ", ps.l.io.size)))
                end
            end
        end
        ret = EXPR(:errortoken, EXPR[ret, nextarg], nothing)
    end
    return ret
end

"""
    parse_paren(ps, ret)

Parses an expression starting with a `(`.
"""
function parse_paren(ps::ParseState)
    args = EXPR[]
    trivia = EXPR[EXPR(ps)]
    @closeparen ps @default ps @nocloser ps :inwhere parse_comma_sep(ps, args, trivia, false, true, true, insert_params_at = 1)
    if length(args) == 1 && length(trivia) == 1 && ((kindof(ps.ws) !== SemiColonWS || headof(args[1]) === :block) && headof(args[1]) !== :parameters)
        accept_rparen(ps, trivia)
        ret = EXPR(:brackets, args, trivia)
    elseif VERSION < v"1.5" && length(args) == 1 && args[1].head === :parameters && isempty(args[1].args)
        accept_rparen(ps, trivia)
        pop!(args)
        push!(args, EXPR(:block, EXPR[], nothing))
        ret = EXPR(:brackets, args, trivia)
    else
        accept_rparen(ps, trivia)
        ret = EXPR(:tuple, args, trivia)
    end
    return ret
end

"""
    parse(str, cont = false)

Parses the passed string. If `cont` is true then will continue parsing until the end of the string returning the resulting expressions in a TOPLEVEL block.
"""
function parse(str::String, cont=false)
    ps = ParseState(str)
    x, _ = parse(ps, cont)
    return x
end

"""
    parse_doc(ps::ParseState)

Used for top-level parsing - attaches documentation (such as this) to expressions.
"""
function parse_doc(ps::ParseState)
    if (kindof(ps.nt) === Tokens.STRING || kindof(ps.nt) === Tokens.TRIPLE_STRING) && !isemptyws(ps.nws)
        doc = mLITERAL(next(ps))
        if kindof(ps.nt) === Tokens.ENDMARKER || kindof(ps.nt) === Tokens.END || ps.t.endpos[1] + 1 < ps.nt.startpos[1]
            ret = doc
        elseif isbinaryop(ps.nt) && !closer(ps)
            ret = parse_compound_recur(ps, doc)
        else
            ret = parse_expression(ps)
            ret = EXPR(:macrocall, EXPR[EXPR(:globalrefdoc, 0, 0), EXPR(:NOTHING, 0, 0), doc, ret], nothing)
        end
    else
        ret = parse_expression(ps)
    end
    if _continue_doc_parse(ps, ret)
        push!(ret, parse_expression(ps))
    end
    return ret
end

function parse(ps::ParseState, cont=false)
    if ps.l.io.size == 0
        return (cont ? EXPR(:file, EXPR[]) : nothing), ps
    end
    last_line = 0
    curr_line = 0

    if cont
        top = EXPR(:file, EXPR[], nothing)
        if kindof(ps.nt) === Tokens.WHITESPACE || kindof(ps.nt) === Tokens.COMMENT
            next(ps)
            push!(top, EXPR(:NOTHING, ps.nt.startbyte, ps.nt.startbyte, ""))
        elseif kindof(ps.nt) === Tokens.SEMICOLON
            next(ps)
            push!(top, EXPR(:toplevel, EXPR[EXPR(:NOTHING, ps.nt.startbyte, ps.nt.startbyte, "")]))
        end

        prevpos = position(ps)
        while kindof(ps.nt) !== Tokens.ENDMARKER
            curr_line = ps.nt.startpos[1]
            ret = parse_doc(ps)
            # join semicolon sep items
            if curr_line == last_line && headof(last(top.args)) === :toplevel
                push!(last(top.args), ret)
                top.fullspan += ret.fullspan
                top.span = top.fullspan - (ret.fullspan - ret.span)
            elseif kindof(ps.ws) == SemiColonWS
                push!(top, EXPR(:toplevel, EXPR[ret]))
            else
                push!(top, ret)
            end
            last_line = curr_line
            kindof(ps.nt) === Tokens.ENDMARKER && break # don't do loop check if eof
            prevpos = loop_check(ps, prevpos)
        end
    else
        if kindof(ps.nt) === Tokens.WHITESPACE || kindof(ps.nt) === Tokens.COMMENT
            next(ps)
            top = EXPR(:NOTHING, ps.nt.startbyte, ps.nt.startbyte, "")
        elseif !(ps.done || kindof(ps.nt) === Tokens.ENDMARKER)
            last_line = current_line(ps)
            if ps.nt.kind === Tokens.SEMICOLON
                next(ps)
                top = EXPR(:toplevel, EXPR[EXPR(:NOTHING, ps.nt.startbyte, ps.nt.startbyte, "")])
            else
                top = parse_doc(ps)
            end
            if kindof(ps.ws) == SemiColonWS# && curr_line == last_line
                top = EXPR(:toplevel, EXPR[top], nothing)
                prevpos = position(ps)
                while kindof(ps.ws) == SemiColonWS && current_line(ps) == last_line && kindof(ps.nt) != Tokens.ENDMARKER
                    last_line = current_line(ps)
                    ret = parse_doc(ps)
                    push!(top, ret)
                    prevpos = loop_check(ps, prevpos)
                end
            end
        else
            top = EXPR(:errortoken, EXPR[], nothing, 0, 0)
        end
    end

    return top, ps
end

function _continue_doc_parse(ps::ParseState, x::EXPR)
    headof(x) === :macrocall &&
    valof(x.args[1]) == "@doc" &&
    length(x.args) < 4 &&
    ps.t.endpos[1] + 1 == ps.nt.startpos[1]
end

include("precompile.jl")
_precompile()
end


text:
using HDF5
using Test

@testset "filter" begin

H5Z_FILTER_DEFLATE = 1
H5Z_FILTER_SHUFFLE = 2

# Create a new file
fn = tempname()

# Create test data
data = rand(1000, 1000)

# Open temp file for writing
f = h5open(fn, "w")

# Create datasets
dsdeflate = create_dataset(f, "deflate", datatype(data), dataspace(data),
                           chunk=(100, 100), compress=3)

dsshufdef = create_dataset(f, "shufdef", datatype(data), dataspace(data),
                           chunk=(100, 100), shuffle=(), compress=3)

dsfiltdef = create_dataset(f, "filtdef", datatype(data), dataspace(data),
                           chunk=(100, 100), filter=(H5Z_FILTER_DEFLATE,0,3))

dsfiltshufdef = create_dataset(f, "filtshufdef", datatype(data), dataspace(data),
                               chunk=(100, 100), filter=((H5Z_FILTER_SHUFFLE,0),
                                                         (H5Z_FILTER_DEFLATE,0,3)))

# Write data
write(dsdeflate, data)
write(dsshufdef, data)
write(dsfiltdef, data)
write(dsfiltshufdef, data)

# Close and re-open file for reading
close(f)
f = h5open(fn)

# Read dataseta
datadeflate = f["deflate"][]
datashufdef = f["shufdef"][]
datafiltdef = f["filtdef"][]
datafiltshufdef = f["filtshufdef"][]

close(f)

# Test for equality
@test datadeflate == data
@test datashufdef == data
@test datafiltdef == data
@test datafiltshufdef == data

end # @testset "filter"


text:
# This runs the forward model to make an ice thickness and then inverts for it.

# TODO:
# - run several tests
# - make into a test suite
# - figure out whether dropping temp or fsl is better

error("not working")

println("Loading BITEModel.jl")
using Optim
@time using BITEModel
const BM=BITEModel
import NLopt

plotyes = false
verbose = false
gid = SyntheticGlacier(:bench)

# Run it
########
println("Time to to make load-paras")

# parameters to play with
@show meas_dens = [:low, :normal, :high][2]
@show fsl0 = [0.1, 0.5, 0.9][1]
@show temp0 = [-10, -0.5][2]
@show sigma_small = 0.5

# TODO: This is bad as gl,pp,pm,pn,pl get loaded again below!
# I am not sure how the parameters flow...
# However, I also need gl in the
pl = BM.LoadPara(gid)
dt, pl = BM.make_datatable(gid,pl)
# Change defaults for :temp and/or :fsl
dt.fsl.opts[:fsl] = fsl0
dt.temp.opts[:temp] = temp0
gl,pp,pm,pn,pl = load_glacier(dt,pl)

# add the parameter to the options so the
# forward model is used to produce the thickness:
error_fac = 0

h_opts = pl.dataset_opts[BM.ThicknessData]
h_opts[:use_forward_model] = true
if meas_dens==:low
    h_opts[:profile_xlocs] = [2000.0]
    h_opts[:profile_ylocs] = [0.0]
elseif meas_dens==:high
    h_opts[:profile_xlocs] = [100:100:5900;]
    h_opts[:profile_ylocs] = [-250:250:250;]
end
h_opts[:fwd_para] = (gl, pp, pm, pn)
h_opts[:err_sigma] = 10.0 * error_fac
h_opts[:err_cor_length] = 600.0
h_opts[:profile_margin_dist] = 10.0

iv_opts = pl.dataset_opts[BM.IVData]
iv_opts[:use_forward_model] = true
iv_opts[:fwd_para] = (gl, pp, pm, pn)
iv_opts[:err_sigma] = 1.0 * error_fac
iv_opts[:err_cor_length] = 500.0
if meas_dens==:high
    iv_opts[:iv_grid] = (gl.dem.x, gl.dem.y)
elseif meas_dens==:low
    iv_opts[:iv_grid] = (gl.dem.x[1:10:end], gl.dem.y[1:10:end])
end
# for convenience:
h_opts_dt = dt.thickness.opts
iv_opts_dt = dt.iv.opts

# Now update all
dt, pl = BM.make_datatable(gid,pl)
# Again change defaults for :temp and/or :fsl
dt.fsl.opts[:fsl] = fsl0
dt.temp.opts[:temp] = temp0
gl,pp,pm,pn,pl = load_glacier(dt,pl)
gb,pm = BM.make_bands(gl, pp, pm, pn);

# Now, finally update pm (only benign things!!!)
## set model uncertainty to almost 0
pm = BM.MPara(pm, sigma_h_model=sigma_small, sigma_iv_model=sigma_small)

print("Time to run forward model:"); gc()
@time fwdsol = BM.fwdm(gb, pp, pm, pn, gb.bdot-gb.dhdt, gb.fsl, gb.temp)

qtot1d, qd1d = BM.postproc1d(gb, pp, pm)

if plotyes
    print("Time to plot:")
    @time display(BM.plot_run(fwdsol,
                              reuse=false, plot2d=[:h,:iv,:ivlog][1]))
    # @time display(BM.plot_run(gb, hs1d, taus1d, ivs1d, hs2d, taus2d,
    #                           ivs2d, taus1d_l, pp, pm,
    #                           reuse=false, plot2d=[:h,:iv,:ivlog][2]))
end

# Inversion
###########
n_1d_theta = 3
run=[:h, :iv, :h_iv][3]
(theta0, logposterior, logposterior1d, logprior, logprior_debug,
 loglikelihood, loglikelihood1d, fwdm_fn, fwdm1d_fn,
 pmcmc_defaults, fit_target) =
     BM.init_inverse(gb, pp, pm, pn, n_1d_theta=n_1d_theta, run=run,
                     theta0_dict = Dict(:btilde=> zeros(BM.F,n_1d_theta), # NOTE in m/y (error)
                                        :fsl=> zeros(BM.F,n_1d_theta), # error
                                        :temp=> zeros(BM.F,n_1d_theta), # error
                                        :dist_exp=> pm.dist_exp,
                                        :iv_h_exp=> pm.iv_h_exp,
                                        :iv_dist_exp1=> pm.iv_dist_exp1,
                                        #:iv_dist_exp2=> pm.iv_dist_exp2,
                                        )
                     )

th0_true = deepcopy(theta0.th0)

if h_opts[:err_sigma]==0 && iv_opts[:err_sigma]==0
    # check that th0_true is indeed minimal
    good = true
    maxll = loglikelihood(th0_true)[1]
    for delta = [1e-8, 1e-4, 1e-3,1e-2,1,10,100,1000]
        for i=1:length(th0_true)
            th = copy(th0_true)
            th[i] += delta
            if logprior(th)>-Inf && loglikelihood(th)[1]>maxll
                good = false
                warn("loglikelihood < for $th, delta=$delta, i=$i")
                println(loglikelihood(th)[1]-maxll)
            end
            th[i] -= 2delta
            if logprior(th)>-Inf && loglikelihood(th)[1]>maxll
                good = false
                warn("loglikelihood < for $th, delta=-$delta, i=$i")
                println(loglikelihood(th)[1]-maxll)
            end
        end
    end
    if good
        println("Log-likelihood indeed minimal at th0_true = $maxll")
    else
        println("Log-likelihood NOT minimal at th0_true = $maxll")
    end
else
    println("Log-likelihood at th0_true = $(loglikelihood(th0_true)[1])")
end

# some tests:
fwdm_fn(theta0.th0);
print("Time to run forward model:")
@time hs1d, taus1d, ivs1d, hs2d, taus2d, ivs2d, taus1d_l = fwdm_fn(theta0.th0);
print("Time to calculate posterior:")
@time logposterior(theta0.th0)[1]

print("Prior value:")
@show logprior(th0_true)
print("Likelihood value:")
@show loglikelihood(th0_true)[1]
print("Posterior value:")
@show logposterior(th0_true)[1]

########
# max posterior estimation
#######
println("Optimizing for maximum posterior...")

# # Using Optim.jl
# #
# # Work only reliably when close to the maximum.  Also the tolerances need to be reasonably low.
# #
# # - all converge if started at the minimum.  However Nealder-Mead takes 999 objective calls!
# # - Nealder-Mead claims to not be converged even though it gets closest to the max
# #
# opts = Optim.Options(x_tol=1e-8,
#                      f_tol=1e-8,
#                      g_tol=1e-7,
#                      iterations=1000)

# srand(1)
# ic = th0_true + rand(length(th0_true)) * 0.
# #ic = theta0.th0
# od = theta_vec -> logprior(theta_vec)[1]>-Inf ? -loglikelihood(theta_vec)[1] : Inf
# @show @time theta_max_likelihood1 = optimize(od, ic, opts)
# @show @time theta_max_likelihood2 = optimize(od, ic, GradientDescent(), opts)
# @show @time theta_max_likelihood2 = optimize(od, ic, AcceleratedGradientDescent(), opts)
# @show @time theta_max_likelihood3 = optimize(od, ic,  BFGS(), opts)

# println("Optimizing for maximum posterior...")
# @show @time theta_max_posterior = optimize(theta_vec -> -logposterior(theta_vec)[1], ic, BFGS(), opts)
# # these don't work so well:
# #   theta_max2 = optimize(theta_vec -> -logposterior(theta_vec)[1], float(theta0.th0), GradientDescent())
# #theta_max2 = optimize(theta_vec -> -logposterior(theta_vec)[1], float(theta0.th0),  LBFGS())

# using NLopt
#
# Much faster than Optim.jl. LN_SBPLX seems better than LN_NELDERMEAD

ic = th0_true + rand(length(th0_true)) * 0.01

od_ll = (th, grad) -> (logprior(th)[1]>-Inf ? loglikelihood(th)[1] : -Inf)
od_post = (th, grad) -> logposterior(th)[1]

# https://nlopt.readthedocs.io/en/latest/NLopt_Algorithms/#local-derivative-free-optimization
alg = [:LN_COBYLA, :LN_NELDERMEAD, :LN_BOBYQA, :LN_PRAXIS, :LN_SBPLX][end]
opt = NLopt.Opt(alg, length(ic))

NLopt.max_objective!(opt, od_post)
tol = 1e-6
NLopt.ftol_rel!(opt, tol)
NLopt.ftol_abs!(opt, tol)
NLopt.xtol_abs!(opt, tol)
NLopt.xtol_rel!(opt, tol)
NLopt.maxeval!(opt, 5000)

@time maxf, theta_start, ret = NLopt.optimize(opt, ic)
good_returns = [:SUCCESS, :STOPVAL_REACHED, :FTOL_REACHED, :XTOL_REACHED]
if !(ret in good_returns)
    println("Optimization not successful: $ret\n   Best max found=$maxf vs true $(logposterior(th0_true)[1])")
else
    println("Optimization successful.\n     Max found=$maxf vs true $(logposterior(th0_true)[1])")
end
theta0.th0[:] = theta_start

#############
# MCMC
############
pmcmc = BM.MCMCNum(niter = 10^3,
                   nthin = 20,
                   sigma_ppdf = ones(length(theta0.th0))*3e-2,
                   alg=[:emcee, :metro][1]) # metro needs many more iterations:
                                            # - emcee good with 10^5
                                            # - metro needs 10^?

# sample posterior
thetas, theta_expect, theta_mode, accept_ratio, blobs, lpost, mc_hs1d, mc_ivs1d,
  mh2d, sh2d, miv2d, siv2d, min_h2d, max_h2d, min_iv2d, max_iv2d = mcmc(logposterior,
                                                                        theta0, pmcmc; verbose=true)

# sample the prior to get that distribution too:
# pmcmc_prior = BM.MCMCNum(niter=10^4,
#                          nthin = 10)
# thetas_prior, theta_expect_prior, accept_ratio_prior = mcmc(logprior, theta0, pmcmc_prior; verbose=true)

# pmcmc_prior2 = BM.MCMCNum(niter=10^6,
#                           nthin = 1, alg=:metro,
#                           sigma_ppdf=ones(length(theta0.th0))/5  )
# thetas_prior2, theta_expect_prior, accept_ratio_prior = mcmc(logprior, theta0, pmcmc_prior2; verbose=true)

varnames = BM.get_varnames(theta0)
using KissMCMC
print_results(thetas, accept_ratio, names=varnames, theta_true=th0_true)

if plotyes
    display(BM.plotinv1d(gb, blobs, reuse=false))
    display(BM.plotinv1d_iverr(gb, blobs, reuse=false))

    display(BM.plotinv2d(gl, blobs, reuse=false))
    display(BM.plotinv2d_h(gl, blobs, reuse=false))
    display(BM.plotinv2d_iverr(gl, blobs, reuse=false))

    display(BM.plottheta(thetas, theta0, toplot=:all))
    # display(BM.plottheta(thetas, theta0, toplot=:btilde, reuse=false))
    # display(BM.plottheta(thetas, theta0, toplot=:fsl, reuse=false))
    # display(BM.plottheta(thetas, theta0, toplot=:temp, reuse=false))

    display(BM.plottheta_violin((thetas, thetas_prior), theta0, :btilde, reuse=false, width=1))
    display(BM.plottheta_violin((thetas, thetas_prior), theta0, :fsl, reuse=false, width=1))
    display(BM.plottheta_violin((thetas, thetas_prior), theta0, :temp, reuse=false, width=1))
end

nothing


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".

"""
    mean_obliquity(jd) -> m_eps

### Purpose ###

Return the mean obliquity of the ecliptic for a given Julian date

### Explanation ###

The function is used by the [`co_nutate`](@ref) procedure.

### Arguments ###

* `jd`: julian date

### Output ###

* `m_eps`: mean obliquity of the ecliptic, in radians

### Example ###

```jldoctest
julia> using AstroLib

julia> mean_obliquity(jdcnv(1978,01,7,11, 01))
0.4091425159336512
```

### Notes ###

The algorithm used to find the mean obliquity(`eps0`) is mentioned in
USNO Circular 179, but the canonical reference for the IAU adoption
is apparently Hilton et al., 2006, Celest.Mech.Dyn.Astron. 94, 351. 2000
"""
function mean_obliquity(jd::Real)
    t = (jd - J2000) / JULIANCENTURY
    eps0 = @evalpoly t 84381.406 -46.836769 -0.0001831 0.0020034 -0.576e-6 -4.34e-8
    return sec2rad(eps0)
end


text:
requiredFields = ["byr", "iyr", "eyr", "hgt", "hcl", "ecl", "pid"]

mutable struct Passport
    byr:: Int
    iyr:: Int
    eyr:: Int
    hgt:: Int
    hcl:: Int
    ecl:: Int
    pid:: Int
    isValid:: Bool
    Passport() = new()
end

function addLineToPassport(line, passport)
    fieldes = split(line," ")
    for field in fieldes
        (key,value) = split(field, ":")
        if key == "byr" passport.byr = 1 end
        if key == "iyr" passport.iyr = 1 end
        if key == "eyr" passport.eyr = 1 end
        if key == "hgt" passport.hgt = 1 end
        if key == "hcl" passport.hcl = 1 end
        if key == "ecl" passport.ecl = 1 end
        if key == "pid" passport.pid = 1 end
    end
    return passport
end

function validatePassport(p::Passport)
    if(p.byr + p.iyr + p.eyr + p.hgt + p.hcl + p.ecl + p.pid) == 7
        p.isValid = true
    end
end

function collectPassports(fileName)
    passports = []
    passportCounter = 1
    currentPassport = Passport()
    open(fileName) do file
        for line in eachline(file)
            if length(line) == 0
                validatePassport(currentPassport)
                push!(passports, currentPassport)
                currentPassport = Passport()
                passportCounter += 1
            else
                addLineToPassport(line, currentPassport)
            end
        end
    end
    return passports
end

allPassports = collectPassports("data.txt")
numberOfValidPassports = 0
for passport in allPassports
    if passport.isValid
        global numberOfValidPassports += 1
    end
end
println(numberOfValidPassports)

# passport = Passport()
# passport.byr = 1
# passport.iyr = 1
# passport.eyr = 1
# passport.hgt = 1
# passport.hcl = 1
# validatePassport(passport)
# println(passport.isValid)

# line = "hcl:1f7352 iyr:2014 hgt:70cm eyr:1983"
# passport = Passport()
# println(passport)
# pasport = addLineToPassport(line, passport)
# println(passport)

# passports = collectPassports("data_test.txt")
# for passport in passports
#     println(passport)
# end