text:
using DataDeps

register(DataDep("OpenOffice-MyThes-1.0",
    """
    Dataset: OpenOffice MyThes-1.0
    Website: https://www.openoffice.org/lingucomponent/thesaurus.html
    Author: Kevin Hendricks <kevin.hendricks@sympatico.ca>

    Standalone version of the MyThes thesaurus code. This includes a thesaurus for en_US in its new format for OOo 2.0.

    License:
    /*
     * Copyright 2003 Kevin B. Hendricks, Stratford, Ontario, Canada
     * And Contributors.  All rights reserved.
     *
     * Redistribution and use in source and binary forms, with or without
     * modification, are permitted provided that the following conditions
     * are met:
     *
     * 1. Redistributions of source code must retain the above copyright
     *    notice, this list of conditions and the following disclaimer.
     *
     * 2. Redistributions in binary form must reproduce the above copyright
     *    notice, this list of conditions and the following disclaimer in the
     *    documentation and/or other materials provided with the distribution.
     *
     * 3. All modifications to the source code must be clearly marked as
     *    such.  Binary redistributions based on modified source code
     *    must be clearly marked as modified versions in the documentation
     *    and/or other materials provided with the distribution.
     *
     * THIS SOFTWARE IS PROVIDED BY KEVIN B. HENDRICKS AND CONTRIBUTORS
     * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
     * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
     * FOR A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL
     * KEVIN B. HENDRICKS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
     * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
     * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
     * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
     * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
     * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
     * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
     * SUCH DAMAGE.
     *
     */
    """,
    "https://github.com/rdeits/openoffice-mythes-mirror/archive/v1.0.zip",
    "5dddb4ee551bb9f850932650da039189501e70a20dfa1663519e46bda78706ab",
    post_fetch_method = unpack
))

register(DataDep("SCOWL-wordlist-en_US-large",
    """
    Dataset: SCOWL word list (en_US dialect, large size)
    Website: http://wordlist.aspell.net/
    Author: Kevin Atkinson <kevina@gnu.org>

    COPYRIGHT, SOURCES, and CREDITS:
    The English dictionaries come directly from SCOWL
    and is thus under the same copyright of SCOWL.  The affix file is
    a heavily modified version of the original english.aff file which was
    released as part of Geoff Kuenning's Ispell and as such is covered by
    his BSD license.  Part of SCOWL is also based on Ispell thus the
    Ispell copyright is included with the SCOWL copyright.

    The collective work is Copyright 2000-2016 by Kevin Atkinson as well
    as any of the copyrights mentioned below:

      Copyright 2000-2016 by Kevin Atkinson

      Permission to use, copy, modify, distribute and sell these word
      lists, the associated scripts, the output created from the scripts,
      and its documentation for any purpose is hereby granted without fee,
      provided that the above copyright notice appears in all copies and
      that both that copyright notice and this permission notice appear in
      supporting documentation. Kevin Atkinson makes no representations
      about the suitability of this array for any purpose. It is provided
      "as is" without express or implied warranty.

    Alan Beale <biljir@pobox.com> also deserves special credit as he has,
    in addition to providing the 12Dicts package and being a major
    contributor to the ENABLE word list, given me an incredible amount of
    feedback and created a number of special lists (those found in the
    Supplement) in order to help improve the overall quality of SCOWL.

    The 10 level includes the 1000 most common English words (according to
    the Moby (TM) Words II [MWords] package), a subset of the 1000 most
    common words on the Internet (again, according to Moby Words II), and
    frequently class 16 from Brian Kelk's "UK English Wordlist
    with Frequency Classification".

    The MWords package was explicitly placed in the public domain:

        The Moby lexicon project is complete and has
        been place into the public domain. Use, sell,
        rework, excerpt and use in any way on any platform.

        Placing this material on internal or public servers is
        also encouraged. The compiler is not aware of any
        export restrictions so freely distribute world-wide.

        You can verify the public domain status by contacting

        Grady Ward
        3449 Martha Ct.
        Arcata, CA  95521-4884

        grady@netcom.com
        grady@northcoast.com

    The "UK English Wordlist With Frequency Classification" is also in the
    Public Domain:

      Date: Sat, 08 Jul 2000 20:27:21 +0100
      From: Brian Kelk <Brian.Kelk@cl.cam.ac.uk>

      > I was wondering what the copyright status of your "UK English
      > Wordlist With Frequency Classification" word list as it seems to
      > be lacking any copyright notice.

      There were many many sources in total, but any text marked
      "copyright" was avoided. Locally-written documentation was one
      source. An earlier version of the list resided in a filespace called
      PUBLIC on the University mainframe, because it was considered public
      domain.

      Date: Tue, 11 Jul 2000 19:31:34 +0100

      > So are you saying your word list is also in the public domain?

      That is the intention.

    The 20 level includes frequency classes 7-15 from Brian's word list.

    The 35 level includes frequency classes 2-6 and words appearing in at
    least 11 of 12 dictionaries as indicated in the 12Dicts package.  All
    words from the 12Dicts package have had likely inflections added via
    my inflection database.

    The 12Dicts package and Supplement is in the Public Domain.

    The WordNet database, which was used in the creation of the
    Inflections database, is under the following copyright:

      This software and database is being provided to you, the LICENSEE,
      by Princeton University under the following license.  By obtaining,
      using and/or copying this software and database, you agree that you
      have read, understood, and will comply with these terms and
      conditions.:

      Permission to use, copy, modify and distribute this software and
      database and its documentation for any purpose and without fee or
      royalty is hereby granted, provided that you agree to comply with
      the following copyright notice and statements, including the
      disclaimer, and that the same appear on ALL copies of the software,
      database and documentation, including modifications that you make
      for internal use or for distribution.

      WordNet 1.6 Copyright 1997 by Princeton University.  All rights
      reserved.

      THIS SOFTWARE AND DATABASE IS PROVIDED "AS IS" AND PRINCETON
      UNIVERSITY MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR
      IMPLIED.  BY WAY OF EXAMPLE, BUT NOT LIMITATION, PRINCETON
      UNIVERSITY MAKES NO REPRESENTATIONS OR WARRANTIES OF MERCHANT-
      ABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE OF THE
      LICENSED SOFTWARE, DATABASE OR DOCUMENTATION WILL NOT INFRINGE ANY
      THIRD PARTY PATENTS, COPYRIGHTS, TRADEMARKS OR OTHER RIGHTS.

      The name of Princeton University or Princeton may not be used in
      advertising or publicity pertaining to distribution of the software
      and/or database.  Title to copyright in this software, database and
      any associated documentation shall at all times remain with
      Princeton University and LICENSEE agrees to preserve same.

    The 40 level includes words from Alan's 3esl list found in version 4.0
    of his 12dicts package.  Like his other stuff the 3esl list is also in the
    public domain.

    The 50 level includes Brian's frequency class 1, words appearing
    in at least 5 of 12 of the dictionaries as indicated in the 12Dicts
    package, and uppercase words in at least 4 of the previous 12
    dictionaries.  A decent number of proper names is also included: The
    top 1000 male, female, and Last names from the 1990 Census report; a
    list of names sent to me by Alan Beale; and a few names that I added
    myself.  Finally a small list of abbreviations not commonly found in
    other word lists is included.

    The name files form the Census report is a government document which I
    don't think can be copyrighted.

    The file special-jargon.50 uses common.lst and word.lst from the
    "Unofficial Jargon File Word Lists" which is derived from "The Jargon
    File".  All of which is in the Public Domain.  This file also contain
    a few extra UNIX terms which are found in the file "unix-terms" in the
    special/ directory.

    The 55 level includes words from Alan's 2of4brif list found in version
    4.0 of his 12dicts package.  Like his other stuff the 2of4brif is also
    in the public domain.

    The 60 level includes all words appearing in at least 2 of the 12
    dictionaries as indicated by the 12Dicts package.

    The 70 level includes Brian's frequency class 0 and the 74,550 common
    dictionary words from the MWords package.  The common dictionary words,
    like those from the 12Dicts package, have had all likely inflections
    added.  The 70 level also included the 5desk list from version 4.0 of
    the 12Dics package which is in the public domain.

    The 80 level includes the ENABLE word list, all the lists in the
    ENABLE supplement package (except for ABLE), the "UK Advanced Cryptics
    Dictionary" (UKACD), the list of signature words from the YAWL package,
    and the 10,196 places list from the MWords package.

    The ENABLE package, mainted by M Cooper <thegrendel@theriver.com>,
    is in the Public Domain:

      The ENABLE master word list, WORD.LST, is herewith formally released
      into the Public Domain. Anyone is free to use it or distribute it in
      any manner they see fit. No fee or registration is required for its
      use nor are "contributions" solicited (if you feel you absolutely
      must contribute something for your own peace of mind, the authors of
      the ENABLE list ask that you make a donation on their behalf to your
      favorite charity). This word list is our gift to the Scrabble
      community, as an alternate to "official" word lists. Game designers
      may feel free to incorporate the WORD.LST into their games. Please
      mention the source and credit us as originators of the list. Note
      that if you, as a game designer, use the WORD.LST in your product,
      you may still copyright and protect your product, but you may *not*
      legally copyright or in any way restrict redistribution of the
      WORD.LST portion of your product. This *may* under law restrict your
      rights to restrict your users' rights, but that is only fair.

    UKACD, by J Ross Beresford <ross@bryson.demon.co.uk>, is under the
    following copyright:

      Copyright (c) J Ross Beresford 1993-1999. All Rights Reserved.

      The following restriction is placed on the use of this publication:
      if The UK Advanced Cryptics Dictionary is used in a software package
      or redistributed in any form, the copyright notice must be
      prominently displayed and the text of this document must be included
      verbatim.

      There are no other restrictions: I would like to see the list
      distributed as widely as possible.

    The 95 level includes the 354,984 single words, 256,772 compound
    words, 4,946 female names and the 3,897 male names, and 21,986 names
    from the MWords package, ABLE.LST from the ENABLE Supplement, and some
    additional words found in my part-of-speech database that were not
    found anywhere else.

    Accent information was taken from UKACD.

    The VarCon package was used to create the American, British, Canadian,
    and Australian word list.  It is under the following copyright:

      Copyright 2000-2016 by Kevin Atkinson

      Permission to use, copy, modify, distribute and sell this array, the
      associated software, and its documentation for any purpose is hereby
      granted without fee, provided that the above copyright notice appears
      in all copies and that both that copyright notice and this permission
      notice appear in supporting documentation. Kevin Atkinson makes no
      representations about the suitability of this array for any
      purpose. It is provided "as is" without express or implied warranty.

      Copyright 2016 by Benjamin Titze

      Permission to use, copy, modify, distribute and sell this array, the
      associated software, and its documentation for any purpose is hereby
      granted without fee, provided that the above copyright notice appears
      in all copies and that both that copyright notice and this permission
      notice appear in supporting documentation. Benjamin Titze makes no
      representations about the suitability of this array for any
      purpose. It is provided "as is" without express or implied warranty.

      Since the original words lists come from the Ispell distribution:

      Copyright 1993, Geoff Kuenning, Granada Hills, CA
      All rights reserved.

      Redistribution and use in source and binary forms, with or without
      modification, are permitted provided that the following conditions
      are met:

      1. Redistributions of source code must retain the above copyright
         notice, this list of conditions and the following disclaimer.
      2. Redistributions in binary form must reproduce the above copyright
         notice, this list of conditions and the following disclaimer in the
         documentation and/or other materials provided with the distribution.
      3. All modifications to the source code must be clearly marked as
         such.  Binary redistributions based on modified source code
         must be clearly marked as modified versions in the documentation
         and/or other materials provided with the distribution.
      (clause 4 removed with permission from Geoff Kuenning)
      5. The name of Geoff Kuenning may not be used to endorse or promote
         products derived from this software without specific prior
         written permission.

      THIS SOFTWARE IS PROVIDED BY GEOFF KUENNING AND CONTRIBUTORS ``AS IS'' AND
      ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
      IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
      ARE DISCLAIMED.  IN NO EVENT SHALL GEOFF KUENNING OR CONTRIBUTORS BE LIABLE
      FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
      DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
      OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
      HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
      LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
      OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
      SUCH DAMAGE.

    Build Date: Mon Apr 16 22:14:07 EDT 2018
    Wordlist Command: mk-list -v1 --accents=both en_US 70
    """,
    "https://github.com/rdeits/speller-wordlist-en-US-large-mirror/raw/8d85394ac1b38383c6fc352d5cd5f4994ed068ca/wordlist-en_US-large-2018.04.16.zip",
    "136da8bf589f26a04c45c36ef204b983c4941ad74a44006de74dd235ca3a311b",
    post_fetch_method = unpack
))


text:
using Test
using ASE
using LinearAlgebra
import Random
using JuLIP
using DataDeps
using Optim

data = joinpath(dirname(pathof(JuLIP)), "..", "data") * "/"

##
# @testset "EAM with ASE" begin

test_pots = joinpath(datadep"JuLIP_testpots", "JuLIP_data") * "/"
alloy = test_pots * "PdAgH_HybridPd3Ag.eam.alloy"
eam_fs = test_pots * "Fe-P.eam.fs"
Ni = test_pots * "Ni.eam.fs"

@test EAM(alloy) isa EAM{T, 1} where {T}
@test EAM(eam_fs) isa EAM{T, 2} where {T}

eam = pyimport("ase.calculators.eam")
ase_calc = ASECalculator(eam.EAM(potential=Ni))
# make sure we get a perfect fit ...
julip_calc = EAM(Ni)

atoms = bulk(:Ni) * 3
Random.seed!(0)
rattle!(atoms, 0.1)
E_jl = energy(julip_calc, atoms)
E_py = energy(ase_calc, atoms)
# ... but even then we the evaluation codes aren't the
#     same so we only get ca 1e-6 to 1e-7 match.
@test abs(E_jl - E_py) < 1e-10

# Test again with unordered species in parameter file
ase_calc = ASECalculator(eam.EAM(potential=alloy))
julip_calc = EAM(alloy)

atoms = bulk(:Pd) * 3
atoms.Z[1:3:end] .= AtomicNumber(:Ag)
atoms.Z[1:5:end] .= AtomicNumber(:H)
rattle!(atoms, 0.1)
E_jl = energy(julip_calc, atoms)
E_py = energy(ase_calc, atoms)
@test E_jl ≈ E_py rtol=1e-6
# end

# Test a few quantities of interest with respect to reference values obtained with LAMMPS
# You can set `got_lammps`to true if you have ASE and LAMMPS Python libraries installed

got_lammps = false
pot_file = test_pots * "w_eam4.fs"
eam_julip = EAM(pot_file)

if got_lammps
    lammpslib = pyimport("ase.calculators.lammpslib")
    eam_lammps = ASECalculator(lammpslib.LAMMPSlib(lmpcmds=["pair_style eam/fs",
                                "pair_coeff * * $pot_file W"],
                                atom_types=Dict("W" => 1), keep_alive=true))
end

# simple energy test
unitcell = bulk(:W, cubic=true)
@test energy(eam_julip, unitcell) ≈ -17.792788114962622
if got_lammps
    @show energy(eam_julip, unitcell)
    @show energy(eam_lammps, unitcell)
    @test energy(eam_julip, unitcell) ≈ energy(eam_lammps, unitcell)
end

# lattice constant
eam = eam_julip # or eam_lammps
variablecell!(unitcell)
set_calculator!(unitcell, eam)
res = optimize(x -> energy(unitcell, x), x -> gradient(unitcell, x),
                dofs(unitcell), inplace=false)
alat = unitcell.cell[1, 1]
@test alat ≈ 3.1433900286308045

# (110) surface energy
shift = 2.0
ase_lattice_cubic = pyimport("ase.lattice.cubic")
make_bulk() = Atoms(ASEAtoms(ase_lattice_cubic.BodyCenteredCubic(symbol="W",
                             latticeconstant=alat,
                             directions=[[1,-1,0],[1,1,0],[0,0,1]]) * (1, 1, 10)))
bulk_at = make_bulk()
set_calculator!(bulk_at, eam)
surface = make_bulk() #copy(bulk_at) # is there a better way to do this in JuLIP?
X = positions(surface) |> mat
X[3, :] .+= shift
set_positions!(surface, X)
wrap_pbc!(surface)
c = Matrix(surface.cell)
c[3, :] += [0.0, 0.0, 10.0]
set_cell!(surface, c)
fixedcell!(surface)
set_calculator!(surface, eam)
area = norm(cross(bulk_at.cell[:, 1], bulk_at.cell[:, 2]))
γ = (energy(surface) - energy(bulk_at)) / (2 * area)
@show γ
@test γ ≈ 0.18415972172748984

ext: jl


text:
"""
    QINT{N, AM} <: AbstractContinuousPost

Reachability method for one-dimensional quadratic ODEs with non-deterministic inputs.

## Fields

The type fields are:

- `N`   -- number type of the step-size
- `AM`  -- type of the approximation model

## Notes

"""
struct QINT{N, AM} <: AbstractContinuousPost
    Δ::N
    δ::N
    θ::N
    maxiter::Int
    approx_model::AM
end

# convenience constructor using symbols
function QINT(; Δ::N, δ::N, θ::N, maxiter::Int,
               approx_model::AM=Forward(sih=:concrete, exp=:base, setops=:interval)) where {N, AM}
    return QINT(Δ, δ, θ, maxiter, approx_model)
end

step_size(alg::QINT) = alg.δ
numtype(::QINT{N}) where {N} = N
setrep(::QINT{N}) where {N} = Interval{N, IA.Interval{N}}
rsetrep(::QINT{N}) where {N} = ReachSet{N, Interval{N, IA.Interval{N}}}

include("post.jl")
include("reach_homog.jl")
#include("reach_inhomog.jl")


text: println(sum(filter(x -> x%3 == 0 || x%5 == 0, 1:999)))

max_stars_repo_name: sujimodern/project-euler-jl

max_issues_repo_name: sujimodern/project-euler-jl

max_forks_repo_name: sujimodern/project-euler-jl


text:
using NASAPrecipitation
using Documenter

makedocs(;
    modules  = [NASAPrecipitation],
    doctest  = false,
    format   = Documenter.HTML(collapselevel=1,prettyurls=false),
    authors  = "Nathanael Wong <natgeo.wong@outlook.com>",
    sitename = "NASAPrecipitation.jl",
    pages    = [
        "Home"                        => "index.md",
        "NASAPrecipitation Datasets" => [
            "The NASAPrecipitation Dataset Type" => "datasets/npd.md",
            "GPM IMERG Datasets"             => "datasets/imerg.md",
            "TRMM TMPA Datasets"             => "datasets/trmm.md",
            "Land-Sea Mask Datasets"         => "datasets/landseamask.md",
        ],
        "A Basic Primer to `GeoRegion`s" => "georegions/basics.md",
        "Downloading Datasets" => "download.md",
        "Reading Downloaded Data" => "read.md",
        "Examples" => [
            "Downloading IMERG Data"  => "examples/download.md",
            "Land-Sea Mask Filtering" => "examples/landseamask.md",
        ],
    ],
)

deploydocs(
    repo = "github.com/natgeo-wong/NASAPrecipitation.jl.git",
    devbranch = "main"
)

max_stars_repo_name: natgeo-wong/NASAPrecipitation

max_issues_repo_name: natgeo-wong/NASAPrecipitation

max_forks_repo_name: natgeo-wong/NASAPrecipitation


text: int(x) = (Int ∘ round)(x)


text:
origin(prim::HyperCube) = prim.origin
width(prim::HyperCube)  = prim.width
widths(prim::HyperCube{N,T}) where {N,T} = Vec{N,T}(prim.width)
maximum(prim::HyperCube{N,T}) where {N,T} = origin(prim)+widths(prim)
minimum(prim::HyperCube{N,T}) where {N,T} = origin(prim)
centered(C::Type{HyperCube{N,T}}) where {N,T} = C(Vec{N,T}(-0.5), T(1))
function centered(::Type{T}) where T <: HyperCube
    centered(HyperCube{ndims_or(T, 3), eltype_or(T, Float32)})
end


text:
using ModelingToolkit, Catalyst #, LinearAlgebra
#using DiffEqBase, DiffEqJump
#using Plots, SpecialFunctions

## Parameter
N = 10                       # maximum cluster size
Vₒ = (4π/3)*(10e-06*100)^3   # volume of a monomers in cm³
Nₒ = 1e-06/Vₒ                # initial conc. = (No. of init. monomers) / bulk volume
uₒ = 10000                   # No. of monomers initially
V = uₒ/Nₒ                    # Bulk volume of system in cm³

integ(x) = Int(floor(x))
n        = integ(N/2)
nr       = N%2 == 0 ? (n*(n + 1) - n) : (n*(n + 1)) # No. of forward reactions

# possible pairs of reactant multimers
pair = []
for i = 2:N
    push!(pair,[1:integ(i/2)  i .- (1:integ(i/2))])
end
pair = vcat(pair...)
vᵢ = @view pair[:,1]   # Reactant 1 indices
vⱼ = @view pair[:,2]   # Reactant 2 indices
volᵢ = Vₒ*vᵢ           # cm⁻³
volⱼ = Vₒ*vⱼ           # cm⁻³
sum_vᵢvⱼ = @. vᵢ + vⱼ  # Product index

# set i to  1 for additive kernel, 2  for constant
i = 1
if i==1
    B = 1.53e03                # s⁻¹
    kv = @. B*(volᵢ + volⱼ)/V  # dividing by volume as its a bi-molecular reaction chain
elseif i==2
    C = 1.84e-04               # cm³ s⁻¹
    kv = fill(C/V, nr)
end

# state variables are X, pars stores rate parameters for each rx
@parameters t
@variables k[1:nr]  X[1:N](t)
pars = Pair.(collect(k), kv)

# time-span
if i == 1
    tspan = (0. ,2000.)
elseif i == 2
    tspan = (0. ,350.)
end

 # initial condition of monomers
u₀    = zeros(Int64, N)
u₀[1] = uₒ
u₀map = Pair.(collect(X), u₀)   # map variable to its initial value

# vector to store the Reactions in
rx = []
for n = 1:nr
    # for clusters of the same size, double the rate
    if (vᵢ[n] == vⱼ[n])
        push!(rx, Reaction(k[n], [X[vᵢ[n]]], [X[sum_vᵢvⱼ[n]]], [2], [1]))
    else
        push!(rx, Reaction(k[n], [X[vᵢ[n]], X[vⱼ[n]]], [X[sum_vᵢvⱼ[n]]],
                           [1, 1], [1]))
    end
end
@named rs = ReactionSystem(rx, t, collect(X), collect(k))
Catalyst.Graph(rs)


text:
update_loss_log(run_logs, loss) = append!(run_logs["loss"], loss)
Zygote.@nograd update_loss_log

function train(prob, t, u0, p3, ps, data, epochs, opt, avg_meter,
               plen; cont_train = false)
    log_infected = log.(abs.(data["infected"]) .+ 1.0)
    log_recovered = log.(abs.(data["recovered"] .+ data["deceased"]) .+ 1.0)
    weightage = data["time_train"] ./ length(data["time_train"])

    run_logs = Dict(
        "loss" => [],
        "p1"   => [],
        "p2"   => [],
        "p3"   => [],
    )

    dat = Iterators.repeated((), epochs)

    function predict_adjoint()
        Array(concrete_solve(prob, Rosenbrock23(autodiff = false), u0, p3, saveat=t))
    end

    if !cont_train
        reset!(avg_meter)
    end

    function loss_adjoint()
        prediction = predict_adjoint()
        infect_loss = 3.0 * sum(weightage .* abs2.(log_infected .- log.(abs.(prediction[2, :]) .+ 1.0)))
        recover_loss = sum(weightage .* abs2.(log_recovered .- log.(abs.(prediction[3, :]) .+ 1.0)))
        loss = infect_loss + recover_loss

        update(avg_meter, loss)
        update_loss_log(run_logs, loss)

        return loss
    end

    loss_adjoint()

    cb = function()
        println("Running Mean Loss : $(avg_meter.val)")

        append!(run_logs["p1"], p3[plen + 1])
        append!(run_logs["p2"], p3[plen + 2])
        append!(run_logs["p3"], p3)
    end

    Flux.train!(loss_adjoint, ps, dat, opt, cb = cb)

    # min_loss, idx = findmin(run_logs["loss"])
    # idx1 = (idx - 1) * (plen + 2) + 1
    # idx2 = idx * (plen + 2)
    # p3 = run_logs["p3"][idx1:idx2]

    return u0, p3
end

function inference(prob, u0, p, t)
    prediction = Array(concrete_solve(prob, Rosenbrock23(autodiff = false), u0, p, saveat=t))

    S_NN = prediction[1, :]
    I_NN = prediction[2, :]
    R_NN = prediction[3, :]
    T_NN = prediction[4, :]

    return S_NN, I_NN, R_NN, T_NN
end

function train_named_locations(file::String, districts, states; each_epochs = 200, retrain = false)
    data_all = clean_data(file)
    train_named_locations(data_all, districts, states; each_epochs = each_epochs, retrain = retrain)
end

function train_named_locations(data_all, districts, states; each_epochs = 200, retrain = false)
    avg_meter = RunningAverageMeter()
    districts = format_string.(districts)
    for state in format_string.(states)
        data_state = data_all[state]
        @info "Starting to train models for State: $state"
        for district in keys(data_state)
            if length(districts) > 0 && district ∉ districts
                continue
            end
            df = data_state[district]
            if isnothing(df)
                continue
            end
            filepath = joinpath(@__DIR__, "../data/" * state * "_" * district * ".bson")
            if isfile(filepath) && !retrain
                @info "Model already trained for District $district"
                continue
            end
            @info "District $district"
            t_train = df["time_train"]
            t_test = df["time_test"]

            u0_1 = df["infected"][1]
            u0_2 = df["recovered"][1] + df["deceased"][1]
            ps, re, u0, prob, p3, plen =
                construct_model(t_train, t_test, u0_1, u0_2, hidden_dim = 30,
                                test_model = false)
            opt = ADAM(1e-2)
            @show df["infected"], df["recovered"] .+ df["deceased"], df["time_train"]
            u0, p = train(prob, t_train, u0, p3, ps, df, each_epochs, opt, avg_meter, plen,
                          cont_train = false)

            @show inference(prob, u0, p, t_train)[2:3]
            @show inference(prob, u0, p, t_test)[2:3]
            BSON.@save filepath u0 p t_train t_test
            @info "Saved model to $filepath"
        end
    end
end

function train_all_districts(file; each_epochs = 200, retrain = false)
    data_all = clean_data(file)
    train_named_locations(data_all, [], collect(keys(data_all)); each_epochs = each_epochs, retrain = retrain)
end

prediction(state::String, district::String, day_num::Int64) =
    prediction(state, district, [day_num])

function prediction(state::String, district::String, day_num::Vector{Int64})
    filepath = joinpath(@__DIR__, "../data/" * format_string(state) * "_" * format_string(district) * ".bson")
    if !isfile(filepath)
        return nothing
    end
    BSON.@load filepath u0 p t_train t_test
#     @show t_train, t_test
    _, _, _, prob, _, _ =
        construct_model(t_train, t_test, u0[2], u0[3], hidden_dim = 30, test_model = false)
    s_nn, i_nn, r_nn, t_nn = inference(prob, u0, p, t_train)
#     @show i_nn, r_nn
    s_nn, i_nn, r_nn, t_nn = inference(prob, u0, p, t_test)
    return i_nn[day_num], r_nn[day_num]
end


text:
module MetaGen

using Gen
using GenRFS
using Distributions
using Combinatorics
using Random
#using Plots
using PyPlot
const plt = PyPlot
#import PyPlot; const plt = PyPlot

include("declaring_structs.jl")
include("custom_distributions.jl")
include("geometry_optics.jl")
include("inverse_optics.jl")
include("receptive_fields.jl")
include("video.jl")
include("involution.jl")
include("particle_filter_helper_functions.jl")
include("pf_inference.jl")
include("main.jl")
include("printing.jl")
include("visualizations.jl")

function __init__()
    @load_generated_functions
end

end # module

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
#=
Note: This file was taken from the TimeWarp.jl package

The TimeWarp.jl package is licensed under the MIT "Expat" License:

> Copyright (c) 2016: Alex Williams and Joseph Fowler
>
> Permission is hereby granted, free of charge, to any person obtaining a copy
> of this software and associated documentation files (the "Software"), to deal
> in the Software without restriction, including without limitation the rights
> to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
> copies of the Software, and to permit persons to whom the Software is
> furnished to do so, subject to the following conditions:
>
> The above copyright notice and this permission notice shall be included in all
> copies or substantial portions of the Software.
>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
> IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
> FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
> AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
> LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
> OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
> SOFTWARE.
=#

using Base.Cartesian

"""
    Sequence{N,T}(::AbstractArray)

A `Sequence` is a thin wrapper around a multi-dimensional array that
allows you to index it like a vector. For example, the following
code shows a 10-dimensional time series, sampled at 100 time points:

```
data = randn(10,100)
seq = Sequence(data)
seq[1] == data[:,1] # true
```

This generalizes to higher-order arrays. Consider a time series
where a 10x10 matrix of data is collected at 100 time points:

```
data = randn(10,10,100)
seq = Sequence(data)
seq[1] == data[:,:,1] # true
```
"""
immutable Sequence{N,T} <: AbstractArray{T,1}
    val::AbstractArray{T,N}
end

Base.size{N}(x::Sequence{N}) = (size(x.val,N),)
Base.eltype{N,T}(x::Sequence{N,T}) = T

@generated function Base.getindex{N}(x::Sequence{N}, i)
    :( x.val[@ntuple($N, (n-> n==$N ? i : Colon()))...] )
end

@generated function Base.setindex!{N}(x::Sequence{N}, val, i)
    :( x.val[@ntuple($N, (n-> n==$N ? i : Colon()))...] = val )
end


text:
"""
	Implements all needed parts for constructing GANomaly model
		(https://arxiv.org/pdf/1805.06725.pdf) for detecting anomalies.
	Code is inspired by original Pytorch implementation https://github.com/samet-akcay/ganomaly.

"""

struct Encoder
	layers
end

"""
	function ConvEncoder(isize::Int, in_ch::Int, out_ch::Int, nf::Int, extra_layers::Int)

Create the convolutional encoder with parameters
	isize           - size of input image (must be divisible by 16 i.e isize%16==0)
	in_ch           - number of input channels
	out_ch          - number of output channels / dimension of latent space
	nf              - number of covnolutional masks/filters
	extra_layers    - number of additional conv blocks (conv2d + batch norm + leaky relu)
"""
function ConvEncoder(isize::Int, in_ch::Int, out_ch::Int, nf::Int, extra_layers::Int)
	# ověřit jestli isize%16==0
	layers = []
	push!(layers, Conv((4, 4), in_ch => nf; stride = 2, pad = 1))
	push!(layers, BatchNorm(nf))
	push!(layers, x->leakyrelu.(x, 0.2f0))

	csize, cnf = isize/2, nf
	for i=1:extra_layers
		push!(layers, Conv((3, 3), cnf => cnf; stride = 1, pad = 1))
		push!(layers, BatchNorm(cnf))
		push!(layers, x->leakyrelu.(x, 0.2f0))
	end

	while csize > 4
		in_feat = cnf
		out_feat = cnf*2
		push!(layers, Conv((4, 4), in_feat => out_feat; stride = 2, pad = 1))
		push!(layers, BatchNorm(out_feat))
		push!(layers, x->leakyrelu.(x, 0.2f0))
		cnf = cnf*2
		csize=csize/2
	end
	push!(layers, Conv((4, 4), cnf => out_ch; stride = 1, pad = 0))
	return Encoder(Chain(layers...))
end

Flux.@functor Encoder

""" Definition of Encoder's forward pass """
function (e::Encoder)(x)
	e.layers(x)
end

struct Decoder
	layers
end

"""
	function ConvDecoder(isize::Int, in_ch::Int, out_ch::Int, nf::Int, extra_layers::Int)

Create the convolutional decoder with parameters
	isize           - size of output image (must be divisible by 16 i.e isize%16==0)
	in_ch           - number of input channels / dimension of latent space
	out_ch          - number of output channels / same as original input channels
	nf              - number of covnolutional masks/filters
	extra_layers    - number of additional conv blocks (conv2d + batch norm + relu)
"""
function ConvDecoder(isize::Int, in_ch::Int, out_ch::Int, nf::Int, extra_layers::Int)
	layers = []

	cnf, tisize = div(nf,2), 4
	while tisize != isize
		cnf = cnf*2
		tisize = tisize*2
	end
	push!(layers, ConvTranspose((4, 4), in_ch => cnf; stride = 1, pad = 0))
	push!(layers, BatchNorm(cnf, relu))

	csize, _ = 4, cnf
	while csize < div(isize,2)
		push!(layers, ConvTranspose((4, 4), cnf => div(cnf,2); stride = 2, pad = 1))
		push!(layers, BatchNorm(div(cnf,2), relu))
		#push!(layers, x->leakyrelu.(x, 0.2f0))
		cnf = div(cnf,2)
		csize=csize*2
	end
	for i=1:extra_layers
		push!(layers, ConvTranspose((3, 3), cnf => cnf; stride = 1, pad = 1))
		push!(layers, BatchNorm(cnf, relu))
	end
	push!(layers, ConvTranspose((4, 4), cnf => out_ch, tanh; stride = 2, pad = 1))
	return Decoder(Chain(layers...))
end

Flux.@functor Decoder

""" Definition of Decoder's forward pass """
function (d::Decoder)(x)
	d.layers(x)
end

struct Discriminator
	features
	classifier
end

"""
	function Discriminator(isize::Int, in_ch::Int, out_ch::Int, nf::Int, extra_layers::Int)

Create the convolutional discriminator with parameters (almost the same structure as Encoder)
	isize           - size of output image (must be divisible by 16 i.e isize%16==0)
	in_ch           - number of input channels / dimension of latent space
	out_ch          - number of output channels / same as original input channels
	nf              - number of covnolutional masks/filters
	extra_layers    - number of additional conv blocks (conv2d + batch norm + leaky relu)
"""
function ConvDiscriminator(isize::Int, in_ch::Int, out_ch::Int, nf::Int, extra_layers::Int)
	features = []
	push!(features, Conv((4, 4), in_ch => nf; stride = 2, pad = 1))
	push!(features, BatchNorm(nf))
	push!(features, x->leakyrelu.(x, 0.2f0))

	csize, cnf = isize/2, nf
	for i=1:extra_layers
		push!(features, Conv((3, 3), cnf => cnf; stride = 1, pad = 1))
		push!(features, BatchNorm(cnf))
		push!(features, x->leakyrelu.(x, 0.2f0))
	end

	while csize > 4
		in_feat = cnf
		out_feat = cnf*2
		push!(features, Conv((4, 4), in_feat => out_feat; stride = 2, pad = 1))
		push!(features, BatchNorm(out_feat))
		push!(features, x->leakyrelu.(x, 0.2f0))
		cnf = cnf*2
		csize=csize/2
	end
	cls = Conv((4, 4), cnf => out_ch, sigmoid; stride = 1, pad = 0)
	return Discriminator(Chain(features...), cls)
end

Flux.@functor Discriminator

""" Definition of Discriminator's forward pass """
function (d::Discriminator)(x)
	feat = d.features(x)
	class = d.classifier(feat)
	return class, feat
end

"""
	Implementation of Generator, the main core of GANomaly.
	Generator consists of two encoders (E1, E2) and one decoder (D)
"""
struct Generator
	encoder1
	decoder
	encoder2
end

"""
	function Generator(isize::Int, latent_dim::Int, in_ch::Int, num_filters::Int, extra_layers::Int)

Create the convolutional discriminator with parameters (almost the same structure as Encoder)
	isize           - size of output image (must be divisible by 16 i.e isize%16==0)
	latent_dim      - dimension of latent space
	in_ch           - number of input (output) channels
	num_filters     - number of covnolutional masks/filters
	extra_layers    - number of additional conv blocks (conv2d + batch norm + leaky relu (resp. relu))
"""
function ConvGenerator(isize, latent_dim, in_ch, num_filters, extra_layers)
	encoder1 = ConvEncoder(isize, in_ch, latent_dim, num_filters, extra_layers)
	decoder = ConvDecoder(isize, latent_dim, in_ch, num_filters, extra_layers)
	encoder2 = ConvEncoder(isize, in_ch, latent_dim, num_filters, extra_layers)
	return Generator(encoder1, decoder, encoder2)
end

Flux.@functor Generator

""" Definition of Generators's forward pass """
function (g::Generator)(x)
	latent_i = g.encoder1(x)
	gen_imag = g.decoder(latent_i)
	latent_o = g.encoder2(gen_imag)
	return gen_imag, latent_i, latent_o
end

#    Implementation of loss functions for both generator and discriminator
"""
	function generator_loss(g::Generator, d::Discriminator, real_input, weights=[1,50,1])

will perform forward pass of generator and returns NTuple{4, Float32} of losses.
	Generator loss:     L_gen = w_1 * L_adv + w_2 * L_con + w_3 * L_enc
	Adversarial loss:   L_adv = || f(x) - f(D(E1(x)) ||_2
	Contextual loss:    L_con = || x - D(E1(x)) ||_1
	Encoder loss:       L_enc = || E1(x) - E2(D(E1(x))) ||_2

"""
function generator_loss(g::Generator, d::Discriminator, real_input; weights=[1,50,1])
	fake, latent_i, latent_o = g(real_input)
	pred_real, feat_real = d(real_input)
	pred_fake, feat_fake = d(fake)

	adversarial_loss = Flux.mse(feat_real, feat_fake) # l2_loss
	contextual_loss = Flux.mae(fake, real_input) # l1_loss
	encoder_loss = Flux.mse(latent_o, latent_i)
	return adversarial_loss*weights[1]+contextual_loss*weights[2]+encoder_loss*weights[3],
		adversarial_loss, contextual_loss, encoder_loss
end

"""
	function discriminator_loss(g::Generator, d::Discriminator, real_input)
		&
	function discriminator_loss(d::Discriminator, real_input, fake_input)

are equivalent version functions which perform forward pass of Discriminator
and computes discriminator loss.
"""
function discriminator_loss(g::Generator, d::Discriminator, real_input)
	fake, latent_i, latent_o = g(real_input)
	pred_real, feat_real = d(real_input)
	pred_fake, feat_fake = d(fake)

	loss_for_real = Flux.crossentropy(pred_real, 1f0)
	loss_for_fake = Flux.crossentropy(1f0.-pred_fake, 1f0)
	return 0.5f0*(loss_for_real+loss_for_fake)
end

function discriminator_loss(d::Discriminator, real_input, fake_input)
	pred_real, feat_real = d(real_input)
	pred_fake, feat_fake = d(fake_input)

	loss_for_real = Flux.crossentropy(pred_real, 1f0) # ones(typeof(pred_real), size(pred_real)) has same speed
	loss_for_fake = Flux.crossentropy(1f0.-pred_fake, 1f0)
	return 0.5f0*(loss_for_real+loss_for_fake)
end

"""
	function validation_loss(g::Generator, d::Discriminator, real_input, weights=[1,50,1])
computes generator and discriminator loss without additional pass through model
"""
function validation_loss(g::Generator, d::Discriminator, real_input; weights=[1,50,1], to_testmode::Bool=true)
	(to_testmode == true) ? begin
			Flux.testmode!(g)
			Flux.testmode!(d)
		end : nothing

	fake, latent_i, latent_o = g(real_input)
	pred_real, feat_real = d(real_input)
	pred_fake, feat_fake = d(fake)

	adversarial_loss = Flux.mse(feat_real, feat_fake) # l2_loss
	contextual_loss = Flux.mae(fake, real_input) # l1_loss
	encoder_loss = Flux.mse(latent_o, latent_i)

	loss_for_real = Flux.crossentropy(pred_real, 1f0) # ones(typeof(pred_real), size(pred_real)) has same speed
	loss_for_fake = Flux.crossentropy(1f0.-pred_fake, 1f0)

	(to_testmode == true) ? begin
			Flux.testmode!(g, false)
			Flux.testmode!(d, false)
		end : nothing

	return adversarial_loss*weights[1]+contextual_loss*weights[2]+encoder_loss*weights[3],
		0.5f0*(loss_for_real+loss_for_fake)
end

"""
	function anomaly_score(generator::Generator, real_input;dims=3)

computes unscaled anomaly score A(x) = || E1(x) - E2(D(E1(x))) ||_1
"""
function anomaly_score(generator::Generator, real_input; dims=3, to_testmode::Bool=true)
	(to_testmode == true) ? Flux.testmode!(generator) : nothing
	_, latent_i, latent_o = generator(real_input)
	(to_testmode == true) ? Flux.testmode!(generator, false) : nothing
	return vec(Flux.mae(latent_i, latent_o, agg=x->mean(x, dims=dims)))'
end

function anomaly_score_gpu(generator::Generator, real_input; dims=3, batch_size=64, to_testmode::Bool=true)
	real_input = Flux.Data.DataLoader(real_input, batchsize=batch_size)
	(to_testmode == true) ? Flux.testmode!(generator) : nothing
	generator = generator|>gpu
	output = Array{Float32}([])
	for X in real_input
		_, latent_i, latent_o = generator(X|>gpu)
		output = cat(output,vec(Flux.mae(latent_i, latent_o, agg=x->mean(x, dims=dims))) |> cpu, dims=1)
	end
	(to_testmode == true) ? Flux.testmode!(generator, false) : nothing
	return output'
end

"""
	Model's training and inference
"""

"""
	fit!(generator::Generator, discriminator::Discriminator, data, params)
"""
function StatsBase.fit!(generator::Generator, discriminator::Discriminator, data, params)
	# prepare batches & loaders
	train_loader, val_loader = prepare_dataloaders(data, params)
	# training info logger
	history = GANomalyHistory()
	# prepare for early stopping
	best_generator = deepcopy(generator)
	best_discriminator = deepcopy(discriminator)
	patience = params.patience
	best_val_loss = 1e10
	val_batches = length(val_loader)

	# ADAMW(η = 0.001, β = (0.9, 0.999), decay = 0) = Optimiser(ADAM(η, β), WeightDecay(decay))
	opt_g = haskey(params, :decay) ? ADAMW(params.lr, (0.9, 0.999), params.decay) : ADAM(params.lr)
	opt_d = haskey(params, :decay) ? ADAMW(params.lr, (0.9, 0.999), params.decay) : ADAM(params.lr)

	ps_g = Flux.params(generator)
	ps_d = Flux.params(discriminator)

	loss_weights = haskey(params, :weights)  ? params.weights : [1, 50, 1]

	progress = Progress(length(train_loader))
	for (iter, X) in enumerate(train_loader)
		#generator update
		loss1, back = Flux.pullback(ps_g) do
			generator_loss(generator, discriminator, getobs(X)|>gpu, weights=loss_weights)
		end
		grad = back((1f0, 0f0, 0f0, 0f0))
		Flux.Optimise.update!(opt_g, ps_g, grad)

		# discriminator update
		loss2, back = Flux.pullback(ps_d) do
			discriminator_loss(generator, discriminator, getobs(X)|>gpu)
		end
		grad = back(1f0)
		Flux.Optimise.update!(opt_d, ps_d, grad)

		history = update_history(history, loss1, loss2)
		next!(progress; showvalues=[
			(:iters, "$(iter)/$(params.iters)"),
			(:generator_loss, loss1[1]),
			(:discriminator_loss, loss2)
			])
		#TODO optionaly add discriminator restrart if its loss drops under 1e-5
		if mod(iter, params.check_every) == 0
			total_val_loss_g = 0
			total_val_loss_d = 0
			for X_val in val_loader
				vgl, vdl = validation_loss(generator, discriminator, X_val |> gpu)
				total_val_loss_g += vgl
				total_val_loss_d += vdl
			end
			history = update_val_history(history, total_val_loss_g/val_batches, total_val_loss_d/val_batches)
			if total_val_loss_g < best_val_loss
				best_val_loss = total_val_loss_g
				patience = params.patience
				best_generator = deepcopy(generator)
				best_discriminator = deepcopy(discriminator)
			else
				patience -= 1
				if patience == 0
					@info "Stopped training after $(iter) iterations"
					global iters = iter - params.check_every*params.patience
					break
				end
			end
		end
		if iter == params.iters
			global iters = params.iters
		end
	end
	return history, best_generator, best_discriminator, sum(map(p->length(p), ps_g))+sum(map(p->length(p), ps_d)), iters
end

"""
	StatsBase.predict(generator::Generator, data; dims=3)

computes scaled anomaly score (interval [0,1]). Data should contain both anomal and normal datapoints.
"""

function StatsBase.predict(generator::Generator, data; dims=3)
	_, latent_i, latent_o = generator(data)
	anomaly_score = anomaly_score(generator, data; dims=dims)
	s_anomaly_score = (anomaly_score .- minimum(anomaly_score))./(maximum(anomaly_score)-minimum(anomaly_score))
	return s_anomaly_score
end

"""
	conv GANomaly constructor for parameters
"""
function ganomaly_constructor(kwargs)
	generator_params = (isize=kwargs.isize,
						latent_dim=kwargs.hdim,
						in_ch = kwargs.in_ch,
						nf = kwargs.num_filters,
						extra_layers = kwargs.extra_layers)

	discriminator_params = (isize=kwargs.isize,
							in_ch = kwargs.in_ch,
							out_ch = 1,
							nf = kwargs.num_filters,
							extra_layers = kwargs.extra_layers)

	# little control to random initialization
	(kwargs.init_seed !== nothing) ? Random.seed!(kwargs.init_seed) : nothing
	generator = ConvGenerator(generator_params...)
	discriminator = ConvDiscriminator(discriminator_params...)
	(kwargs.init_seed !== nothing) ? Random.seed!() : nothing

	return generator, discriminator, generator_params, discriminator_params
end

function tabular_ganomaly_constructor(kwargs)
	(kwargs.init_seed !== nothing) ? Random.seed!(kwargs.init_seed) : nothing
	encoder1 = build_mlp(
		kwargs.idim,
		kwargs.hdim,
		kwargs.zdim,
		kwargs.nlayers,
		activation=kwargs.activation,
		lastlayer="linear"
	)
	encoder2 = build_mlp(
		kwargs.idim,
		kwargs.hdim,
		kwargs.zdim,
		kwargs.nlayers,
		activation=kwargs.activation,
		lastlayer="linear"
	)
	decoder = build_mlp(
		kwargs.zdim,
		kwargs.hdim,
		kwargs.idim,
		kwargs.nlayers,
		activation=kwargs.activation,
		lastlayer="tanh"
	)

	features = build_mlp(
		kwargs.idim,
		kwargs.hdim,
		kwargs.hdim,
		Int(kwargs.nlayers - 1),
		activation=kwargs.activation
	)

	classsifier = Flux.Dense(kwargs.hdim, 1, σ)

	generator = Generator(encoder1, decoder, encoder2)
	discriminator = Discriminator(features, classsifier)
	(kwargs.init_seed !== nothing) ? Random.seed!() : nothing

	return generator, discriminator
end


text:
function displaz_version()
    displaz_cmd = get(ENV, "DISPLAZ_CMD", "displaz")
    try
        verstring = read(`$displaz_cmd -version`, String)
        m = match(r"version *([\d\.]*)-?", verstring)
        VersionNumber(m.captures[1])
    catch err
        @warn(err)
        nothing
    end
end

ver = displaz_version()

if ver === nothing
    @warn(
        """
        displaz could not be found - please check that you have it installed
        (for now, you will need to build it from source - see
        https://github.com/c42f/displaz).  To ensure julia can find it, you
        should add its location to your PATH, or optionally set the DISPLAZ_CMD
        environment variable to the full path of the executable.
        """
    )
elseif ver < v"0.3.1-317"
    @warn(
        """
        You have an older version ($ver) of displaz.  You should upgrade
        to the latest version to ensure the Displaz.jl bindings work correctly.
        """
    )
end


text:
#!/usr/bin/julia

include("vis.jl")

img_name = ARGS[1]
out_name = ARGS[2]

braille_chars = image_to_braille(img_name)

io = open(out_name,"w")

for word in braille_chars
	for char in word
		for c in transpose(char)
			write(io, string(c))
			write(io, ",")
		end
		write(io,"\n")
	end
end

close(io)


text:
# File   : spherical_cylinder_approx.jl
# License: GNU v3.0
# Author : Andrei Leonard Nicusan <a.l.nicusan@bham.ac.uk>
# Date   : 28.08.2021

using LinearAlgebra

using GLMakie
using GeometryBasics

# Vertices and face indices generated using MeshLab; two spherical caps were created, then the
# convex hull was saved in the "spherical_cylinder.ply" file
verts = [
    2.468838e-09 -8.270091e-08 -1
    2.468838e-09 -4.721871e-09 2
    0.258819 0.9659259 1
    -0.9659259 0.258819 1
    0.9659259 -0.258819 1
    -0.4062855 0.703707 -0.5828626
    -0.4062854 0.703707 1.582862
    -0.258819 -0.9659259 1
    0.8660254 0.5 -4.371139e-08
    0.6957635 0.2464413 1.67467
    -0.6957635 -0.2464414 -0.6746702
    -0.8125709 -2.265369e-09 1.582862
    -0.7071068 0.7071068 1.000345
    0.4062854 -0.703707 1.582862
    0.4062855 -0.7037071 -0.5828624
    0.258819 0.9659259 -8.444392e-08
    0.7071068 0.7071068 1.000345
    0.2761805 0.4783586 1.833605
    -0.4167868 0.2406319 -0.8765758
    -0.9659259 0.258819 -2.262668e-08
    -0.4167868 0.240632 1.876576
    -0.8660254 -0.5 1
    -0.258819 0.9659259 1
    -0.6957635 0.2464413 1.67467
    0.4167868 -0.240632 1.876576
    -0.4062855 -0.703707 1.582862
    -0.2761805 -0.4783587 -0.8336051
    0.9659259 -0.258819 2.262668e-08
    0.258819 -0.9659259 1
    0.5613061 0.4793281 -0.6746704
    0.9659259 0.258819 1
    0.5 0.8660254 1
    0.1344573 0.7257695 1.67467
    -0.7071068 0.7071068 -0.0003454111
    -0.6957635 0.2464413 -0.6746703
    -0.1385453 0.2399675 1.960844
    -0.7071068 -0.7071068 -0.0003452875
    -0.9659259 -0.258819 1
    -0.5 0.8660254 1
    0 1 1
    -0.5613061 0.4793282 1.67467
    -0.8660254 0.5 1
    -0.552361 2.718866e-10 1.833605
    0.6957635 -0.2464413 1.67467
    0.2770906 -1.700542e-10 1.960844
    -0.1344573 -0.7257695 1.67467
    -0.7071068 -0.7071068 1.000345
    -0.6957635 -0.2464413 1.67467
    0.2761805 -0.4783587 -0.8336051
    -0.2770906 -8.282016e-08 -0.9608439
    0 -1 8.742278e-08
    0.7071068 -0.7071068 1.000345
    0.7071068 0.7071068 -0.0003454111
    0.812571 -4.863583e-08 -0.5828623
    0.1344573 0.7257694 -0.6746703
    0.8660254 0.5 1
    0.9659259 0.258819 -2.262668e-08
    0.4062855 0.703707 1.582862
    -0.258819 0.9659259 -8.444392e-08
    -0.552361 -7.314797e-08 -0.8336051
    -0.8125709 -4.869008e-08 -0.5828624
    -0.5613062 0.4793281 -0.6746703
    -0.8660254 0.5 -4.371139e-08
    4.927248e-09 0.481264 1.876576
    -0.4062854 -0.7037071 -0.5828624
    -1 1.224647e-16 1
    -0.9659259 -0.258819 2.262668e-08
    -0.1344573 0.7257695 1.67467
    -0.2770906 -1.17948e-09 1.960844
    0.5613061 -0.4793282 1.67467
    0.552361 -2.306788e-10 1.833605
    0.812571 -2.31961e-09 1.582862
    0.1385453 -0.2399675 1.960844
    0.1385453 0.2399675 1.960844
    0.1344573 -0.7257695 1.67467
    -0.5613062 -0.4793282 1.67467
    -0.5 -0.8660254 1
    0.1385453 -0.2399676 -0.9608439
    0.1344573 -0.7257695 -0.6746702
    -0.1385453 0.2399674 -0.9608439
    0 -1 1
    0.5 -0.8660254 7.571035e-08
    -0.258819 -0.9659259 8.444392e-08
    0.5 -0.8660254 1
    0.8660254 -0.5 1
    0.5 0.8660254 -7.571035e-08
    0.6957635 0.2464413 -0.6746702
    0.6957635 -0.2464414 -0.6746702
    1 0 1
    1 0 0
    0.5613062 0.4793282 1.67467
    -0.5 0.8660254 -7.571035e-08
    0 1 -8.742278e-08
    -0.4167868 -0.2406321 -0.8765758
    -0.2761805 0.4783585 -0.8336052
    -0.2761805 0.4783586 1.833605
    -0.5613061 -0.4793282 -0.6746703
    -0.8660254 -0.5 4.371139e-08
    -1 -1.224647e-16 1.07062e-23
    0.4167869 0.240632 1.876576
    -0.1385453 -0.2399675 1.960844
    0.2761805 -0.4783586 1.833605
    -0.4167868 -0.240632 1.876576
    -0.2761805 -0.4783586 1.833605
    4.927248e-09 -0.4812641 -0.8765757
    -0.1385453 -0.2399676 -0.9608439
    -0.1344573 -0.7257695 -0.6746702
    0.258819 -0.9659259 8.444392e-08
    0.2761805 0.4783586 -0.8336052
    0.7071068 -0.7071068 -0.0003452875
    -0.5 -0.8660254 7.571035e-08
    0.4062854 0.703707 -0.5828626
    0.2770906 -8.382958e-08 -0.9608439
    0.5613062 -0.4793282 -0.6746702
    0.552361 -7.26454e-08 -0.8336051
    -0.1344573 0.7257694 -0.6746703
    -4.927248e-09 -0.481264 1.876576
    -4.927248e-09 0.4812639 -0.8765758
    0.8660254 -0.5 4.371139e-08
    0.1385453 0.2399674 -0.9608439
    0.4167868 0.2406319 -0.8765758
    0.4167869 -0.2406321 -0.8765757
]

faces = [
    23 3 11
    38 12 6
    38 6 22
    39 32 2
    40 6 12
    40 23 20
    41 3 23
    41 40 12
    41 23 40
    42 23 11
    42 20 23
    45 25 7
    47 42 11
    47 37 21
    47 11 37
    52 8 29
    55 9 30
    55 52 16
    55 8 52
    56 55 30
    56 8 55
    57 31 2
    57 2 32
    57 16 31
    57 32 17
    59 34 18
    59 18 49
    60 19 34
    60 59 10
    60 34 59
    61 18 34
    61 33 5
    62 34 19
    62 61 34
    62 33 61
    62 41 12
    62 12 33
    62 19 3
    62 3 41
    63 17 32
    65 37 11
    65 11 3
    65 3 19
    66 60 10
    67 39 22
    67 22 6
    67 32 39
    67 63 32
    68 35 20
    68 20 42
    68 1 35
    69 43 24
    69 13 51
    70 24 43
    70 44 24
    71 30 9
    71 43 4
    71 70 43
    71 9 70
    72 44 1
    72 24 44
    73 35 1
    73 1 44
    73 63 35
    73 17 63
    74 28 13
    75 46 25
    75 21 46
    75 47 21
    76 7 25
    76 25 46
    78 48 14
    79 49 18
    79 0 49
    80 28 74
    80 74 45
    80 45 7
    82 80 7
    82 50 80
    82 7 76
    83 51 13
    83 13 28
    84 4 43
    84 69 51
    84 43 69
    85 15 2
    85 2 31
    85 31 16
    85 16 52
    86 29 8
    86 56 53
    86 8 56
    87 53 27
    88 71 4
    88 30 71
    89 27 53
    89 53 56
    89 88 4
    89 4 27
    89 56 30
    89 30 88
    90 57 17
    90 9 55
    90 55 16
    90 16 57
    91 58 5
    91 5 33
    91 33 12
    91 12 38
    91 38 22
    91 22 58
    92 58 22
    92 22 39
    92 39 2
    92 2 15
    92 15 54
    93 59 49
    93 10 59
    94 61 5
    94 18 61
    94 79 18
    95 20 35
    95 35 63
    95 40 20
    95 63 67
    95 67 6
    95 6 40
    96 64 36
    96 26 64
    96 93 26
    96 10 93
    97 66 10
    97 46 21
    97 36 46
    97 21 37
    97 37 66
    97 96 36
    97 10 96
    98 65 19
    98 66 37
    98 37 65
    98 19 60
    98 60 66
    99 70 9
    99 90 17
    99 9 90
    99 17 73
    99 73 44
    99 44 70
    100 72 1
    100 1 68
    101 24 72
    101 69 24
    101 74 13
    101 13 69
    102 42 47
    102 47 75
    102 100 68
    102 68 42
    103 75 25
    103 25 45
    103 102 75
    103 100 102
    104 77 48
    104 48 78
    105 49 0
    105 0 77
    105 93 49
    105 77 104
    105 104 26
    105 26 93
    106 64 26
    106 78 50
    106 104 78
    106 26 104
    106 82 64
    106 50 82
    107 78 14
    107 50 78
    107 14 81
    107 83 28
    107 81 83
    107 80 50
    107 28 80
    109 81 14
    109 83 81
    109 51 83
    110 82 76
    110 76 46
    110 46 36
    110 36 64
    110 64 82
    111 54 15
    111 15 85
    111 85 52
    111 52 29
    111 108 54
    111 29 108
    112 77 0
    113 14 48
    113 109 14
    114 86 53
    114 53 87
    115 92 54
    115 94 5
    115 5 58
    115 58 92
    116 101 72
    116 72 100
    116 100 103
    116 103 45
    116 45 74
    116 74 101
    117 79 94
    117 94 115
    117 115 54
    117 54 108
    118 87 27
    118 113 87
    118 109 113
    118 27 4
    118 4 84
    118 84 51
    118 51 109
    119 112 0
    119 0 79
    119 79 117
    119 117 108
    120 114 112
    120 119 108
    120 112 119
    120 108 29
    120 29 86
    120 86 114
    121 113 48
    121 114 87
    121 87 113
    121 48 77
    121 77 112
    121 112 114
]

# Split vertices into the bottom and top caps
half = size(verts)[1] ÷ 2

verts_bottom = zeros(half, 3)
index_bottom = 1

verts_top = zeros(half, 3)
index_top = 1

faces_split = copy(faces)

for (i, v) in enumerate(eachrow(verts))
    if v[3] < 0.5
        verts_bottom[index_bottom, :] .= v
        faces_split[faces .== i - 1] .= index_bottom - 1

        global index_bottom += 1
    else
        verts_top[index_top, :] .= v
        faces_split[faces .== i - 1] .= index_top - 1 + half

        global index_top += 1
    end
end

verts_split = [verts_bottom ; verts_top]

# Move top cap back to (0, 0, 0) and flip it
verts_split[half + 1:end, 3] .-= 1
verts_split[half + 1:end, 3] .*= -1

# Create spherical cylinder spanning two points p1 and p2 on a particle trajectory
p1 = [5, 4, 3]
r1 = 2

p2 = [6, 7, 8]
r2 = 4

# Expand vertices to given particles' radii
verts_split[1:half, :] .*= r1
verts_split[half + 1:end, :] .*= r2

# Create rotation matrix from vertical +Z to |p2 - p1|
a = [0., 0., 1.]
b = [0., 0., 0.]

# Compute b = |p2 - p1|
s = 0.
for i in 1:3
    b[i] = p2[i] - p1[i]
    global s += b[i] * b[i]
end

# Divide by ||p2 - p1|| to get unity vector
s = sqrt(s)
for i in 1:3
    b[i] /= s
end

# Compute cross product between a and b
v = [0., 0., 0.]
v[1] = a[2] * b[3] - a[3] * b[2]
v[2] = a[3] * b[1] - a[1] * b[3]
v[3] = a[1] * b[2] - a[2] * b[1]

# Compute dot product between a and b
c = a[1] * b[1] + a[2] * b[2] + a[3] * b[3]

# Create skew-symmetric cross-product matrix of v
vx = [0., 0., 0., 0., 0., 0., 0., 0., 0.]
vx[2] = -v[3]
vx[3] = v[2]
vx[4] = v[3]
vx[6] = -v[1]
vx[7] = -v[2]
vx[8] = v[1]

# Multiply vx by vx
vx2 = [0., 0., 0., 0., 0., 0., 0., 0., 0.]
for i in 1:3
    for j in 1:3
        for k in 1:3
            vx2[(i - 1) * 3 + j] += vx[(i - 1) * 3 + k] * vx[(k - 1) * 3 + j]
        end
    end
end

# Create rotation matrix as I(3) + vx + vx.vx / (1 + c)
rot = [1., 0., 0., 0., 1., 0., 0., 0., 1.]

if abs(1 + c) > 1e-6
    for i in 1:9
        rot[i] += vx[i] + vx2[i] / (1 + c)
    end
end

# Rotate all vertices
temp = [0., 0., 0.]
for i in 1:size(verts_split)[1]
    for j in 1:3
        temp[j] = verts_split[i, j]
    end

    for j in 1:3
        verts_split[i, j] = rot[(j - 1) * 3 + 1] * temp[1] +
                            rot[(j - 1) * 3 + 2] * temp[2] +
                            rot[(j - 1) * 3 + 3] * temp[3]
    end
end

# Translate vertices
verts_split[1:half, :] .+= p1'
verts_split[half + 1:end, :] .+= p2'

# Plot mesh and direction from p1 to p2
m = GeometryBasics.Mesh(
    Point3.(eachrow(verts_split)),
    TriangleFace.(eachrow(faces_split .+ 1)),
)

wireframe(m)
lines!([p1 p2], color = :red)

current_figure()

max_stars_repo_name: anicusan/KonigCell

max_issues_repo_name: anicusan/KonigCell

max_forks_repo_name: anicusan/KonigCell


text:
using LinearAlgebra
using LatticeDiracOperators

function test_staggered()
    NX=4
    NY=4
    NZ=4
    NT=4
    Nwing = 1
    NC=3

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    U2 = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "hot")
    x = Initialize_pseudofermion_fields(U[1],"staggered")

    gauss_distribution_fermion!(x)
    #x[1,1,1,1,1,1] = 4
    println(x[1,1,1,1,1,1])

    params = Dict()
    params["Dirac_operator"] = "staggered"
    params["mass"] = 0.1
    params["verbose_level"] = 3
    D = Dirac_operator(U,x,params)
    D2 = D(U2)

    for Nf in [8,4,2]
        println("Nf = $Nf")
        parameters_action = Dict()
        parameters_action["Nf"] = Nf
        fermi_action = FermiAction(D,parameters_action)
        gauss_sampling_in_action!(x,U,fermi_action)
        println("Sfold = ", dot(x,x))
        y = similar(x)
        sample_pseudofermions!(y,U,fermi_action,x)

        UdSfdU = calc_UdSfdU(fermi_action,U,y)

        Sf = evaluate_FermiAction(fermi_action,U,y)
        println("Sfnew = ", Sf)

    end

    y = similar(x)
    mul!(y,D,x)
    mul!(y,D2,x)
    mul!(y,D(U2),x)

    println("BICG method")
    @time solve_DinvX!(y,D,x)
    #@time bicg(y,D,x,verbose = Verbose_3())

    DdagD = DdagD_operator(U,x,params)
    mul!(y,DdagD,x)

    @time solve_DinvX!(y,DdagD,x)

    #@time cg(y,DdagD,x,verbose = Verbose_3())

end

function test_wilson()
    NX=4
    NY=4
    NZ=4
    NT=4
    Nwing = 1
    NC=3

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    x = Initialize_pseudofermion_fields(U[1],"Wilson")
    gauss_distribution_fermion!(x)
    println(x[1,1,1,1,1,1])

    params = Dict()
    params["Dirac_operator"] = "Wilson"
    params["κ"] = 0.1
    params["verbose_level"] = 3
    D = Dirac_operator(U,x,params)
    D2 = D(U)

    y = similar(x)

    for i=1:3
        y = similar(x)
        y1 = deepcopy(y)
        y2 = deepcopy(y)

        println("BICG")
        @time bicg(y1,D,x, verbose = Verbose_print(3))
        println(dot(y1,y1))
        println("BICGstab")
        @time bicgstab(y2,D,x,verbose = Verbose_print(3))
        println(dot(y2,y2))
    end

    parameters_action = Dict()
    fermi_action = FermiAction(D,parameters_action)
    gauss_sampling_in_action!(x,U,fermi_action)
    println("Sfold = ", dot(x,x))
    y = similar(x)
    sample_pseudofermions!(y,U,fermi_action,x)

    UdSfdU = calc_UdSfdU(fermi_action,U,y)

    Sf = evaluate_FermiAction(fermi_action,U,y)
    println("Sfnew = ", Sf)

    #=

    y = similar(x)
    mul!(y,D,x)

    println("BICG method")
    @time solve_DinvX!(y,D,x)

    #@time bicg(y,D,x,verbose = Verbose_3())

    DdagD = DdagD_operator(U,x,params)
    mul!(y,DdagD,x)
    @time solve_DinvX!(y,DdagD,x)
    #@time cg(y,DdagD,x,verbose = Verbose_3())
    =#

    return

end
#test_staggered()
test_wilson()


text:
using StatsBase
using JLD2

daily(times::AbstractArray{T} where T, max_days::Integer) = fit(Histogram, times, 0:(max_days) ).weights
daily(times::AbstractVector{T} where T) = daily(times, maximum(times))

mutable struct DailyTrajectories <: Output
  file::JLD2.JLDFile
  last_trajectory::Int
end

DailyTrajectories(fname::AbstractString, ::Integer) = DailyTrajectories(JLD2.jldopen(fname, "w+", compress=true), 0)

function pushtrajectory!(d::DailyTrajectories, state::Simulation.SimState, params::Simulation.SimParams, cb::DetectionCallback)
  d.last_trajectory += 1
  trajectory_group = JLD2.Group(d.file, string(d.last_trajectory))
  save_daily_trajectories(trajectory_group, state, params, cb)
end

aftertrajectories(d::DailyTrajectories) = close(d.file)

function save_daily_trajectories(dict, state::Simulation.SimState, params::Simulation.SimParams, cb::DetectionCallback)
  max_days = Simulation.time(state) |> floor |> Int
  num_individuals = Simulation.numindividuals(state)

  infection_times = Vector{OptTimePoint}(missing, num_individuals)
  contact_kinds = Vector{Simulation.ContactKind}(undef, num_individuals)

  for i in 1:num_individuals
    event = Simulation.backwardinfection(state, i)
    kind = contactkind(event)

    contact_kinds[i] = kind
    infection_times[i] = ifelse(kind == Simulation.NoContact, time(event), missing)
  end
  hospitalization_progressions = getproperty.(params.progressions, :severe_symptoms_time)
  death_progressions = getproperty.(params.progressions, :death_time)

  dict["infection_times"] = daily(filter(!ismissing,infection_times), max_days)
  dict["detection_tiems"] = daily(filter(!ismissing, cb.detection_times), max_days)
  dict["daily_deaths"] = daily(filter(!ismissing, infection_times.+death_progressions), max_days)
  dict["daily_hospitalizations"] = daily(filter(!ismissing, infection_times.+hospitalization_progressions), max_days)
  for kind in instances(ContactKind)
    if kind != NoContact
      dict["daily_" * lowercase(string(kind))] = daily(infection_times[contact_kinds.==Int(kind)], max_days)
    end
  end
end


text:
using Test
using POMDPs
using POMDPGym
using POMDPSimulators, POMDPPolicies

# Test MDP construction
mdp = GymPOMDP(:CartPole, pixel_observations = true)

@test mdp.pixel_observations == true
@test mdp.γ == 0.99
@test mdp.actions == [1,2]

# Test basic POMDP interface
@test length(rand(initialstate(mdp))) == 4
@test actions(mdp) == [1,2]
@test actionindex(mdp, 1) == 1
@test actionindex(mdp, 2) == 2
@test isterminal(mdp, rand(initialstate(mdp))) == false
@test discount(mdp) == mdp.γ

# Test gen function
render(mdp.env)

sp, o, r = gen(mdp, mdp.env.state, 1)
close(mdp.env)
@test length(sp) == 4
@test r isa Real
@test size(o) == (400,600,3)

mdp_nopix = GymPOMDP(:CartPole, pixel_observations = false)
sp, o, r = gen(mdp_nopix, mdp_nopix.env.state, 1)
@test all(o .≈ sp)
render(mdp_nopix)
close(mdp_nopix.env)

# Do 1 episode
h = simulate(HistoryRecorder(), mdp_nopix, RandomPolicy(mdp_nopix))
@test isterminal(mdp_nopix, mdp_nopix.env.state)

# Test out prescribed actions
mdp = GymPOMDP(:Pendulum, pixel_observations = true, actions = [[-1.], [1.]])
@test mdp.actions == [[-1.], [1.]]
h = simulate(HistoryRecorder(max_steps = 200), mdp, RandomPolicy(mdp))
close(mdp.env)
@test !isterminal(mdp, mdp.env.state)

# Test prescribed render function
my_render(sp) = ones(20,20)
mdp = GymPOMDP(:Pendulum, pixel_observations = true, actions = [[-1.], [1.]], special_render = my_render)
sp, o, r = gen(mdp, mdp.env.state, [1.])
close(mdp.env)
@test size(o) == (20,20)

max_forks_repo_name: ancorso/POMDPGym


text:
using YAML

function WeaveDoc(source, informat = nothing)
    path = abspath(source)
    _, fname = splitdir(path)
    basename = splitext(fname)[1]

    isnothing(informat) && (informat = detect_informat(source))
    header, chunks = parse_doc(read(source, String), informat)

    # update default chunk options from header
    chunk_defaults = deepcopy(get_chunk_defaults())
    if (weave_options = get(header, WEAVE_OPTION_NAME, nothing)) !== nothing
        for key in keys(chunk_defaults)
            if (val = get(weave_options, string(key), nothing)) !== nothing
                chunk_defaults[key] = val
            end
        end
    end
    if haskey(header, WEAVE_OPTION_NAME_DEPRECATED)
        @warn "Weave: `options` key is deprecated. Use `weave_options` key instead." _id = WEAVE_OPTION_DEPRECATE_ID maxlog = 1
        for key in keys(chunk_defaults)
            if (val = get(header[WEAVE_OPTION_NAME_DEPRECATED], string(key), nothing)) !== nothing
                chunk_defaults[key] = val
            end
        end
    end

    return WeaveDoc(
        source,
        basename,
        path,
        chunks,
        "",
        nothing,
        "",
        "",
        header,
        chunk_defaults,
    )
end

"""
    detect_informat(path)

Detect Weave input format based on file extension of `path`.
"""
function detect_informat(path)
    ext = lowercase(last(splitext(path)))

    ext == ".jl" && return "script"
    ext == ".jmd" && return "markdown"
    ext == ".ipynb" && return "notebook"
    return "noweb"
end

function parse_doc(document, informat)
    document = replace(document, "\r\n" => "\n") # normalize line ending

    return informat == "markdown" ? parse_markdown(document) :
        informat == "noweb" ? parse_markdown(document; is_pandoc = true) :
        informat == "script" ? parse_script(document) :
        informat == "notebook" ? parse_notebook(document) :
        error("unsupported input format given: $informat")
end

# inline
# ------

function DocChunk(text::AbstractString, number, start_line; notebook = false)
    # don't parse inline code in notebook
    content = notebook ? parse_inline(text) : parse_inlines(text)
    return DocChunk(content, number, start_line)
end

const   INLINE_REGEX = r"`j\s+(.*?)`"
const INLINE_REGEXES = r"`j\s+(.*?)`|^!\s(.*)$"m

# handle code units correctly !
function parse_inlines(str)
    ret = Inline[]
    s = 1
    code_no = text_no = 0
    for m in eachmatch(INLINE_REGEXES, str)
        e = m.offset
        push!(ret, InlineText((str[s:prevind(str,e)]), text_no += 1))
        i = findfirst(!isnothing, m.captures)
        push!(ret, InlineCode(m.captures[i], code_no += 1, isone(i) ? :inline : :line))
        s = e + ncodeunits(m.match)
    end
    push!(ret, InlineText(str[s:end], text_no += 1))
    return ret
end

parse_inline(str) = Inline[InlineText(str, 1)]

# options
# -------

const OptionDict = Dict{Symbol,Any}

function parse_options(str)::OptionDict
    str = string('(', str, ')')
    ex = Meta.parse(str)
    nt = if Meta.isexpr(ex, (
        :block, # "(k1 = v1; k2 = v2, ...)"
        :tuple, # "(k1 = v1, k2 = v2, ...)"
    ))
        eval(Expr(:tuple, filter(is_valid_kv, ex.args)...))
    elseif is_valid_kv(ex) # "(k = v)"
        eval(Expr(:tuple, ex))
    else
        NamedTuple{}()
    end
    return dict(nt)
end

is_valid_kv(x) = Meta.isexpr(x, :(=))
dict(nt) = Dict((k => v for (k,v) in zip(keys(nt), values(nt)))...)
nt(dict) = NamedTuple{(Symbol.(keys(dict))...,)}((collect(values(dict))...,))

# each input format
# -----------------

include("markdown.jl")
include("script.jl")
include("notebook.jl")


text:
using DifferentiableStateSpaceModels
using Test, LinearAlgebra

# The BLAS threads is still an issue in Julia 1.7
# This has no effect with MKL
DifferentiableStateSpaceModels.set_blas_threads()

println("Running Testsuite with Threads.nthreads() = $(Threads.nthreads()) BLAS.vendor = $(BLAS.vendor()), and BLAS.num_threads = $(BLAS.get_num_threads()) \n")

# Delete the .function_cache
# e.g. ENV["DSSM_TEST_DELETE_CACHE"] = "false" environment variable to turn off, can be global
get(ENV, "DSSM_TEST_DELETE_CACHE", "true") == "true" &&
    rm(default_model_cache_location(); force = true, recursive = true)

include("make_perturbation_model.jl")
include("first_order_perturbation.jl")
include("second_order_perturbation.jl")
# include("first_order_sequence.jl")
# include("second_order_sequence.jl")
# include("rbc_estimation.jl")
include("sgu.jl")
include("FVGQ20.jl")
include("symbolic_utils.jl")
include("utils.jl")


text:
#=
    The file contains implementation of the F4 algorithm,
    described originally by Jean-Charles Faugere in
        "A new efficient algorithm for computing Grobner bases"
=#

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

function select_normal!(
            pairset::Pairset, basis::Basis, matrix::MacaulayMatrix,
            ht::MonomialHashtable, symbol_ht::MonomialHashtable;
            maxpairs::Int=0)

    sort_pairset_by_degree!(pairset, 1, pairset.load - 1)
    # sort by degree

    ps = pairset.pairs
    min_deg = ps[1].deg
    min_idx = 0

    while min_idx < pairset.load && ps[min_idx + 1].deg == min_deg
        min_idx += 1
    end

    # number of selected pairs
    npairs = min_idx

    sort_pairset_by_lcm!(pairset, npairs, ht)

    reinitialize_matrix!(matrix, npairs)

    uprows  = matrix.uprows
    lowrows = matrix.lowrows

    # polynomials from pairs in order (p11, p12)(p21, p21)
    # (future rows of the matrix)
    gens   = Vector{Int}(undef, 2*npairs)

    # buffer !
    etmp = ht.exponents[1]
    i = 1
    while i <= npairs
        matrix.ncols += 1
        load = 1
        lcm = ps[i].lcm
        j = i

        # we collect all generators with same lcm into gens
        while j <= npairs && ps[j].lcm == lcm
            gens[load] = ps[j].poly1
            load += 1
            gens[load] = ps[j].poly2
            load += 1
            j += 1
        end
        load -= 1

        # sort by number in the basis (by=identity)
        sort_generators_by_position!(gens, load)

        # now we collect reducers, and reduced

        # first generator index in groebner basis
        prev = gens[1]
        # first generator in hash table
        poly = basis.gens[prev]
        # first generator lead monomial index in hash data
        vidx = poly[1]

        # first generator exponent
        eidx = ht.exponents[vidx]
        # exponent of lcm corresponding to first generator
        elcm = ht.exponents[lcm]
        for u in 1:ht.explen
            etmp[u] = elcm[u] - eidx[u]
        end
        # now etmp contents complement to eidx in elcm

        # hash of complement
        htmp = ht.hashdata[lcm].hash - ht.hashdata[vidx].hash

        # add row as a reducer
        matrix.nup += 1
        uprows[matrix.nup] = multiplied_poly_to_matrix_row!(symbol_ht, ht, htmp, etmp, poly)
        # map upper row to index in basis
        matrix.up2coef[matrix.nup] = prev

        # mark lcm column as reducer in symbolic hashtable
        symbol_ht.hashdata[uprows[matrix.nup][1]].idx = 2
        # increase number of rows set
        matrix.nrows += 1

        # over all polys with same lcm,
        # add them to the lower part of matrix
        for k in 1:load
            # duplicate generator,
            # we can do so as long as generators are sorted
            if gens[k] == prev
                continue
            end

            # if the table was reallocated
            elcm = ht.exponents[lcm]

            # index in gb
            prev = gens[k]
            # poly of indices of monoms in hash table
            poly = basis.gens[prev]
            vidx = poly[1]
            # leading monom idx
            eidx = ht.exponents[vidx]
            for u in 1:ht.explen
                etmp[u] = elcm[u] - eidx[u]
            end

            htmp = ht.hashdata[lcm].hash - ht.hashdata[vidx].hash

            # add row to be reduced
            matrix.nlow += 1
            lowrows[matrix.nlow] = multiplied_poly_to_matrix_row!(symbol_ht, ht, htmp, etmp, poly)
            # map lower row to index in basis
            matrix.low2coef[matrix.nlow] = prev

            symbol_ht.hashdata[lowrows[matrix.nlow][1]].idx = 2

            matrix.nrows += 1
        end

        i = j
    end

    resize!(matrix.lowrows, matrix.nrows - matrix.ncols)

    # remove selected parirs from pairset
    for i in 1:pairset.load - npairs
        ps[i] = ps[i + npairs]
    end
    pairset.load -= npairs
end

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

# Prepares the given set of polynomials L to be reduced by the ideal G,
# starts the reduction routine,
# and returns reduced polynomials
function reduction!(
            basis::Basis, matrix::MacaulayMatrix,
            ht::MonomialHashtable, symbol_ht::MonomialHashtable;
            linalg::Symbol=:sparse)

     convert_hashes_to_columns!(matrix, symbol_ht)

     sort_matrix_rows_decreasing!(matrix) # for pivots,  AB part
     sort_matrix_rows_increasing!(matrix) # for reduced, CD part

     exact_sparse_linear_algebra!(matrix, basis)

     convert_matrix_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
end

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

function find_multiplied_reducer!(
            basis::Basis, matrix::MacaulayMatrix,
            ht::MonomialHashtable, symbol_ht::MonomialHashtable,
            vidx::Int)

    e = symbol_ht.exponents[vidx]
    etmp = ht.exponents[1]
    divmask = symbol_ht.hashdata[vidx].divmask

    blen = basis.ndone     #
    leaddiv = basis.lead  #

    # searching for a poly from state.basis whose leading monom
    # divides the given exponent e
    i = 1
    @label Letsgo
    @inbounds while i <= basis.nlead && (leaddiv[i] & ~divmask) != 0
        i += 1
    end

    # here found polynomial from basis with leading monom
    # dividing symbol_ht.exponents[vidx]
    if i <= basis.nlead
        # reducers index and exponent in hash table
        @inbounds rpoly = basis.gens[basis.nonred[i]]
        @inbounds rexp  = ht.exponents[rpoly[1]]

        @inbounds for j in 1:ht.explen
            # if it actually does not divide and divmask lies
            if e[j] < rexp[j]
                i += 1
                @goto Letsgo
            end
            etmp[j] = e[j] - rexp[j]
        end
        # now etmp = e // rexp in terms of monomias,
        # hash is linear
        h = symbol_ht.hashdata[vidx].hash - ht.hashdata[rpoly[1]].hash

        matrix.uprows[matrix.nup + 1] = multiplied_poly_to_matrix_row!(symbol_ht, ht, h, etmp, rpoly)
        matrix.up2coef[matrix.nup + 1] = basis.nonred[i]

        # upsize matrix
        symbol_ht.hashdata[vidx].idx = 2
        matrix.nup += 1
        i += 1
    end

end

# Given the set of polynomials L and the basis G,
# extends L to contain possible polys for reduction by G,
# and returns it
function symbolic_preprocessing!(
                basis::Basis, matrix::MacaulayMatrix,
                ht::MonomialHashtable, symbol_ht::MonomialHashtable;
                linalg::Symbol=:sparse)

    # check matrix sizes (we want to omit this I guess)

    symbol_load = symbol_ht.load

    nrr    = matrix.ncols
    onrr   = matrix.ncols

    #=
        TODO: matrix_enlarge!
    =#
    while matrix.size <= nrr + symbol_load
        matrix.size *= 2
        resize!(matrix.uprows, matrix.size)
        resize!(matrix.up2coef, matrix.size)
    end

    # for each lcm present in symbolic_ht set on select stage
    i = symbol_ht.offset
    #= First round, we add multiplied polynomials which divide  =#
    #= a monomial exponent from selected spairs                 =#
    while i <= symbol_load
        # not a reducer already
        if symbol_ht.hashdata[i].idx == 0
            symbol_ht.hashdata[i].idx = 1
            matrix.ncols += 1
            find_multiplied_reducer!(basis, matrix, ht, symbol_ht, i)
        end
        i += 1
    end

    #= Second round, we add multiplied polynomials which divide  =#
    #= lcm added on previous for loop                             =#
    while i <= symbol_ht.load
        if matrix.size == matrix.nup
            matrix.size *= 2
            resize!(matrix.uprows, matrix.size)
            resize!(matrix.up2coef, matrix.size)
        end

        symbol_ht.hashdata[i].idx = 1
        matrix.ncols += 1
        find_multiplied_reducer!(basis, matrix, ht, symbol_ht, i)
        i += 1
    end

    # shrink matrix sizes, set constants
    resize!(matrix.uprows, matrix.nup)
    matrix.nrows += matrix.nup - onrr
    matrix.nlow  = matrix.nrows - matrix.nup
    matrix.size  = matrix.nrows

end

function update_pairset!(
            pairset,
            basis,
            ht,
            update_ht,
            idx)

    pl = pairset.load
    bl = idx
    nl = pl + bl
    ps = pairset.pairs

    new_lead  = basis.gens[idx][1]

    # initialize new critical lcms
    plcm = Vector{Int}(undef, bl + 1)

    # for each combination (new_Lead, basis.gens[i][1])
    # generate a pair
    for i in 1:bl-1
        plcm[i] = get_lcm(basis.gens[i][1], new_lead, ht, update_ht)
        deg = update_ht.hashdata[plcm[i]].deg
        newidx = pl + i
        if basis.isred[i] == 0
            ps[newidx] = SPair(i, idx, plcm[i], deg)
        else
            # lcm == 0 will mark redundancy of spair
            ps[newidx] = SPair(i, idx, 0, deg)
        end
    end

    # traverse existing pairs
    for i in 1:pl
        j = ps[i].poly1
        l = ps[i].poly2
        m = max(ps[pl + l].deg, ps[pl + j].deg)

        # if an existing pair is divisible by lead of new poly
        # and has a greater degree than newly generated one
        if is_monom_divisible(ps[i].lcm, new_lead, ht) && ps[i].deg > m
            # mark lcm as 0
            ps[i] = SPair(ps[i].poly1, ps[i].poly2, 0, ps[i].deg)
        end
    end
    # TODO: this can be done faster is we
    # create only non-redundant pairs at first

    # traverse new pairs to check for redundancy
    j = 1
    for i in 1:bl-1
        if basis.isred[i] == 0
            ps[pl + j] = ps[pl + i]
            j += 1
        end
    end

    sort_pairset_by_degree!(pairset, pl + 1, j - 2)

    for i in 1:j - 1
        plcm[i] = ps[pl + i].lcm
    end
    plcm[j] = 0
    pc = j
    pc -= 1

    # mark redundancy of some pairs from plcm
    for j in 1:pc
        # if is not redundant
        if plcm[j] > 0
            check_monomial_division_in_update(plcm, j + 1, pc, plcm[j], update_ht)
        end
    end

    # remove useless pairs from pairset
    # by moving them to the end
    j = 1
    for i in 1:pairset.load
        ps[i].lcm == 0 && continue
        ps[j] = ps[i]
        j += 1
    end

    # assure that basis hashtable can store new lcms
    if ht.size - ht.load <= pc
        enlarge_hash_table!(ht)
    end

    # add new lcms to the basis hashtable,
    # including j and not including pc
    insert_plcms_in_basis_hash_table!(pairset, pl, ht, update_ht, basis, plcm, j, pc+1)

    # mark redundant elements in masis
    nonred = basis.nonred
    lml = basis.nlead
    for i in 1:lml
        if basis.isred[nonred[i]] == 0
            if is_monom_divisible(basis.gens[nonred[i]][1], new_lead, ht)
                basis.isred[nonred[i]] = 1
            end
        end
    end

end

function update_basis!(
            basis,
            ht::MonomialHashtable,
            update_ht::MonomialHashtable)

    # here we could check overall redundancy and update basis.lead

    k = 1
    lead   = basis.lead
    nonred = basis.nonred

    for i in 1:basis.nlead
        if basis.isred[nonred[i]] == 0
            basis.lead[k]   = lead[i]
            basis.nonred[k] = nonred[i]
            k += 1
        end
    end
    basis.nlead = k - 1

    for i in basis.ndone+1:basis.ntotal
        if basis.isred[i] == 0
            lead[k]   = ht.hashdata[basis.gens[i][1]].divmask
            nonred[k] = i
            k += 1
        end
    end

    basis.nlead = k - 1
    basis.ndone = basis.ntotal
end

# checks if element of basis at position idx is redundant
function is_redundant!(
            pairset, basis, ht, update_ht, idx)

    reinitialize_hash_table!(update_ht, 2*idx)

    ps = pairset.pairs

    # lead of new polynomial
    lead_new = basis.gens[idx][1]
    # degree of lead
    lead_deg = ht.hashdata[lead_new].deg

    for i in idx+1:basis.ntotal
        i == idx && continue
        if basis.isred[i] == 1
            continue
        end

        lead_i = basis.gens[i][1]

        if is_monom_divisible(lead_new, lead_i, ht)
            lcm_new = get_lcm(lead_i, lead_new, ht, ht)

            psidx = pairset.load + 1
            ps[psidx] = SPair(i, idx, lcm_new, ht.hashdata[lcm_new].deg)

            basis.isred[idx] = 1
            pairset.load += 1

            return true
        end
    end

    return false
end

function update!(
        pairset::Pairset,
        basis::Basis,
        ht::MonomialHashtable,
        update_ht::MonomialHashtable)

    #=
        Always check redundancy, for now
    =#

    # number of added elements
    npivs = basis.ntotal

    # number of potential critical pairs to add
    npairs = basis.ndone * npivs
    for i in 1:npivs
        npairs += i
    end

    # make sure pairset and update hashtable have enough
    # space to store new pairs
    # TODO: we create too big array, can be fixed
    check_enlarge_pairset!(pairset, npairs)

    # for each new element in basis
    for i in basis.ndone+1:basis.ntotal
        # check redundancy of new poly
        if is_redundant!(pairset, basis, ht, update_ht, i)
            continue
        end
        update_pairset!(pairset, basis, ht, update_ht, i)
    end

    update_basis!(basis, ht, update_ht)
end

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

function filter_redundant!(basis::Basis)
    j = 1
    for i in 1:basis.nlead
        if basis.isred[basis.nonred[i]] == 0
            basis.lead[j] = basis.lead[i]
            basis.nonred[j] = basis.nonred[i]
            #=
            basis.coeffs[j] = basis.coeffs[basis.nonred[i]]
            basis.gens[j] = basis.gens[basis.nonred[i]]
            =#
            j += 1
        end
    end
    basis.nlead = j - 1
    @assert basis.ndone == basis.ntotal
    # basis.ndone = basis.ntotal = basis.nlead
    basis
end

function export_basis_data(basis::Basis, ht::MonomialHashtable)
    exps   = Vector{Vector{Vector{UInt16}}}(undef, basis.nlead)
    coeffs = Vector{Vector{UInt64}}(undef, basis.nlead)

    for i in 1:basis.nlead
        idx = basis.nonred[i]
        poly = basis.gens[idx]
        exps[i] = Vector{Vector{UInt16}}(undef, length(poly))
        for j in 1:length(poly)
            exps[i][j] = ht.exponents[poly[j]]
        end
        coeffs[i] = basis.coeffs[idx]
    end

    exps, coeffs
end

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

# given input exponent and coefficient vectors hashes exponents into `ht`
# and then constructs hashed polynomial vectors for `basis`
function fill_data!(basis, ht, exponents, coeffs)
    ngens = length(exponents)

    for i in 1:ngens
        while length(exponents[i]) >= ht.size - ht.load
            enlarge_hash_table!(ht)
        end
        # ht.exponents one can be reallocated together with ht,
        # so we need to reset it on each iteration
        etmp = ht.exponents[1]

        nterms = length(coeffs[i])
        basis.coeffs[i] = coeffs[i]
        basis.gens[i] = Vector{Int}(undef, nterms)
        poly = basis.gens[i]
        for j in 1:nterms
            poly[j] = insert_in_hash_table!(ht, exponents[i][j])
        end

        # sort terms (not needed),
        # beautify coefficients (not needed),
        # and all
    end

    # the initial update traverses all elements in
    # basis.ndone+1:basis.ntotal
    # which is 1:ngens
    basis.ntotal = ngens
end

# Initializes Basis and MonomialHashtable structures,
# fills input data from exponents and coeffs
#
# Hashtable initial size is set to tablesize
function initialize_structures(
            ring::PolyRing,
            exponents::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            rng::Random.AbstractRNG,
            tablesize::Int)

    # basis for storing basis elements,
    # pairset for storing critical pairs of basis elements to assess,
    # hashtable for hashing monomials occuring in the basis
    basis    = initialize_basis(ring, length(exponents))
    basis_ht = initialize_basis_hash_table(ring, rng, initial_size=tablesize)

    # filling the basis and hashtable with the given inputs
    fill_data!(basis, basis_ht, exponents, coeffs)

    # every monomial in hashtable is associated with its divmask
    # to perform divisions faster. Filling those
    fill_divmask!(basis_ht)

    # sort input, smaller leading terms first
    sort_gens_by_lead_increasing!(basis, basis_ht)

    # divide each polynomial by leading coefficient
    normalize_basis!(basis)

    basis, basis_ht
end

# Initializes Basis with the given hashtable,
# fills input data from exponents and coeffs
function initialize_structures(
            ring::PolyRing,
            exponents::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            rng::Random.AbstractRNG,
            tablesize::Int,
            present_ht::MonomialHashtable)

    # basis for storing basis elements,
    # pairset for storing critical pairs of basis elements to assess,
    # hashtable for hashing monomials occuring in the basis
    basis    = initialize_basis(ring, length(exponents))

    # filling the basis and hashtable with the given inputs
    fill_data!(basis, present_ht, exponents, coeffs)

    # sort input, smaller leading terms first
    sort_gens_by_lead_increasing!(basis, present_ht)

    # divide each polynomial by leading coefficient
    # We do not need normalization for normal forms
    # normalize_basis!(basis)

    basis, present_ht
end

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

function reducegb_f4!(
            basis::Basis, matrix::MacaulayMatrix,
            ht::MonomialHashtable, symbol_ht::MonomialHashtable)

    etmp  = ht.exponents[1]
    etmp .= UInt16(0)
    # etmp is now set to zero, and has a zero hash

    reinitialize_matrix!(matrix, basis.nlead)
    uprows = matrix.uprows

    #=
    @warn "entering reduce"
    dump(basis, maxdepth=5)
    @warn "ht"
    println(ht.exponents[1:10])
    @warn "symbol"
    println(symbol_ht.exponents[1:10])
    println("NONRED ", basis.nonred)
    =#

    # add all non redundant elements from basis
    # as matrix upper rows
    for i in 1:basis.nlead
        matrix.nrows += 1
        uprows[matrix.nrows] = multiplied_poly_to_matrix_row!(
                                    symbol_ht, ht, UInt32(0), etmp,
                                    basis.gens[basis.nonred[i]])

        matrix.up2coef[matrix.nrows] = basis.nonred[i]
        # set lead index as 1
        symbol_ht.hashdata[uprows[matrix.nrows][1]].idx = 1
    end

    # needed for correct counting in symbol
    matrix.ncols = matrix.nrows
    matrix.nup = matrix.nrows

    #=
    @warn "after multiplied_poly_to_matrix_row"
    dump(basis, maxdepth=5)
    @warn "matrix"
    dump(matrix, maxdepth=5)
    @warn "ht"
    println(ht.exponents[1:10])
    =#

    symbolic_preprocessing!(basis, matrix, ht, symbol_ht)
    # all pivots are unknown
    for i in symbol_ht.offset:symbol_ht.load
        symbol_ht.hashdata[i].idx = 1
    end

    #=
    @warn "after symbolic"
    dump(basis, maxdepth=5)
    @warn "matrix"
    dump(matrix, maxdepth=5)
    @warn "ht"
    println(ht.exponents[1:10])
    =#

    # x1*x2^6 + 413612941*x1*x2^3
    #  x1^2*x2^3 + 174101409*x1*x2^5
    convert_hashes_to_columns!(matrix, symbol_ht)
    matrix.ncols = matrix.nleft + matrix.nright

    #=
    @warn "after convert"
    dump(matrix, maxdepth=5)
    =#

    sort_matrix_rows_decreasing!(matrix)

    interreduce_matrix_rows!(matrix, basis)

    #=
    @warn "after interreduce"
    dump(matrix, maxdepth=5)
    =#

    convert_matrix_rows_to_basis_elements_use_symbol!(matrix, basis)

    # no longer need in two hashtables
    ht = symbol_ht

    basis.ntotal = matrix.npivots + basis.ndone
    basis.ndone = matrix.npivots

    #=
    @warn "basis"
    dump(basis, maxdepth=5)
    =#

    #= we may have added some multiples of reduced basis polynomials
    * from the matrix, so we get rid of them. =#
    k = 0
    i = 1
    @label Letsgo
    while i <= basis.ndone
        for j in 1:k
            if is_monom_divisible(
                    basis.gens[basis.ntotal - i + 1][1],
                    basis.gens[basis.nonred[j]][1],
                    ht)

                i += 1
                @goto Letsgo
            end
        end
        k += 1
        basis.nonred[k] = basis.ntotal - i + 1
        # xd
        basis.lead[k] = ht.hashdata[basis.gens[basis.ntotal - i + 1][1]].divmask
        i += 1
    end
    basis.nlead = k

    # TODO
    sort_gens_by_lead_increasing_in_reduce!(basis, ht)
end

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

"""
Main function to calculate the Groebner basis of the given polynomial ideal
in degrevlex ordering and finite fields coefficients.

Parameters
.

"""
function f4(ring::PolyRing,
            exponents::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            rng::Rng,
            reduced::Bool;
            tablesize::Int=2^16
            ) where {Rng<:Random.AbstractRNG}

    # initialize basis and hashtable with input polynomials,
    # fields are not meaningful yet and will be set during update! step
    basis, ht = initialize_structures(
                        ring, exponents, coeffs, rng, tablesize)

    # matrix storing coefficients in rows
    # wrt columns representing the current monomial basis
    matrix = initialize_matrix(ring)

    # initialize hash tables for update and symbolic preprocessing steps
    update_ht  = initialize_secondary_hash_table(ht)
    symbol_ht  = initialize_secondary_hash_table(ht)

    # a set to store critical pairs of polynomials to be reduced
    pairset = initialize_pairset()

    # makes basis fields valid,
    # does not copy,
    # checks for redundancy of new elems
    update!(pairset, basis, ht, update_ht)

    #=
    @warn "initial update"
    dump(basis, maxdepth = 5)
    =#

    d = 0
    # while there are pairs to be reduced
    while !isempty(pairset)
        d += 1
        @debug "F4 ITER $d"
        @debug "Available $(pairset.load) pairs"

        #=
        @warn "start of $d"
        dump(basis, maxdepth = 5)
        println(ht.exponents[1:10])
        =#

        # selects pairs for reduction from pairset following normal strategy
        # (minimal lcm degrees are selected),
        # and puts these into the matrix rows
        select_normal!(pairset, basis, matrix, ht, symbol_ht)
        @debug "Selected $(divexact(matrix.nrows, 2)) pairs"

        #=
        @warn "noxd $d"
        dump(basis, maxdepth = 5)
        dump(matrix, maxdepth = 5)
        println(ht.exponents[1:10])
        println(symbol_ht.exponents[1:10])
        =#

        symbolic_preprocessing!(basis, matrix, ht, symbol_ht)
        @debug "Matrix of size $((matrix.nrows, matrix.ncols)), density TODO"

        #=
        @warn "xd $d"
        dump(basis, maxdepth = 5)
        dump(matrix, maxdepth = 5)
        println(ht.exponents[1:10])
        println(symbol_ht.exponents[1:10])
        =#

        # reduces polys and obtains new potential basis elements
        reduction!(basis, matrix, ht, symbol_ht)
        @debug "Matrix reduced, density TODO"

        #@warn "after reduction"
        #dump(matrix, maxdepth = 5)

        #=
        printstyled("nice basis\n", color=:red)
        for (_, cfs, poly) in zip(1:basis.ntotal, basis.coeffs, basis.gens)
            for i in 1:length(poly)
                c = cfs[i]
                e = ht.exponents[poly[i]]
                print("$c*$e + ")
            end
            println("")
        end
        println(basis.nonred)
        =#

        # update the current basis with polynomials produced from reduction,
        # does not copy,
        # checks for redundancy
        update!(pairset, basis, ht, update_ht)

        #=
        @warn "after update"
        printstyled("nice basis\n", color=:red)
        for (_, cfs, poly) in zip(1:basis.ntotal, basis.coeffs, basis.gens)
            for i in 1:length(poly)
                c = cfs[i]
                e = ht.exponents[poly[i]]
                print("$c*$e + ")
            end
            println("")
        end
        println(basis.nonred)
        =#

        # TODO: is this okay hm ?
        # to be changed
        matrix    = initialize_matrix(ring)
        symbol_ht = initialize_secondary_hash_table(ht)
        # clear symbolic hashtable
        # clear matrix

        if d > 1000
            @error "Something is probably wrong in f4.."
            break
        end
    end

    # remove redundant elements
    filter_redundant!(basis)

    #=
    @warn "after filter red"
    printstyled("nice basis\n", color=:red)
    for (_, cfs, poly) in zip(1:basis.ntotal, basis.coeffs, basis.gens)
        for i in 1:length(poly)
            c = cfs[i]
            e = ht.exponents[poly[i]]
            print("$c*$e + ")
        end
        println("")
    end
    println("NONRED ", basis.nonred)
    println(basis.isred)
    println(basis.nlead)
    =#

    if reduced
        reducegb_f4!(basis, matrix, ht, symbol_ht)
        ht = symbol_ht
    end

    #=
    @warn "after reduce"
    printstyled("nice basis\n", color=:red)
    for (_, cfs, poly) in zip(1:basis.ntotal, basis.coeffs, basis.gens)
        for i in 1:length(poly)
            c = cfs[i]
            e = ht.exponents[poly[i]]
            print("$c*$e + ")
        end
        println("")
    end
    println(basis.nonred)
    =#

    export_basis_data(basis, ht)
end

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

function select_isgroebner!(
        pairset::Pairset, basis::Basis, matrix::MacaulayMatrix,
        ht::MonomialHashtable, symbol_ht::MonomialHashtable)

    ps = pairset.pairs
    npairs = pairset.load

    sort_pairset_by_lcm!(pairset, npairs, ht)

    reinitialize_matrix!(matrix, npairs)

    uprows  = matrix.uprows
    lowrows = matrix.lowrows

    # polynomials from pairs in order (p11, p12)(p21, p21)
    # (future rows of the matrix)
    gens   = Vector{Int}(undef, 2*npairs)

    # buffer !
    etmp = ht.exponents[1]
    i = 1
    while i <= npairs
        matrix.ncols += 1
        load = 1
        lcm = ps[i].lcm
        j = i

        # we collect all generators with same lcm into gens
        while j <= npairs && ps[j].lcm == lcm
            gens[load] = ps[j].poly1
            load += 1
            gens[load] = ps[j].poly2
            load += 1
            j += 1
        end
        load -= 1

        # sort by number in the basis (by=identity)
        sort_generators_by_position!(gens, load)

        # now we collect reducers, and reduced

        # first generator index in groebner basis
        prev = gens[1]
        # first generator in hash table
        poly = basis.gens[prev]
        # first generator lead monomial index in hash data
        vidx = poly[1]

        # first generator exponent
        eidx = ht.exponents[vidx]
        # exponent of lcm corresponding to first generator
        elcm = ht.exponents[lcm]
        for u in 1:ht.explen
            etmp[u] = elcm[u] - eidx[u]
        end
        # now etmp contents complement to eidx in elcm

        # hash of complement
        htmp = ht.hashdata[lcm].hash - ht.hashdata[vidx].hash

        # add row as a reducer
        matrix.nup += 1
        uprows[matrix.nup] = multiplied_poly_to_matrix_row!(symbol_ht, ht, htmp, etmp, poly)
        # map upper row to index in basis
        matrix.up2coef[matrix.nup] = prev

        # mark lcm column as reducer in symbolic hashtable
        symbol_ht.hashdata[uprows[matrix.nup][1]].idx = 2
        # increase number of rows set
        matrix.nrows += 1

        # over all polys with same lcm,
        # add them to the lower part of matrix
        for k in 1:load
            # duplicate generator,
            # we can do so as long as generators are sorted
            if gens[k] == prev
                continue
            end

            # if the table was reallocated
            elcm = ht.exponents[lcm]

            # index in gb
            prev = gens[k]
            # poly of indices of monoms in hash table
            poly = basis.gens[prev]
            vidx = poly[1]
            # leading monom idx
            eidx = ht.exponents[vidx]
            for u in 1:ht.explen
                etmp[u] = elcm[u] - eidx[u]
            end

            htmp = ht.hashdata[lcm].hash - ht.hashdata[vidx].hash

            # add row to be reduced
            matrix.nlow += 1
            lowrows[matrix.nlow] = multiplied_poly_to_matrix_row!(symbol_ht, ht, htmp, etmp, poly)
            # map lower row to index in basis
            matrix.low2coef[matrix.nlow] = prev

            symbol_ht.hashdata[lowrows[matrix.nlow][1]].idx = 2

            matrix.nrows += 1
        end

        i = j
    end

    resize!(matrix.lowrows, matrix.nrows - matrix.ncols)
end

function isgroebner_f4!(ring::PolyRing,
                        basis::Basis,
                        ht::MonomialHashtable,
                        rng)

    matrix = initialize_matrix(ring)
    symbol_ht = initialize_secondary_hash_table(ht)
    update_ht  = initialize_secondary_hash_table(ht)

    pairset = initialize_pairset()

    update!(pairset, basis, ht, update_ht)

    if pairset.load == 0
        return true
    end

    select_isgroebner!(pairset, basis, matrix, ht, symbol_ht)

    symbolic_preprocessing!(basis, matrix, ht, symbol_ht)

    convert_hashes_to_columns!(matrix, symbol_ht)

    sort_matrix_rows_increasing!(matrix)
    sort_matrix_rows_decreasing!(matrix) # for pivots,  AB part

    return exact_sparse_linear_algebra_isgroebner!(matrix, basis)
end

function isgroebner_f4(
            ring::PolyRing,
            exps::Vector{Vector{Vector{UInt16}}},
            coeffs::Vector{Vector{UInt64}},
            rng)

    basis, ht = initialize_structures(ring, exps,
                                        coeffs, rng, 2^16)

    isgroebner_f4!(ring, basis, ht, rng)
end