text:
module Jarrett
  import Base.+

  struct JarrettNumber <: Number
    y::Number
  end

  Base.:+(x::Number, j::JarrettNumber) = JarrettNumber(x + j.y)

  export JarrettNumber
end

module Peter
  import Base.+

  using Hyperspecialize

  struct PeterNumber <: Number
    x::Number
  end

  @replicable Base.:+(p::PeterNumber, y::@hyperspecialize(Number)) = PeterNumber(p.x + y)

  export PeterNumber
end

using .Peter
using .Jarrett

p = PeterNumber(1.0) + 3
j = 2.0 + JarrettNumber(2.0)
friends = p + j

println("Peter and Jarrett are friends! $friends")


text:
# This file includes the plugin for standard deviation of data

using Statistics

"""
    $TYPEDEF

Constructs a `Std` plugin. The [`process(plg::Std, x)`](@ref) function takes the standard deviation of the input data `x`
along dimension `dims`.

# Fields

    $TYPEDFIELDS
"""
struct Std <: AbstractPlugin
    "Dimension"
    dims::Int
end
Std(;dims::Int=1) = Std(dims)

show(io::IO, plg::Std) = print(io, "Mean(dims:$(plg.dims))")

"""
    $SIGNATURES

Returns the standard deviation of `x` along the dimension `plg.dims`.

# Example
```julia
julia> x = collect(reshape(1:16, 4,4))
4×4 Array{Int64,2}:
 1  5   9  13
 2  6  10  14
 3  7  11  15
 4  8  12  16

julia> plg = Plugins.Std(dims=1)
Mean(dims:1)

julia> process(plg, x)
1×4 Array{Float64,2}:
 1.29099  1.29099  1.29099  1.29099
```
"""
process(plg::Std, x) = std(x, dims=plg.dims)

ext: jl


text:
# Utilities that should probably go into Core.Compiler
using Core.Compiler: CFG, BasicBlock, BBIdxIter

function Base.push!(cfg::CFG, bb::BasicBlock)
    @assert cfg.blocks[end].stmts.stop+1 == bb.stmts.start
    push!(cfg.blocks, bb)
    push!(cfg.index, bb.stmts.start)
end

function Core.Compiler.NewInstruction(node)
    Core.Compiler.NewInstruction(node, Any)
end

function Base.setproperty!(x::Core.Compiler.Instruction, f::Symbol, v)
    Core.Compiler.setindex!(x, v, f)
end

function Base.getproperty(x::Core.Compiler.Instruction, f::Symbol)
    Core.Compiler.getindex(x, f)
end

function Base.setindex!(ir::Core.Compiler.IRCode, ni::NewInstruction, i::Int)
    stmt = ir.stmts[i]
    stmt.inst = ni.stmt
    stmt.type = ni.type
    stmt.flag = ni.flag
    stmt.line = something(ni.line, 0)
    ni
end

function Base.push!(ir::IRCode, ni::NewInstruction)
    # TODO: This should be a check in insert_node!
    @assert length(ir.new_nodes.stmts) == 0
    ir[Core.Compiler.add!(ir.stmts)] = ni
    ir
end

function Base.iterate(it::Iterators.Reverse{BBIdxIter},
        (bb, idx)::Tuple{Int, Int}=(length(it.itr.ir.cfg.blocks), length(it.itr.ir.stmts)+1))
    idx == 1 && return nothing
    active_bb = it.itr.ir.cfg.blocks[bb]
    if idx == first(active_bb.stmts)
        bb -= 1
    end
    return (bb, idx - 1), (bb, idx - 1)
end


text:
# L-infinity norm

export NormLinf

"""
**``L_∞`` norm**

    NormLinf(λ=1.0)

Returns the function
```math
f(x) = λ⋅\\max\\{|x_1|, …, |x_n|\\},
```
for a nonnegative parameter `λ`.
"""
NormLinf(lambda::R=1.0) where {R <: Real} = Conjugate(IndBallL1(lambda))

function (f::Conjugate{IndBallL1{R}})(x::AbstractArray{S}) where {R <: Real, S <: RealOrComplex}
  return (f.f.r)*norm(x, Inf)
end

function gradient!(y::AbstractArray{T}, f::Conjugate{IndBallL1{R}}, x::AbstractArray{T}) where {T <: RealOrComplex, R <: Real}
  absxi, i = findmax(abs.(x)) # Largest absolute value
  y .= 0
  y[i] = f.f.r*sign(x[i])
  return f.f.r*absxi
end

fun_name(f::Postcompose{Conjugate{IndBallL1{R}}, R}) where {R <: Real} = "weighted L-infinity norm"
fun_expr(f::Postcompose{Conjugate{IndBallL1{R}}, R}) where {R <: Real} = "x ↦ λ||x||_∞ = λ⋅max(abs(x))"
fun_params(f::Postcompose{Conjugate{IndBallL1{R}}, R}) where {R <: Real} = "λ = $(f.a*(f.f.f.r))"


text:
type DistinctSorted
    it
    prev
end


text:
# This file is not an active part of the package. This is the code
# that uses the Clang.jl package to wrap Silo using the headers.

# Find all headers
# To run the script from Julia console:
# include(joinpath(Pkg.dir("Silo"), "src", "wrap_Silo.jl"));
push!(Libdl.DL_LOAD_PATH, "/usr/local/Cellar/llvm/4.0.0_1/lib")
using Clang.wrap_c
const wdir = dirname(@__FILE__)
const pkgbasedir = joinpath(wdir, "..")

# `outpath` specifies, where the julian wrappers would be generated.
# If the generated .jl files are ok, they have to be copied to the "src" folder
# overwriting the old ones
const outpath = joinpath(pkgbasedir, "new")
rm(outpath, recursive = true, force = true)
mkpath(outpath)
mkpath(outpath)

const incpath = joinpath(pkgbasedir, "deps", "usr", "include")
if !isdir(incpath)
  error("Run Pkg.build(\"Silo\") before trying to wrap C headers.")
end

info("Scanning Silo headers in $incpath...")
const Silo_header_files = ["silo", "pmpio"]
const Silo_headers =[joinpath(incpath, y) for y in (string(x, ".h") for x in Silo_header_files)]
# map(x->joinpath(incpath, x), readdir(incpath))

const clang_path = "/usr/local/Cellar/llvm/4.0.0_1/lib/clang/4.0.0/" # change to your clang location
const includes = [
    joinpath(clang_path, "include"),
    incpath,
    "/usr/local/Cellar/open-mpi/2.1.0/include/"
]

function find_outfile(s)
    joinpath(outpath, string(first(splitext(basename(s))), ".jl"))
end
find_libfile(s) = "libsilo"

const context = wrap_c.init(
    headers = Silo_headers,
    common_file = joinpath(outpath, "libsilo_types.jl"),
    clang_args = [
        "-D", "__STDC_LIMIT_MACROS",
        "-D", "__STDC_CONSTANT_MACROS",
        # "-v"
    ],
    # clang_diagnostics = true,
    header_library = find_libfile,
    header_outputfile = find_outfile,
    clang_includes = includes,
)

info("Generating .jl wrappers for Silo in $outpath...")
run(context)
info("Done generating .jl wrappers for Silo in $outpath")


text:
function read_grd(filename)
    a = readlines(filename)
    filter!(x->!isempty(x) && !(x[1] == '['), a)
    Dict((c = match(r"([^=]+)=(.*)", st); c.captures[1] => strip(c.captures[2])) for st in a)
end

function write_grd(grdfile::String, vsr::VerySimpleRaster)
    grdfile[end-3:end] ∈ (".grd", ".gri") && (grdfile = grdfile[1:end-4])
    open(grdfile*".gri", "w") do IO
        write(IO, vsr.mat)
    end
    write_grd_header(grdfile, vsr)
end

function write_grd_header(grdfile, vsr::VerySimpleRaster)
    bb = bbox(vsr)
    write_grd_header(grdfile,
        nrows= size(vsr.mat, 2),
        ncols= size(vsr.mat, 1),
        xmin= bb.xmin,
        ymin= bb.ymin,
        xmax= bb.xmax,
        ymax= bb.ymax,
        projection= vsr.projection,
        datatype= eltype(vsr.mat),
        nodata= vsr.nodata,
        minvalue = minimum(filter(x->x!=vsr.nodata, vsr.mat)),
        maxvalue= maximum(filter(x->x!=vsr.nodata, vsr.mat))
    )
end

function write_grd_header(grdfile;
                        nrows = error("unspecified argument nrows"), ncols = error("unspecified argument ncols"),
                        xmin = error("unspecified argument xmin"), ymin = error("unspecified argument ymin"),
                        xmax = error("unspecified argument xmax"), ymax = error("unspecified argument ymax"),
                        projection = error("unspecified argument projection"), datatype = error("unspecified argument datatype"),
                        nodata = error("unspecified argument nodata"),
                        minvalue = error("unspecified argument minvalue"), maxvalue = error("unspecified argument maxvalue"))
    grdfile[end-3:end] == ".grd" || (grdfile = grdfile * ".grd")
    a = open(grdfile, "w") do IO
        write(IO, """
        [general]
        creator=Julia Package VerySimpleRasters
        created= $(string(now()))
        [georeference]
        nrows= $(nrows)
        ncols= $(ncols)
        xmin= $(xmin)
        ymin= $(ymin)
        xmax= $(xmax)
        ymax= $(ymax)
        projection= $(projection)
        [data]
        datatype= $(rev_datatype_translation[datatype])
        nodatavalue= $(nodata)
        byteorder= little
        nbands= 1
        minvalue= $(minvalue)
        maxvalue= $(maxvalue)
        [description]
        layername= layer1""")
    end
end

function VerySimpleRaster(grdfile::String)
    @assert grdfile[end-3:end] ∈ (".grd", ".gri")
    fname = grdfile[1:end-4]
    data = read_grd(fname*".grd")
    ncols = parse(Int, data["nrows"])
    nrows = parse(Int, data["ncols"])
    dt = datatype_translation[data["datatype"]]
    fname = abspath(fname*".gri")
    mat = open(fname, "r") do IO
        Mmap.mmap(IO, Matrix{dt}, (nrows, ncols))
    end
    xmin, xmax = parse(Float64, data["xmin"]), parse(Float64, data["xmax"])
    ymin, ymax = parse(Float64, data["ymin"]), parse(Float64, data["ymax"])
    celly = (ymax - ymin) / ncols
    cellx = (xmax - xmin) / nrows
    VerySimpleRaster(mat, parse(dt, data["nodatavalue"]), xmin, ymin, cellx, celly, data["projection"], fname)
end

const datatype_translation = Dict{String, DataType}(
    "LOG1S" => Bool,
    "INT1S" => Int8,
    "INT2S" => Int16,
    "INT4S" => Int32,
    "INT8S" => Int64,
    "INT1U" => UInt8,
    "INT2U" => UInt16,
    "FLT4S" => Float32,
    "FLT8S" => Float64
)

const rev_datatype_translation = Dict{DataType, String}(v => k for (k,v) in datatype_translation)

"""
    writeraster(file::String, vsr::VerySimpleRaster)

Writes the raster as an R .grd file. Under the hood simply copies the temporary
on-disk version of the raster to `file`.
"""
function writeraster(fname::String, vsr::VerySimpleRaster; force = false)
    fname[end-4:end] ∈ (".gri", ".grd") && (fname = fname[1:end-4])
    Mmap.sync!(vsr.mat)
    cp(vsr.filename, fname*".gri", force = force)
    cp(vsr.filename[1:end-4]*".grd", fname*".grd", force = force)
end


text:
#===============
    PROBLEMS
===============#
mutable struct Column{T}
    ind::Vector{Int}
    nzval::Vector{T}

    Column{T}() where T = new{T}(Int[], T[])
end

abstract type AbstractProblem{T <: Real} end
abstract type AbstractStandardProblem{T <: Real} <: AbstractProblem{T} end

"""
    Problem{T <: Real}

Data structure to represent the linear program:
    minmize     c' * x
    subject to  A * x = b and lo ≦ x ≦ hi
Suppose the constraint matrix A is m x n.
    `nc` is the number of constraints, i.e. nc = m.
    `nv` is the number of variables,   i.e. nv = n.
"""
mutable struct Problem{T <: Real} <: AbstractProblem{T}
    nc::Int
    nv::Int

    cols::Vector{Column{T}} # columns of A
    b::Vector{T}
    c::Vector{T}

    lo::Vector{T}
    hi::Vector{T}

    Problem{T}() where T <: Real = new{T}(
        0, 0,
        Column{T}[], T[], T[],
        T[], T[]
    )
end

function Problem{T}(
    nc::Int, nv::Int,
    cols::Vector{Column{T}}, b::Vector{T},
    c::Vector{T},
    lo::Vector{T}, hi::Vector{T}
) where T <: Real
    lp = Problem{T}()
    lp.nc = nc
    lp.nv = nv
    lp.cols = cols
    lp.c = c
    lp.lo = lo
    lp.hi = hi

    lp
end

"""
    StandardProblem{T <: Real}

Data structure to respresent the standard form of a linear program:
    minimize    c' * x
    subject to  A * x = b and x ≧ 0
"""
mutable struct StandardProblem{T <: Real} <: AbstractStandardProblem{T}
    nc::Int # number of constraints
    nv::Int # number of variables
    nu::Int # number of upper-bounded variables

    upb_ind::Vector{Int} # indices of upper bounded variables
    upb_val::Vector{T}   # value of upper bound

    A::SparseMatrixCSC{T}
    b::Vector{T}
    c::Vector{T}
end

#================
    RESIDUALS
================#
mutable struct Residuals{T <: Real}
    rp::Vector{T} # τ * b - A * x
    ru::Vector{T} # τ * u - v - U * x
    rd::Vector{T} # τ * c - A' * y - s + U * x
    rg::T # c' * x - b' * λ - u * w + κ

    # norms of above residuals
    rpn::T
    run::T
    rdn::T
    rgn::T

    Residuals{T}(nc::Int, nv::Int, nu::Int) where T <: Real = new{T}(
        zeros(T, nc), zeros(T, nu), zeros(T, nv), zero(T),
        zero(T), zero(T), zero(T), zero(T),
    )
end

#==============
    ITERATE
==============#
mutable struct Iterate{T <: Real}
    nc::Int # number of constraints
    nv::Int # number of variables
    nu::Int # number of upper-bounded variables

    # primal
    x::Vector{T} # given variables
    v::Vector{T} # upper bound slacks

    # dual
    λ::Vector{T} # dual variables
    s::Vector{T} # dual slacks
    w::Vector{T} # dual upperbound slacks

    τ::T # primal homogenous constant
    κ::T # dual   homogenous constant

    μ::T # duality measure

    # allocate memory for all vectors
    Iterate{T}(nc::Int, nv::Int, nu::Int) where T <: Real = new{T}(
        nc, nv, nu,
        zeros(T, nv), zeros(T, nu),
        zeros(T, nc), zeros(T, nv), zeros(T, nu),
        zero(T), zero(T), zero(T)
    )
end

#=================
    TOLERANCES
=================#
mutable struct Tolerances{T <: Real}
    εp::T
    εd::T
    εg::T
    εi::T

    Tolerances{T}(tol::T = sqrt(eps(T))) where T <: Real = new{T}(tol, tol, tol, tol)
end

#=============
    STATUS
=============#
@enum SolverStatus begin
    SolverUndefined

    # problem status
    SolverOptimal
    SolverPrimalInfeasible
    SolverDualInfeasible

    # computation status
    SolverExceededIterations
    SolverNumericalInstability
    SolverExceededMemory
end

@enum IterateStatus begin
    IterateUndefined

    IterateOptimal
    IterateFeasible
    IterateInfeasible
end

#================
    SOLUTION
================#
mutable struct Solution{T <: Real}
    x::Vector{T}
    status::Bool

    # standard form
    A_std::Union{Nothing, AbstractMatrix{T}}
    b_std::Vector{T}
    c_std::Vector{T}
    x_std::Vector{T}
    λ_std::Vector{T}
    s_std::Vector{T}

    Solution{T}() where T = new{T}(
        T[], false,
        nothing, T[], T[], T[], T[], T[]
    )
end

#=============
    SOLVER
=============#
mutable struct Solver{T <: Real}
    lp::StandardProblem{T}

    iter::Iterate{T}
    res::Residuals{T}
    tols::Tolerances{T}

    niter::Int # number of iterations

    # status
    status::SolverStatus
    status_primal::IterateStatus
    status_dual::IterateStatus

    # linear solver -- adapted from Tulip.jl
    ls::AbstractLinearSolver{T}
    regP::Vector{T} # primal regularization
    regD::Vector{T} # dual regularization
    regG::T         # gap regularization

    function Solver{T}(
        lp::StandardProblem,
        iter::Iterate{T},
        tols::Tolerances{T}
    ) where T <: Real
        solv = new{T}()

        solv.lp = lp
        solv.iter = iter
        solv.res = Residuals{T}(lp.nc, lp.nv, lp.nu)
        solv.tols = tols

        solv.niter = 0

        solv.status = SolverUndefined
        solv.status_primal = IterateUndefined
        solv.status_dual = IterateUndefined

        # defaults to CHOLMOD backend for Float64 type
        solv.ls = AbstractLinearSolver(
            DefaultBackend(),
            DefaultSystem(),
            lp.A
        )

        # starting regularizations
        solv.regP = ones(T, lp.nv)
        solv.regD = ones(T, lp.nc)
        solv.regG = oneunit(T)

        solv
    end
end


text:
"""
    MutableRange(r::AbstractRange)

Wraps a parent range making it mutable.
"""
mutable struct MutableRange{T,P<:AbstractRange{T}} <:AbstractRange{T}
    parent::P
end

mrange(start; kwargs...) = MutableRange(range(start; kwargs...))
function mrange(start, stop; kwargs...)
    if isempty(kwargs)
        return MutableRange(start:stop)
    else
        return MutableRange(range(start, stop; kwargs...))
    end
end

const UnitMRange{T} = MutableRange{T,UnitRange{T}}
const StepMRange{T,S} = MutableRange{T,StepRange{T,S}}
const LinMRange{T} = MutableRange{T,LinRange{T}}
const StepMRangeLen{T,R,S} = MutableRange{T,StepRangeLen{T,R,S}}

const MRange = Union{DynamicAxis,MutableRange}

ArrayInterface.can_change_size(::Type{T}) where {T<:MutableRange} = true

Base.reverse!(r::MutableRange) = setfield!(r, :parent, reverse(parent(r)))
Base.sort!(r::MutableRange) = setfield!(r, :parent, sort(parent(r)))

Base.isempty(r::MutableRange) = isempty(parent(r))

Base.empty!(r::MutableRange) = setfield!(r, :parent, empty(r))

_empty(r::LinRange) = LinRange(first(r), last(r), 0)
_empty(r::StepRangeLen) = StepRangeLen(r.ref, r.step, 0, r.offset)
_empty(r::StepRange)  = StepRange(r.start, r.step, r.start - step(r))
_empty(r::UnitRange{T}) where {T} = UnitRange(first(r), first(r) - one(T))
_empty(r::OneTo{T}) where {T} = OneTo(zero(T))

Base.parent(x::MutableRange) = getfield(x, :parent)

ArrayInterface.parent_type(::Type{MutableRange{T,R}}) where {T,R} = R

@inline Base.getproperty(x::MutableRange, s::Symbol) = getproperty(parent(x), s)

function Base.promote_rule(::Type{T1}, ::Type{T2}) where {T1<:MutableRange,T2<:MutableRange}
    return promote_rule(parent_type(T1), parent_type(T2))
end


text:
@testset "returnlevel.jl" begin
    @testset "Base.show(io, obj)" begin
        # Print ReturnLevel does not throw
        fm = MaximumLikelihoodEVA(BlockMaxima(Variable("y", [1])), [1.0, 1.0, 0.1])
        rl = ReturnLevel(Extremes.BlockMaximaModel(fm), 10, [1.0])

        buffer = IOBuffer()
        @test_logs Base.show(buffer, rl)

    end

end


text:
"""
    spinoperators(s)

For spin `s` returns the operators (`sz`,`sp`, `sm`). `s` has to be
integer or half-integer
"""
function spinoperators(s;
                       symmetry::Type{<:AbstractCharge}=Trivial)
    d = Int(round(2*s + 1))
    ms = [s-n+1 for n in 1:d]
    amps = [sqrt((s+m)*(s-m+1)) for m in ms[1:d-1]]

    if symmetry == Trivial
        sz = diagm(0 => ms)
        sp = diagm(1 => amps)
        sm = diagm(-1 => amps)
        return sz, sp, sm
    elseif symmetry == U1
        V = U1Space(c=>1 for c in 0:d-1)
        sz_blocks = SortedDict{Sector{U1, 2}, Matrix{Float64}}()
        sp_blocks = SortedDict{Sector{U1, 2}, Matrix{Float64}}()
        sm_blocks = SortedDict{Sector{U1, 2}, Matrix{Float64}}()
        for c in 0:d-1
            sz_blocks[Sector{U1}(c, c)] = reverse(ms)[c+1] * ones(1,1)
        end
        for c in 0:d-2
            sp_blocks[Sector{U1}(c+1, c)] = amps[c+1] * ones(1,1)
            sm_blocks[Sector{U1}(c, c+1)] = amps[c+1] * ones(1,1)
        end
        sz = SymMatrix(zero(U1), (V,dual(V)), sz_blocks)
        sp = SymMatrix(U1(+1), (V,dual(V)), sp_blocks)
        sm = SymMatrix(U1(-1), (V,dual(V)), sm_blocks)
        return sz, sp, sm
    else
        throw(ArgumentError("Unkown symmetry : $symmetry"))
    end
end

######################
### spin-1/2 stuff ###

"""
    permutespins(op, perm)
Given an operator `op` on spins 1...n return the same operator on the
given permutation `perm` of those spins.

"""
function permutespins(op::AbstractMatrix, perm::Vector{Int})
    size(op, 1) == size(op, 2) || error("square matrix please!")
    n = Int(log2(size(op,1)))
    sort(perm) == collect(1:n) || error("not a valid permutation!")

    mask = 2^n - 1
    indexes = Vector{Int}(undef, 2^n)
    for num=0:mask
        new = 0
        for i=1:n
            ni = perm[i]
            new += (num & (1 << (ni-1))) >> (ni - i)
        end
        indexes[num+1] = new+1
    end
    op[indexes, indexes]
end

"""
    ringexchangeoperator(n)

defined based on its action on Ising basis, which is a cricular shift
to right plus the hermitian conjugate.

With this definition heiseberg is 1/4 *(R(2)-I(4))
"""
function ringexchangeoperator(n::Int)
    @assert 1 < n < 7
    Pn = spzeros(2^n, 2^n)
    for j=1:2^n
        i = ((j-1) >> 1) | (((j-1) & 1) << (n-1))
        Pn[i+1, j] = 1
    end
    Pn + Pn'
end

# only for spin-half right now
function nbodyopexpansion(n::Int,
                          H::AbstractMatrix;
                          mode::Symbol=:SYMBOL,
                          verbose::Bool=false)
    1 < n < 7 || error("too large!")
    size(H) == (2^n, 2^n) || error("operator size doesn't match n")

    ops = [spinoperators(1/2)..., I(2)]
    ols = ['Z', 'P', 'M', 'I']
    symbols = [:Z, :P, :M, :I]
    T = eltype(H)
    symops = [spinoperators(1/2, symmetry=U1)..., eye(T, U1Space(0=>1, 1=>1))]

    amps = Vector{T}()
    if mode == :SYMBOL
        MType = Symbol
    elseif mode == :Trivial
        MType = Matrix{T}
    elseif mode == :U1
        MType = SymTensor{U1,T,2}
    else
        error("unrecognized mode!")
    end
    terms = Vector{QAmpTerm{MType, n, T}}()
    ⊗ = kron
    for is in Iterators.product([1:4 for i=1:n]...)
        op = reduce(⊗, [ops[i] for i in is], init=I(1))
        opnorm = norm(op)
        amp = dot(H, op)/norm(op)^2
        if amp != 0
            verbose && println([ols[i] for i in is]..., "  $amp")
            is == Tuple([4 for i=1:n]) &&
                @warn "Operator has identity with amplitude $amp"
            if mode == :SYMBOL
                term = [symbols[i] for i in is]
            elseif mode == :Trivial
                term = [Matrix{T}(ops[i]) for i in is]
            elseif mode == :U1
                term = [symops[i] for i in is]
            else
                error("unrecognized mode : $mode")
            end
            push!(terms, QAmpTerm(amp, Tuple(term)))
        end
    end
    terms
end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using QSimulator
using LinearAlgebra

include("unitary-tests.jl")
include("dissipation-tests.jl")

include("speed_test.jl")
#speed_test(5)

#include("test_dissipators.jl")


text:
module TestSimpleTypePrint

if isdefined(@__MODULE__, :LanguageServer)  # hack to make vscode linter work properly
    include("../src/SimpleTypePrint.jl")
    using .SimpleTypePrint
else
    using SimpleTypePrint: repr_type
end

using Test

module A
    module B
        struct Foo{A} end

        f(x) = x
    end
end

macro type_print_test(ty, result)
    esc(:(@test repr_type($ty) == $result))
end

@testset "type simple print" begin
    @test repr_type(Array{Int64, 2}) == "Array{Int64,2}"
    @test repr_type(Array{Int64}) == "Array{Int64,N} where N"
    @test repr_type(Tuple{A,B,C} where {A, B, C}) == "Tuple{A,B,C} where {A,B,C}"
    let nested = Tuple{Tuple{Tuple{A,B},C},D} where {A,B,C,D}
        expected = [
            "Tuple{...} where {A,B,C,D}",
            "Tuple{Tuple{...},D} where {A,B,C,D}",
            "Tuple{Tuple{Tuple{...},C},D} where {A,B,C,D}",
            "Tuple{Tuple{Tuple{A,B},C},D} where {A,B,C,D}",
            "Tuple{Tuple{Tuple{A,B},C},D} where {A,B,C,D}",
        ]
        @testset "with different depth limit" begin
            for d in 1:5
                @test repr_type(nested, max_depth=d) == expected[d]
            end
        end
    end

    @test repr_type(A.B.Foo) == "Foo{A} where A"
    @test endswith(repr_type(A.B.Foo; short_type_name=false), "A.B.Foo{A} where A")
    @test (repr_type(Tuple{A,B,C} where {A, B >: A, A<:C<:B})
            == "Tuple{A,B,C} where {A,B>:A,A<:C<:B}")
    @test (repr_type(Tuple{A} where {A <: Tuple{Tuple{Tuple{Int64}}}}, max_depth=3)
            == "Tuple{A} where A<:Tuple{Tuple{...}}")

    @type_print_test(Tuple{(Tuple{A} where A), A} where A,
        "Tuple{Tuple{A1} where A1,A} where A")

    union_s = repr_type(Union{A,B} where {A,B})
    @test union_s ∈ ["Union{A,B} where {A,B}", "Union{B,A} where {A,B}"]

    @type_print_test(StridedArray{T,N} where {T, N}, "StridedArray{T,N} where {T,N}")
    @type_print_test(StridedVecOrMat{String}, "StridedVecOrMat{String}")

    @type_print_test NTuple{2, Int64} "Tuple{Int64,Int64}"
    @type_print_test NTuple{5, Int64} "NTuple{5,Int64}"

    f_repr = repr_type(typeof(A.B.f))
    @test startswith(f_repr, "typeof(") && endswith(f_repr, "A.B.f)")
end

end  # module TestSimpleTypePrint

max_forks_repo_name: MrVPlusOne/SimpleTypePrint


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -3.428579430434188 <= q <= 8.991868736935128)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.758903986598098 + 0.8106379052679875*1/(1 + exp(-(0.47487042435595583 + 0.4488775132618512*1/(1 + exp(-(0.1716391672055586 + 0.7301849089692007*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + 0.6276077641275091*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))))) + 0.015398554931854491*1/(1 + exp(-(0.8366427516313624 + 0.8654420505988725*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + -0.07448940361379819*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3])))))))))) + -0.6668215058198093*1/(1 + exp(-(0.8216722410540171 + 0.5667551555087149*1/(1 + exp(-(0.1716391672055586 + 0.7301849089692007*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + 0.6276077641275091*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))))) + 0.5912184082121308*1/(1 + exp(-(0.8366427516313624 + 0.8654420505988725*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + -0.07448940361379819*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))))))))))) + 1/(1 + exp(-(-0.09143860592972564 + 0.3407854070498746*1/(1 + exp(-(0.47487042435595583 + 0.4488775132618512*1/(1 + exp(-(0.1716391672055586 + 0.7301849089692007*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + 0.6276077641275091*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))))) + 0.015398554931854491*1/(1 + exp(-(0.8366427516313624 + 0.8654420505988725*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + -0.07448940361379819*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3])))))))))) + 0.24612879529032217*1/(1 + exp(-(0.8216722410540171 + 0.5667551555087149*1/(1 + exp(-(0.1716391672055586 + 0.7301849089692007*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + 0.6276077641275091*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))))) + 0.5912184082121308*1/(1 + exp(-(0.8366427516313624 + 0.8654420505988725*1/(1 + exp(-(0.33202390613731003 + -0.8931954910636541*$(x[1]) + 0.7142724585602056*$(x[2]) + -0.28852147983722487*$(x[3])))) + -0.07448940361379819*1/(1 + exp(-(-0.8935890604396688 + -0.29936605585485676*$(x[1]) + -0.5391102865466948*$(x[2]) + 0.8633503800224367*$(x[3]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using debspkg
using Test

@testset "debspkg.jl" begin
    # Write your tests here.
    @test f(3,2)==25
    @test f(2,2)==16
    @test f(3,1)==16
    @test f(1,1)==4
    @test f(0,0)==0
    @test f(1,0)==1
    @test f(0,1)==1
end

@testset "derivative tests" begin
    @test derivative_of_f1(2,4)==12
    @test derivative_of_f1(2,3)==10
    @test derivative_of_f1(2,2)==8
end


text:
using ITensors, Test

a = [-0.1, -0.12]
@test ITensors.truncate!(a) == (0., 0.)
@test length(a) == 1
a = [0.1, 0.01, 1e-13]
@test ITensors.truncate!(a,absoluteCutoff=true,cutoff=1e-5) == (1e-13, (0.01 + 1e-13)/2)
@test length(a) == 2

i = Index(2,"i")
j = Index(2,"j")
A = randomITensor(i,j)
@test_throws ArgumentError factorize(A, i, dir="fakedir")

@testset "Spectrum" begin
  i = Index(100,"i")
  j = Index(100,"j")

  U,S,V = svd(rand(100,100))
  S ./= norm(S)
  A = ITensor(U*diagm(0=>S)*V', i,j)

  spec = svd(A,i).spec

  @test eigs(spec) ≈ S .^2
  @test truncerror(spec) == 0.0

  spec = svd(A,i; maxdim=length(S)-3).spec
  @test truncerror(spec) ≈ sum(S[end-2:end].^2)

end


text:
"""
Atomic solution for an DAE.

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `λ`: current solution of λ
* `λ̅`: previous solution of λ
* `v`: vector field of q
* `v̅`: vector field of q̅
* `u`: projective vector field of q
* `u̅`: projective vector field of q̅
"""
mutable struct AtomicSolutionDAE{DT,TT} <: AtomicSolution{DT,TT}
    t::TT
    t̅::TT

    q::Vector{DT}
    q̅::Vector{DT}
    q̃::Vector{DT}

    λ::Vector{DT}
    λ̅::Vector{DT}

    v::Vector{DT}
    v̅::Vector{DT}

    u::Vector{DT}
    u̅::Vector{DT}

    function AtomicSolutionDAE{DT, TT}(nd, nm) where {DT <: Number, TT <: Real}
        new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nm), zeros(DT, nm),
                                zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd))
    end
end

AtomicSolutionDAE(DT, TT, nd, nm) = AtomicSolutionDAE{DT, TT}(nd, nm)

function set_solution!(asol::AtomicSolutionDAE, sol)
    t, q, λ = sol
    asol.t  = t
    asol.q .= q
    asol.λ .= λ
    asol.v .= 0
end

function get_solution(asol::AtomicSolutionDAE)
    (asol.t, asol.q, asol.λ)
end

function CommonFunctions.reset!(asol::AtomicSolutionDAE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.λ̅ .= asol.λ
    asol.v̅ .= asol.v
    asol.u̅ .= asol.u
    asol.t += Δt
end

function update!(asol::AtomicSolutionDAE{DT}, y::Vector{DT}, λ::Vector{DT}) where {DT}
    for k in eachindex(y)
        update!(asol, y[k], k)
    end
    for k in eachindex(λ)
        asol.λ[k] = λ[k]
    end
end

function update!(asol::AtomicSolutionDAE{DT}, y::DT, k::Int) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(y, asol.q[k], asol.q̃[k])
end


text:
test = """
    0: 4 1 5
    1: 2 3 | 3 2
    2: 4 4 | 5 5
    3: 4 5 | 5 4
    4: "a"
    5: "b"

    ababbb
    bababa
    abbbab
    aaabbb
    aaaabbb
    """

function parse_input(input::AbstractString)
    (rules, messages) = split.(split(input, "\n\n"), "\n")
    rules = parse_rule.(rules)
    @assert sort(first.(rules)) == 0:length(rules) - 1
    return (sort(rules; by = first) .|> last, messages)
end

function parse_rule(rule::AbstractString)::Pair
    (i, body) = split(rule, ": ")
    rule = body[1] == '"' ? body[2] :
        [parse.(Int, split(subrule, " ")) .+ 1
         for subrule in split(body, " | ")]
    return parse(Int, i) => rule
end

function compile(rules::AbstractVector)::Regex
    re(i::Int) = re(rules[i])
    re(rule::Char) = rule
    re(rule::Vector) = "(?:" *
        join([join(re.(sub), "") for sub in rule], "|") * ")"
    re(f::Function) = f(re)

    return Regex("^$(re(1))\$")
end

function q1(rules, msgs)::Int
    re = compile(rules)
    return count(msg -> match(re, msg) !== nothing, msgs)
end

function q2(rules, msgs, a, b)::Int
    rules = copy(rules)
    rules[9] = re -> "(?:$(re(a))+)"
    rules[12] = re -> "(?<r11>$(re(a))\\g<r11>?$(re(b)))"
    return q1(rules, msgs)
end

@assert q1(parse_input(test)...) == 2
let input = parse_input(read("day19_test.in", String))
    @assert q1(input...) == 3
    @assert q2(input..., 31, 30) == 12
end

input = parse_input(read("day19.in", String))
println("Q1: ", q1(input...))
println("Q2: ", q2(input..., 43, 32))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
using GR
using LaTeXStrings

𝒩(μ, σ) = 1 / (σ * √(2π)) * exp.(-0.5 * ((x .- μ) / σ) .^ 2)

x = LinRange(-5, 5, 500);
y = hcat(𝒩(0, √0.2), 𝒩(0, √1), 𝒩(0, √5), 𝒩(-2, √0.5));

plot(x, y, xlabel=L"\mathcal{X}", ylabel=L"\mathcal{N}(\mu,\,\sigma^{2})", title=L"\frac{1}{\sigma\sqrt{2\pi}} e^{-\frac{1}{2} \left({\frac{x-\mu}{\sigma}}\right)^2}", labels=(L"\mu=0, \sigma^2=0.2", L"\mu=0, \sigma^2=1", L"\mu=0, \sigma^2=5", L"\mu={-2}, \sigma^2=0.5"), xlim=(-5.2, 5.2), ylim=(-0.05, 1.05), linewidth=3)


text:
mutable struct CreditCard
    name :: AbstractString
    prefixes :: AbstractString
    length :: Int64
    security_code :: AbstractString
    security_code_length :: Int64
end

    prefix_mastercard = ["51", "52", "53", "54", "55"]
    prefix_visa = ["4"]
    prefix_amex = ["34", "37"]
    prefix_discover = ["6011",]
    prefix_diners = ["300", "301", "302", "303", "304", "305"]
    prefix_jcb16 = ["3088", "3096", "3112", "3158", "3337", "3528"]
    prefix_jcb15 = ["2100", "1800"]
    prefix_voyager = ["8699"]

    credit_card_types = Dict(
        "mastercard"=>   CreditCard("Mastercard",      prefix_mastercard[rand(1:length(prefix_mastercard))], 16, "CVV",3),
        "visa16"=>       CreditCard("VISA 16 digit",   prefix_visa[rand(1:length(prefix_visa))], 16, "CVC",3),
        "visa13"=>       CreditCard("VISA 13 digit",   prefix_visa[rand(1:length(prefix_visa))] , 13, "CVC",3),
        "amex"=>         CreditCard("American Express",prefix_amex[rand(1:length(prefix_amex))], 15, "CID",4),
        "discover"=>     CreditCard("Discover",        prefix_discover[rand(1:length(prefix_discover))], 16, "CVC",3),
        "diners"=>       CreditCard("Diners Club / Carte Blanche", prefix_diners[rand(1:length(prefix_diners))] , 14, "CVC",3),
        "jcb15"=>        CreditCard("JCB 15 digit",    prefix_jcb15[rand(1:length(prefix_jcb15))] , 15, "CVC",3),
        "jcb16"=>        CreditCard("JCB 16 digit",    prefix_jcb16[rand(1:length(prefix_jcb16))], 16, "CVC",3),
        "voyager"=>      CreditCard("Voyager",         prefix_voyager[rand(1:length(prefix_voyager))] , 15, "CVC",3),
    )
    credit_card_types["visa"] = credit_card_types["visa16"]
    credit_card_types["jcb"] = credit_card_types["jcb16"]

    luhn_lookup = Dict("0"=> 0, "1"=> 2, "2"=> 4, "3"=> 6, "4"=> 8,
                   "5"=> 1, "6"=> 3, "7"=> 5, "8"=> 7, "9"=> 9)

    function credit_card_provider()
        executor(data["en"]["faker"]["business"]["credit_card_types"])
    end

    function credit_card_number(card="None")
         if card == "None"
            card = [key for key in  keys(credit_card_types)][rand(1:9)]
            card = credit_card_types[card]
        end
        prefix = card.prefixes
        number = generate_number(prefix, card.length)
        return number
    end

    #=function credit_card_expire(start="now", end="+10y", date_format="%m/%y")

        expire_date = DateTime.date_time_between(start, end)
        return expire_date.strftime(date_format)
    end=#

    function credit_card_full(card_types="None")
        if card_types == "None"
            card = [key for key in  keys(credit_card_types)][rand(1:9)]
        end
        card = credit_card_types[card]

        tpl = ("$(card.name) $(first_name()) $(last_name()) $(generate_number(card.prefixes,card.length)) $(card.security_code) $(credit_card_security_code())")

        tpl
    end

    function credit_card_security_code()
        return numerify("###")
    end

    function credit_card_type()

        card_type = [key for key in  keys(credit_card_types)][rand(1:9)]

        credit_card_types[card_type]
    end

    function generate_number(prefix, len)
        """
        "prefix" is the start of the CC number as a AbstractString, any number of digits.
        "len" is the len of the CC number to generate. Typically 13 or 16
        """
        number = prefix
        # Generate random char digits
        number *= "#"^(len - length(prefix) - 1)

        number = numerify(number)
        reverse = number[end:-1:1]
        # Calculate sum
        tot = 1
        pos = 1
        while pos < (len - 1)
            tot += luhn_lookup["$(reverse[pos])"]
            if pos != (len - 2)
                tot += Int(reverse[pos+1])
            end
        pos += 2
        end
        # Calculate check digit
        check_digit = (10 - (tot % 10)) % 10
        number *= "$(check_digit)"
        return number
    end


text:
export synsets, antonyms, hypernyms, hyponyms, expanded_hypernyms, expanded_hyponyms

synsets(db::DB, lemma::Lemma) = map(lemma.synset_offsets) do offset
    db.synsets[lemma.pos][offset]
end

relation(db::DB, synset::Synset, pointer_sym) = map(
    ptr -> db.synsets[ptr.pos][ptr.offset],
    filter(ptr -> ptr.sym == pointer_sym, synset.pointers)
)

function expanded_relation(db::DB, synset::Synset, pointer_sym)
     # Inititally include self, so as to avoid cycles,
     # but we will delete it at the end
    all_related = Set{Synset}([synset])
    queue = relation(db, synset, pointer_sym)

    while !isempty(queue)
        node = pop!(queue)
        node in all_related && continue  # done this one already
        push!(all_related, node)

        union!(queue, relation(db, node, pointer_sym))
    end

    delete!(all_related, synset)
    return all_related
end

antonyms(db::DB, synset::Synset)  = relation(db, synset, ANTONYM)
hyponyms(db::DB, synset::Synset)  = relation(db, synset, HYPONYM)
hypernyms(db::DB, synset::Synset) = relation(db, synset, HYPERNYM)

expanded_hyponyms(db::DB, synset::Synset) = expanded_relation(db, synset, HYPONYM)
expanded_hypernyms(db::DB, synset::Synset) = expanded_relation(db, synset, HYPERNYM)


text:
#### n-stream kernels
#  - delta-scaling
#  - adding two sets of properties
#  - extracting subsets
#  - validity checking

# ---------------------------------
# increment 1scalar by nstream
function increment_by_gpoint!(op_1::OneScalar{FT}, op_2::NStream{FT}) where {FT}
    #    integer,                             intent(in   ) :: ncol, nlay, ngpt
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau
    #    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa

    #    integer  :: icol, ilay, igpt

    op_1.tau = op_1.tau .+ op_2.tau * (FT(1) - op_2.ssa)
end

# ---------------------------------
# increment 2stream by nstream
function increment_by_gpoint!(op_1::TwoStream{FT}, op_2::NStream{FT}) where {FT}
    #    integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom2
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa, op_1.g
    #    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa
    #    real(FT), dimension(nmom2,
    #                        ncol,nlay,ngpt), intent(in   ) :: op_2.p

    #    integer  :: icol, ilay, igpt
    #    real(FT) :: tau12, tauscat12

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)

    for igpt in 1:ngpt
        for ilay in 1:nlay
            for icol in 1:ncol
                # t=op_1.tau + op_2.tau
                tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, igpt]
                # w=(op_1.tau*op_1.ssa + op_2.tau*op_2.ssa) / t
                tauscat12 =
                    op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] +
                    op_2.tau[icol, ilay, igpt] * op_2.ssa[icol, ilay, igpt]
                op_1.g[icol, ilay, igpt] =
                    (
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] * op_1.g[icol, ilay, igpt] +
                        op_2.tau[icol, ilay, igpt] * op_2.ssa[icol, ilay, igpt] * op_2.p[1, icol, ilay, igpt]
                    ) / max(ε_machine(FT), tauscat12)
                op_1.ssa[icol, ilay, igpt] = tauscat12 / max(ε_machine(FT), tau12)
                op_1.tau[icol, ilay, igpt] = tau12
            end
        end
    end
end

# increment nstream by 1scalar
function increment_by_gpoint!(op_1::NStream{FT}, op_2::OneScalar{FT}) where {FT}
    #    integer,                              intent(in   ) :: ncol, nlay, ngpt
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
    #    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau

    #    integer  :: icol, ilay, igpt
    #    real(FT) :: tau12

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)

    for igpt in 1:ngpt
        for ilay in 1:nlay
            for icol in 1:ncol
                tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, igpt]
                op_1.ssa[icol, ilay, igpt] =
                    op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] / max(ε_machine(FT), tau12)
                op_1.tau[icol, ilay, igpt] = tau12
                # p is unchanged
            end
        end
    end
end

# ---------------------------------
# increment nstream by 2stream
function increment_by_gpoint!(op_1::NStream{FT}, op_2::TwoStream{FT}) where {FT}
    #    integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom1
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
    #    real(FT), dimension(nmom1,
    #                        ncol,nlay,ngpt), intent(inout) :: op_1.p
    #    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa, op_2.g

    #    integer  :: icol, ilay, igpt
    #    real(FT) :: tau12, tauscat12
    #    real(FT), dimension(nmom1) :: temp_moms # TK
    #    integer  :: imom  #TK

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)

    for igpt in 1:ngpt
        for ilay in 1:nlay
            for icol in 1:ncol
                tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, igpt]
                tauscat12 =
                    op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] +
                    op_2.tau[icol, ilay, igpt] * op_2.ssa[icol, ilay, igpt]
                #
                # Here assume Henyey-Greenstein
                #
                temp_moms[1] = op_2.g[icol, ilay, igpt]
                for imom in 2:nmom1
                    temp_moms[imom] = temp_moms[imom - 1] * op_2.g[icol, ilay, igpt]
                end
                op_1.p[1:nmom1, icol, ilay, igpt] =
                    (
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] * op_1.p[1:nmom1, icol, ilay, igpt] +
                        op_2.tau[icol, ilay, igpt] * op_2.ssa[icol, ilay, igpt] * temp_moms[1:nmom1]
                    ) / max(ε_machine(FT), tauscat12)
                op_1.ssa[icol, ilay, igpt] = tauscat12 / max(ε_machine(FT), tau12)
                op_1.tau[icol, ilay, igpt] = tau12
            end
        end
    end
end

# ---------------------------------
# increment nstream by nstream
function increment_by_gpoint!(op_1::NStream{FT}, op_2::NStream{FT}) where {FT}
    #    integer,                              intent(in   ) :: ncol, nlay, ngpt, nmom1, nmom2
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
    #    real(FT), dimension(nmom1,
    #                        ncol,nlay,ngpt), intent(inout) :: op_1.p
    #    real(FT), dimension(ncol,nlay,ngpt), intent(in   ) :: op_2.tau, op_2.ssa
    #    real(FT), dimension(nmom2,
    #                        ncol,nlay,ngpt), intent(in   ) :: op_2.p

    #    integer  :: icol, ilay, igpt, mom_lim
    #    real(FT) :: tau12, tauscat12

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)

    mom_lim = min(nmom1, nmom2)
    for igpt in 1:ngpt
        for ilay in 1:nlay
            for icol in 1:ncol
                tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, igpt]
                tauscat12 =
                    op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] +
                    op_2.tau[icol, ilay, igpt] * op_2.ssa[icol, ilay, igpt]
                #
                # If op2 has more moments than op1 these are ignored;
                #   if it has fewer moments the higher orders are assumed to be 0
                #
                op_1.p[1:mom_lim, icol, ilay, igpt] =
                    (
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] * op_1.p[1:mom_lim, icol, ilay, igpt] +
                        op_2.tau[icol, ilay, igpt] * op_2.ssa[icol, ilay, igpt] * op_2.p[1:mom_lim, icol, ilay, igpt]
                    ) / max(ε_machine(FT), tauscat12)
                op_1.ssa[icol, ilay, igpt] = tauscat12 / max(eps, tau12)
                op_1.tau[icol, ilay, igpt] = tau12
            end
        end
    end
end

# increment 1scalar by nstream
function increment_bybnd!(op_1::OneScalar{FT}, op_2::NStream{FT}) where {FT}
    #    integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: tau1
    #    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: tau2, ssa2
    #    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

    #    integer :: ibnd, igpt

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)
    gpt_lims = get_band_lims_gpoint(op_1)
    nbnd = get_nband(op_1)

    for ibnd in 1:nbnd
        for igpt in gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
            op_1.tau[:, :, igpt] .= op_1.tau[:, :, igpt] .+ op_2.tau[:, :, ibnd] .* (FT(1) .- op_2.ssa[:, :, ibnd])
        end
    end
end

"""
    increment_bybnd!(op_1::TwoStream{FT}, op_2::NStream{FT}) where FT

integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom2, nbnd
real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa, op_1.g
real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau, op_2.ssa
real(FT), dimension(nmom2,
                    ncol,nlay,nbnd), intent(in   ) :: op_2.p
integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

integer  :: icol, ilay, igpt, ibnd
real(FT) :: tau12, tauscat12
"""
function increment_bybnd!(op_1::TwoStream{FT}, op_2::NStream{FT}) where {FT}

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)
    gpt_lims = get_band_lims_gpoint(op_1)
    nbnd = get_nband(op_1)

    for ibnd in 1:nbnd
        for igpt in gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
            for ilay in 1:nlay
                for icol in 1:ncol
                    # t=op_1.tau + op_2.tau
                    tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, ibnd]
                    # w=(op_1.tau*op_1.ssa + op_2.tau*op_2.ssa) / t
                    tauscat12 =
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] +
                        op_2.tau[icol, ilay, ibnd] * op_2.ssa[icol, ilay, ibnd]
                    op_1.g[icol, ilay, igpt] =
                        (
                            op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] * op_1.g[icol, ilay, igpt] +
                            op_2.tau[icol, ilay, ibnd] * op_2.ssa[icol, ilay, ibnd] * op_2.p[1, icol, ilay, ibnd]
                        ) / max(ε_machine(FT), tauscat12)
                    op_1.ssa[icol, ilay, igpt] = tauscat12 / max(ε_machine(FT), tau12)
                    op_1.tau[icol, ilay, igpt] = tau12
                end
            end
        end
    end
end

# increment nstream by 1scalar
function increment_bybnd!(op_1::NStream{FT}, op_2::OneScalar{FT}) where {FT}
    #    integer,                             intent(in   ) :: ncol, nlay, ngpt, nbnd
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
    #    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau
    #    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

    #    integer  :: icol, ilay, igpt, ibnd
    #    real(FT) :: tau12

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)
    gpt_lims = get_band_lims_gpoint(op_1)
    nbnd = get_nband(op_1)

    for ibnd in 1:nbnd
        for igpt in gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
            for ilay in 1:nlay
                for icol in 1:ncol
                    tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, ibnd]
                    op_1.ssa[icol, ilay, igpt] =
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] / max(ε_machine(FT), tau12)
                    op_1.tau[icol, ilay, igpt] = tau12
                    # p is unchanged
                end
            end
        end
    end
end

# ---------------------------------
# increment nstream by 2stream
function increment_bybnd!(op_1::NStream{FT}, op_2::TwoStream{FT}) where {FT}
    #    integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom1, nbnd
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
    #    real(FT), dimension(nmom1,
    #                        ncol,nlay,ngpt), intent(inout) :: op_1.p
    #    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau, op_2.ssa, op_2.g
    #    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

    #    integer  :: icol, ilay, igpt, ibnd
    #    real(FT) :: tau12, tauscat12
    #    real(FT), dimension(nmom1) :: temp_moms # TK
    #    integer  :: imom  #TK

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)
    gpt_lims = get_band_lims_gpoint(op_1)
    nbnd = get_nband(op_1)

    for ibnd in 1:nbnd
        for igpt in gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
            for ilay in 1:nlay
                for icol in 1:ncol
                    tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, ibnd]
                    tauscat12 =
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] +
                        op_2.tau[icol, ilay, ibnd] * op_2.ssa[icol, ilay, ibnd]
                    #
                    # Here assume Henyey-Greenstein
                    #
                    temp_moms[1] = op_2.g[icol, ilay, ibnd]
                    for imom in 2:nmom1
                        temp_moms[imom] = temp_moms[imom - 1] * op_2.g[icol, ilay, ibnd]
                    end
                    op_1.p[1:nmom1, icol, ilay, igpt] =
                        (
                            op_1.tau[icol, ilay, igpt] *
                            op_1.ssa[icol, ilay, igpt] *
                            op_1.p[1:nmom1, icol, ilay, igpt] +
                            op_2.tau[icol, ilay, ibnd] * op_2.ssa[icol, ilay, ibnd] * temp_moms[1:nmom1]
                        ) / max(ε_machine(FT), tauscat12)
                    op_1.ssa[icol, ilay, igpt] = tauscat12 / max(ε_machine(FT), tau12)
                    op_1.tau[icol, ilay, igpt] = tau12
                end
            end
        end
    end
end

# ---------------------------------
# increment nstream by nstream
function increment_bybnd!(op_1::NStream{FT}, op_2::NStream{FT}) where {FT}
    #    integer,                             intent(in   ) :: ncol, nlay, ngpt, nmom1, nmom2, nbnd
    #    real(FT), dimension(ncol,nlay,ngpt), intent(inout) :: op_1.tau, op_1.ssa
    #    real(FT), dimension(nmom1,
    #                        ncol,nlay,ngpt), intent(inout) :: op_1.p
    #    real(FT), dimension(ncol,nlay,nbnd), intent(in   ) :: op_2.tau, op_2.ssa
    #    real(FT), dimension(nmom2,
    #                        ncol,nlay,nbnd), intent(in   ) :: op_2.p
    #    integer,  dimension(2,nbnd),         intent(in   ) :: gpt_lims # Starting and ending gpoint for each band

    #    integer  :: icol, ilay, igpt, ibnd, mom_lim
    #    real(FT) :: tau12, tauscat12

    ngpt = get_ngpt(op_1)
    nlay = get_nlay(op_1)
    ncol = get_ncol(op_1)
    gpt_lims = get_band_lims_gpoint(op_1)
    nbnd = get_nband(op_1)

    mom_lim = min(nmom1, nmom2)
    for ibnd in 1:nbnd
        for igpt in gpt_lims[1, ibnd]:gpt_lims[2, ibnd]
            for ilay in 1:nlay
                for icol in 1:ncol
                    tau12 = op_1.tau[icol, ilay, igpt] + op_2.tau[icol, ilay, ibnd]
                    tauscat12 =
                        op_1.tau[icol, ilay, igpt] * op_1.ssa[icol, ilay, igpt] +
                        op_2.tau[icol, ilay, ibnd] * op_2.ssa[icol, ilay, ibnd]
                    #
                    # If op2 has more moments than op1 these are ignored;
                    #   if it has fewer moments the higher orders are assumed to be 0
                    #
                    op_1.p[1:mom_lim, icol, ilay, igpt] =
                        (
                            op_1.tau[icol, ilay, igpt] *
                            op_1.ssa[icol, ilay, igpt] *
                            op_1.p[1:mom_lim, icol, ilay, igpt] +
                            op_2.tau[icol, ilay, ibnd] *
                            op_2.ssa[icol, ilay, ibnd] *
                            op_2.p[1:mom_lim, icol, ilay, ibnd]
                        ) / max(ε_machine(FT), tauscat12)
                    op_1.ssa[icol, ilay, igpt] = tauscat12 / max(eps, tau12)
                    op_1.tau[icol, ilay, igpt] = tau12
                end
            end
        end
    end
end

max_stars_repo_name: climate-machine/JRRTMGP

max_issues_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
using DataDeps
using QuantumStateBase

@testset "gen_data" begin
    n = 2

    points, 𝛒s, args, σs = gen_data_sqth_th(n_data=n)
    @test all([
        size(points) == (2, 4096, n),
        size(𝛒s) == (100, 100, n),
        size(args) == (6, n)
    ])

    points, 𝛒s, args, σs = gen_data_sqth(n_data=n)
    @test all([
        size(points) == (2, 4096, n),
        size(𝛒s) == (100, 100, n),
        size(args) == (3, n)
    ])
end


text:
#!/usr/bin/env julia -

using Luxor

using Test

using Random
Random.seed!(42)

function polysampletest(fname)
    Drawing(400, 400, fname)
    origin()
    background("white")

    p = ngon(O, 100, 6, vertices=true)
    setline(12)
    sethue("grey70")
    prettypoly(p, :stroke, close=true, ()->circle(O, 10, :fill))

    upsample = polysample(p, 12)
    setline(3)
    sethue("cyan")
    prettypoly(upsample, :stroke, close=true, ()->circle(O, 7, :fill))

    resample = polysample(p, 6)
    setline(2)
    sethue("magenta")
    prettypoly(upsample, :stroke, close=true, ()->circle(O, 4, :fill))

    @test p == circshift(resample, 1)

    downsample = polysample(p, 3)
    setline(1)
    sethue("orange")
    prettypoly(downsample, :stroke, close=true, ()->circle(O, 2, :fill))

    @test length(findall(pt -> pt == downsample[1], p)) == 1

    @test finish() == true
    println("...finished test: output in $(fname)")

end

polysampletest("polysample.svg")


text:
using MakieGallery, AbstractPlotting, GLMakie, Test, Makie
using MakieGallery: @block, @cell

push!(MakieGallery.plotting_backends, "Makie")
database = MakieGallery.load_database()
tested_diff_path = joinpath(@__DIR__, "tested_different")
test_record_path = joinpath(@__DIR__, "test_recordings")
for path in (tested_diff_path, test_record_path)
    rm(path, force = true, recursive = true)
    mkpath(path)
end
recordings = MakieGallery.record_examples(test_record_path)
@test length(recordings) == length(database)
MakieGallery.run_comparison(test_record_path, tested_diff_path)


text:
export SumFunction

struct SumFunction{T, N} <: AbstractMathFun{T, N}
    terms::MFCollection

    SumFunction{T, N}() where {T, N} = new{T, N}(MFCollection(undef, 0))
    SumFunction{T, N}(terms::MFCollection) where {T, N} = new{T, N}(terms)

    Base.insert!(p::SumFunction, i::Integer, term::AbstractMathFun) = insert!(p.terms, i, term)
    Base.push!(p::SumFunction, term::AbstractMathFun) = push!(p.terms, term)
    Base.deleteat!(p::SumFunction, range) = deleteat!(p.terms, range)
    Base.length(p::SumFunction) = length(p.terms)
    Base.getindex(p::SumFunction, i::Integer) = getindex(p.terms, i)
    Base.setindex!(p::SumFunction, val::AbstractMathFun, i::Integer) = setindex!(p.terms, val, i)
    Base.iterate(p::SumFunction) = Base.iterate(p.terms)
    Base.iterate(p::SumFunction, state) = Base.iterate(p.terms, state)

    function Functions.string(p::SumFunction, emitter::FunctionEmitter)
        str = "∑("
        first_print = true
        for term in p.terms
            term_str = string(term, emitter)
            length(term_str) != 0 && (str *= (first_print ? term_str : ", " * term_str); first_print = false)
        end
        str * ")"
    end

    function (p::SumFunction)(call::MFCall)
        Sum = 0.0
        for term in p.terms
            Sum += term(call)
        end
        call(Sum)
    end

    function Functions.emit(f::SumFunction{T, N}, emitter::FunctionEmitter{T}) where {T, N}
        sum_expr = Expr(:call, :∑)
        for term in f.terms
            push!(sum_expr.args, emit(term, emitter))
        end
        push!(sum_expr.args, Expr(:kw, :T, Symbol(string(T))))
        sum_expr
    end

    function Functions.simplify(f::SumFunction{T, N}, emit_constants::Bool) where {T, N}
        constsum = 0
        termcpy = MFCollection(undef, 0)

        for item in f.terms
            item = simplify(item, emit_constants)
            if is_constant(item, emit_constants)
                constsum += item.value
            else
                if item isa SumFunction
                    append!(termcpy, item.terms)
                else
                    push!(termcpy, item)
                end
            end
        end

        (constsum != 0) && push!(termcpy,  ConstFunction{T}(constsum))
        (length(termcpy) == 0) && return ConstFunction{T}(constsum)
        length(termcpy) == 1 && return termcpy[1]
        return SumFunction{T, N}(termcpy)
    end
end


text:
module CellFieldsTests

#This is a first steps to decouple CellFields from the Geometry

using Test
using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.Fields
using Gridap.Fields: IntKernel
using Gridap.Fields: FieldOpKernel
using Gridap.Fields: trialize_basis_value
using Gridap.TensorValues
using Gridap.Fields: MockField, MockBasis, OtherMockBasis

# In this furure all this shoul be exported from Fields
using Gridap.Geometry: GenericCellField, get_cell_axes
using Gridap.Geometry: convert_to_cell_field, get_metasize
using Gridap.Geometry: CellField, similar_object, trialize_cell_basis
using Gridap.Geometry: merge_skeleton_dof_ids, merge_skeleton_cell_fields
using Gridap.Geometry: is_basis, is_test, is_trial

np = 3
ndofs = 4

p = Point(1,2)
x = fill(p,np)
z = 2.0

v = VectorValue(3.0,1.5)
w = VectorValue(3.4,3.5)
a = MockBasis{2}(v,ndofs)
b = MockBasis{2}(w,ndofs)
c = fill(1.0,ndofs)
f = OtherMockBasis{2}(ndofs)

g = MockField{2}(v)

l = 10
xl = Fill(x,l)
zl = [ z for  i in 1:l]
cl = fill(c,l)
fl = Fill(f,l)
ϕl = lincomb(fl,cl)
gl = fill(g,l)
al = Fill(a,l)
bl = fill(b,l)

gf = GenericCellField(gl,ϕl,Val(true))
af = GenericCellField(al,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
bf = GenericCellField(bl,ϕl,Val(true),Fill((Base.OneTo(ndofs),),l),Val((:,)))
zf = convert_to_cell_field(zl,ϕl)
df = af*zf
dft = trialize_cell_basis(df)

# Check memoization
df_x1 = evaluate(df,xl)
df_x2 = evaluate(df,xl)
@test df_x1 === df_x2
∇gf1 = ∇(gf)
∇gf2 = ∇(gf)
@test ∇gf1 === ∇gf2
@test evaluate(∇gf1,xl) === evaluate(∇gf2,xl)
εgf1 = ε(gf)
εgf2 = ε(gf)
@test εgf1 === εgf2
@test ∇×gf === ∇×gf
@test evaluate(∇×gf,xl) === evaluate(∇×gf,xl)

@test is_test(af)
@test is_trial(dft)
@test is_basis(af)
@test is_basis(dft)
mf = af⋅dft
@test get_metasize(mf) == (:,:)
mf_x = evaluate(mf,xl)
@test size(mf_x[1]) == (np,ndofs,ndofs)

@test get_cell_axes(mf) == Fill((Base.OneTo(ndofs),Base.OneTo(ndofs)),l)

idsL = [ i*collect(1:ndofs)  for i in 1:l]
idsR = [ 2*i*collect(1:ndofs)  for i in 1:l]
axesL = Fill((Base.OneTo(ndofs),),l)
axesR = Fill((Base.OneTo(ndofs),),l)

idsS = merge_skeleton_dof_ids(idsL,idsR,axesL,axesR)
@test isa(idsS,VectorOfBlockArrayCoo)

afL, afR = merge_skeleton_cell_fields(af,2*af)

afL_x = evaluate(afL,xl)
afR_x = evaluate(afR,xl)
@test isa(afL_x,VectorOfBlockArrayCoo)
@test isa(afR_x,VectorOfBlockArrayCoo)

end # module


text:
using Pkg
Pkg.instantiate()
Pkg.develop(PackageSpec(path=pwd()))
Pkg.build("JuSwalbe")

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
module LossFactory

function

end


text:
## ------------------------------------------------------------------
# USER CONFIG
_rbook_toml_name() = ".RemAssistant.toml"
config_file() = joinpath(homedir(), _rbook_toml_name())

function read_configfile()
    cfile = config_file()
    isfile(cfile) ? TOML.parsefile(cfile) : Dict()
end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(CUDArt.init, (Array{Int64, 1},))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Float64, 2}, Int64,))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{UInt16, 2}, Int64,))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Int64, 1}, Int64,))
    precompile(CUDArt.attribute, (Int64, UInt32,))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{Int64, 1}, Int64,))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Int32, 2}, Int64,))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{Float64, 2}, Int64,))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Int32, 2},))
    precompile(CUDArt._copy!, (Array{Float64, 2}, CUDArt.CudaArray{Float64, 2}, CUDArt.Stream,))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{UInt16, 2}, Int64,))
    precompile(CUDArt._copy!, (Array{UInt16, 2}, CUDArt.CudaArray{UInt16, 2}, CUDArt.Stream,))
    precompile(CUDArt._copy!, (Array{Int32, 2}, CUDArt.CudaArray{Int32, 2}, CUDArt.Stream,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray{Int64, 1}}, CUDArt.CudaPtr{Int64}, Tuple{Int64}, Int64,))
    precompile(CUDArt._copy!, (Array{Int64, 1}, CUDArt.CudaArray{Int64, 1}, CUDArt.Stream,))
    precompile(CUDArt._copy!, (CUDArt.HostArray{Int32, 2}, CUDArt.CudaArray{Int32, 2}, CUDArt.Stream,))
    precompile(CUDArt.pitchedptr, (Array{Int32, 2},))
    precompile(CUDArt.reinterpret, (Type{Float64}, CUDArt.CudaArray{Float64, 1}, Tuple{Int64, Int64, Int64, Int64},))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Float32, 1},))
    precompile(CUDArt.pitchedptr, (CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt._copy!, (Array{Float64, 1}, CUDArt.CudaArray{Float64, 1}, CUDArt.Stream,))
    precompile(CUDArt.CUDArt_gen.cudaMalloc, (Base.RefValue{Ptr{Void}}, Int64,))
    precompile(CUDArt.pitchedptr, (Array{UInt16, 2},))
    precompile(CUDArt.stride, (CUDArt.CudaPitchedArray{Int32, 2}, Int64,))
    precompile(CUDArt.pitchedptr, (Array{Float64, 2},))
    precompile(CUDArt._copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Base.SubArray{Float32, 2, Array{Float32, 2}, Tuple{Base.StepRange{Int64, Int64}, Base.Colon}, false}, CUDArt.Stream,))
    precompile(CUDArt._copy!, (CUDArt.CudaArray{Int32, 2}, CUDArt.HostArray{Int32, 2}, CUDArt.Stream,))
    precompile(CUDArt.pitchedptr, (Array{Int64, 1},))
    precompile(CUDArt._copy!, (CUDArt.CudaArray{Float64, 2}, Array{Float64, 2}, CUDArt.Stream,))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt._copy!, (CUDArt.CudaArray{UInt16, 2}, Array{UInt16, 2}, CUDArt.Stream,))
    precompile(CUDArt._copy!, (CUDArt.CudaArray{Int64, 1}, Array{Int64, 1}, CUDArt.Stream,))
    precompile(CUDArt._copy!, (CUDArt.CudaArray{Int32, 2}, Array{Int32, 2}, CUDArt.Stream,))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.device, ())
    precompile(CUDArt._copy!, (CUDArt.CudaArray{Float64, 1}, Array{Float64, 1}, CUDArt.Stream,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray}, Base.UnitRange{Int64},))
    precompile(CUDArt.close, (Array{Int64, 1},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray}, Array{Float64, 2},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPtr{UInt8}}, Ptr{UInt8}, CUDAdrv.CuContext,))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Array{Float32, 2},))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{Int32, 2},))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt.size, (CUDArt.CudaArray{UInt16, 2}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Float64, 2}, Float64,))
    precompile(CUDArt.checkstrides_pitched, (CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt.similar, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.similar, (CUDArt.CudaArray{Int64, 1}, Tuple{Int64},))
    precompile(CUDArt.get!, (CUDArt.CudaPitchedArray{Float64, 2}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Int64, Base.UnitRange{Int64}}, Float64,))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Float32},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPitchedArray}, Base.UnitRange{Int64},))
    precompile(CUDArt.similar, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.checkstrides_pitched, (Array{Float64, 2},))
    precompile(CUDArt.Type, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{UInt16}, Tuple{Int64, Int64},))
    precompile(CUDArt.CUDArt_gen.checkerror, (UInt32,))
    precompile(CUDArt.Type, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Int64}, Tuple{Int64},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Int32},))
    precompile(CUDArt.Type, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Int32}, Tuple{Int64, Int64},))
    precompile(CUDArt.Type, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Float64}, Tuple{Int64, Int64},))
    precompile(CUDArt.device, (Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CUDArt_gen.cudaExtent}, Type{Float32}, Tuple{Int64},))
    precompile(CUDArt.vec, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.free, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Int64},))
    precompile(CUDArt._copy!, (Array{Int32, 2}, CUDArt.CudaPitchedArray{Int32, 2}, CUDArt.Stream,))
    precompile(CUDArt.pointer, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{Float64},))
    precompile(CUDArt.free, (CUDArt.CudaPtr{UInt16},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPitchedArray{Float64, 2}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.pointer, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.size, (CUDArt.CudaArray{Float64, 4},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray{Int32, 2}}, CUDArt.CudaPtr{Int32}, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.copy, (CUDArt.CudaPtr{Float64},))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Int64, 1}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray{Float64, 1}}, CUDArt.CudaPtr{Float64}, Tuple{Int64}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Int64, 1}, UInt8,))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Int32, 2}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Int32, 2}, Int64,))
    precompile(CUDArt.devices, (Function,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray{UInt16, 2}}, CUDArt.CudaPtr{UInt16}, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPitchedArray{Int32, 2}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray{Float32, 1}}, CUDArt.CudaPtr{Float32}, Tuple{Int64}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Float64, 2}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPitchedArray{Float32, 1}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64}, Int64,))
    precompile(CUDArt.free, (CUDArt.CudaPtr{UInt8},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPitchedArray{UInt16, 2}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPitchedArray{Int64, 1}}, CUDArt.CUDArt_gen.cudaPitchedPtr, Tuple{Int64}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaArray{Float64, 2}}, CUDArt.CudaPtr{Float64}, Tuple{Int64, Int64}, Int64,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Int64, 1}, Int64,))
    precompile(CUDArt.size, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Int64, 1}, UInt8,))
    precompile(CUDArt.copy!, (Array{Int32, 2}, CUDArt.CudaArray{Int32, 2},))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPtr{Float32}}, Ptr{Float32}, CUDAdrv.CuContext,))
    precompile(CUDArt.fill!, (CUDArt.CudaArray{Float64, 2}, Float64,))
    precompile(CUDArt.similar, (CUDArt.CudaPitchedArray{Int64, 1}, Tuple{Int64},))
    precompile(CUDArt.fill!, (CUDArt.CudaPitchedArray{Float64, 2}, Int64,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPtr{Int32}}, Ptr{Int32}, CUDAdrv.CuContext,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPtr{Int64}}, Ptr{Int64}, CUDAdrv.CuContext,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPtr{UInt16}}, Ptr{UInt16}, CUDAdrv.CuContext,))
    precompile(CUDArt.Type, (Type{CUDArt.CudaPtr{Float64}}, Ptr{Float64}, CUDAdrv.CuContext,))
    precompile(CUDArt.size, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt.size, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.init, (Int64,))
    precompile(CUDArt.checkstrides_pitched, (Array{UInt16, 2},))
    precompile(CUDArt.checkstrides_pitched, (Array{Int32, 2},))
    precompile(CUDArt.pointer, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.vec, (CUDArt.CudaArray{Float64, 1},))
    precompile(CUDArt.checkstrides_pitched, (Array{Int64, 1},))
    precompile(CUDArt.device_reset, (Int64,))
    precompile(CUDArt.free, (CUDArt.HostArray{Int32, 2},))
    precompile(CUDArt.close, (Int64,))
    precompile(CUDArt._copy!, (CUDArt.CudaArray{Float64, 2}, Base.SubArray{Float32, 2, Array{Float32, 2}, Tuple{Base.StepRange{Int64, Int64}, Base.Colon}, false}, CUDArt.Stream,))
    precompile(CUDArt.copy!, (CUDArt.CudaArray{Float64, 2}, Array{Float32, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}}, CUDArt.CudaPitchedArray{Float64, 2}, Tuple{Base.UnitRange{Int64}, Base.UnitRange{Int64}},))
    precompile(CUDArt.devcount, ())
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{Int64, 1},))
    precompile(CUDArt._copy!, (CUDArt.CudaPitchedArray{UInt16, 2}, Array{UInt16, 2}, CUDArt.Stream,))
    precompile(CUDArt._copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Array{Float64, 2}, CUDArt.Stream,))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{UInt16, 2}, CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt._copy!, (CUDArt.CudaPitchedArray{Int64, 1}, Array{Int64, 1}, CUDArt.Stream,))
    precompile(CUDArt.copy!, (CUDArt.CudaArray{UInt16, 2}, CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{Float64, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{UInt16, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaPitchedArray{Float64, 2}, Base.SubArray{Float32, 2, Array{Float32, 2}, Tuple{Base.StepRange{Int64, Int64}, Base.Colon}, false},))
    precompile(CUDArt._copy!, (CUDArt.CudaPitchedArray{Int32, 2}, Array{Int32, 2}, CUDArt.Stream,))
    precompile(CUDArt.to_host, (CUDArt.CudaPitchedArray{Int32, 2},))
    precompile(CUDArt.copy!, (CUDArt.CudaArray{Float64, 2}, Base.SubArray{Float32, 2, Array{Float32, 2}, Tuple{Base.StepRange{Int64, Int64}, Base.Colon}, false},))
    precompile(CUDArt.copy!, (CUDArt.HostArray{Int32, 2}, CUDArt.CudaArray{Int32, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Float64, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Int32, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Int64, 1},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{UInt16, 2},))
    precompile(CUDArt.to_host, (CUDArt.CudaArray{Float64, 1},))
    precompile(CUDArt.free, (CUDArt.CudaArray{Float64, 2},))
    nothing
end

_precompile_()


text:
using Base: @kwdef

abstract type SOperator end

@kwdef mutable struct SDecorrelator
    stepval::Int = 16
    rngrange::Int = 255
    buffer::CircularBuffer{Bool} = CircularBuffer{Bool}(2)
    counter::Int = 0
    rng::MersenneTwister = MersenneTwister(rand(UInt))

    function SDecorrelator(stepval, rngrange, buffer, counter, rng)
        decorr = new(stepval, rngrange, buffer, counter, rng)
        fill!(decorr.buffer, false)

        return decorr
    end
end
function (op::SDecorrelator)(x::Bool)
    # Increment counters
    op.counter += op.stepval * x

    # Decide output
    z = pop!(op.buffer)
    r = rand(op.rng, 0:op.rngrange)
    if r <= op.counter
        push!(op.buffer, true)
    else
        push!(op.buffer, false)
    end

    # Decrement counters
    op.counter = max(op.counter - op.stepval * z, 0)

    return z
end

"""
    SSignedDecorrelator

A stochastic bitstream decorrelator.
"""
@kwdef struct SSignedDecorrelator <: SOperator
    pdecorr::SDecorrelator = SDecorrelator()
    ndecorr::SDecorrelator = SDecorrelator()
end
(op::SSignedDecorrelator)(x::SBit) = (op.pdecorr(pos(x)), op.ndecorr(neg(x)))

@kwdef mutable struct SAdder
    counter::Int = 0
end
function (op::SAdder)(x::Bool, y::Bool)
    # Increment counter
    op.counter += x + y

    # Decide output
    z = (op.counter >= 1)

    # Decrement counter
    op.counter = max(op.counter - z, 0)

    return z
end

"""
    SSignedAdder

A signed stochastic bitstream add operator.
"""
@kwdef struct SSignedAdder <: SOperator
    padder::SAdder = SAdder()
    nadder::SAdder = SAdder()
end
function (op::SSignedAdder)(x::SBit, y::SBit)
    pbit = op.padder(pos(x), pos(y))
    nbit = op.nadder(neg(x), neg(y))

    return (pbit, nbit)
end

@kwdef mutable struct SSaturatingSubtractor
    counter::Int = 0
end
function (op::SSaturatingSubtractor)(x::Bool, y::Bool)
    # # Increment counter
    # op.counter += (xor(x, y) && y) - (xor(x, y) && x)

    # # Saturate counters
    # op.counter = max(op.counter, 0)

    # # Decide output
    # z = (op.counter == 0 && x && !y)

    op.counter += x - y

    z = (op.counter >= 1)

    op.counter = max(op.counter - z, 0)

    return z
end

"""
    SSignedSubtractor

A signed stochastic bitstream subtract operator.
"""
@kwdef struct SSignedSubtractor <: SOperator
    padder::SAdder = SAdder()
    nadder::SAdder = SAdder()
    ppsub::SSaturatingSubtractor = SSaturatingSubtractor()
    pnsub::SSaturatingSubtractor = SSaturatingSubtractor()
    npsub::SSaturatingSubtractor = SSaturatingSubtractor()
    nnsub::SSaturatingSubtractor = SSaturatingSubtractor()
end
function (op::SSignedSubtractor)(x::SBit, y::SBit)
    pp = op.ppsub(pos(x), pos(y))
    pn = op.pnsub(pos(y), pos(x))
    np = op.npsub(neg(y), neg(x))
    nn = op.nnsub(neg(x), neg(y))
    z = (op.padder(pp, np), op.nadder(pn, nn))

    return z
end

struct SMultiplier end
(op::SMultiplier)(x::Bool, y::Bool) = x & y

"""
    SSignedMultiplier

A signed stochastic bitstream multiply operator.
"""
@kwdef struct SSignedMultiplier <: SOperator
    ppmult::SMultiplier = SMultiplier()
    pnmult::SMultiplier = SMultiplier()
    npmult::SMultiplier = SMultiplier()
    nnmult::SMultiplier = SMultiplier()
    sub11::SSaturatingSubtractor = SSaturatingSubtractor()
    sub12::SSaturatingSubtractor = SSaturatingSubtractor()
    sub13::SSaturatingSubtractor = SSaturatingSubtractor()
    sub14::SSaturatingSubtractor = SSaturatingSubtractor()
    padder::SAdder = SAdder()
    nadder::SAdder = SAdder()
end
function (op::SSignedMultiplier)(x::SBit, y::SBit)
    pp = op.ppmult(pos(x), pos(y))
    pn = op.pnmult(pos(x), neg(y))
    np = op.npmult(neg(x), pos(y))
    nn = op.nnmult(neg(x), neg(y))

    s₁₁ = op.sub11(pp, pn)
    s₁₂ = op.sub12(pn, pp)
    s₁₃ = op.sub13(np, nn)
    s₁₄ = op.sub14(nn, np)

    z = (op.padder(s₁₁, s₁₄), op.nadder(s₁₂, s₁₃))

    return z
end

@kwdef mutable struct SDivider
    counter::Int = 0
    zand::Bool = false
    rng::MersenneTwister = MersenneTwister(rand(UInt))
end
function (op::SDivider)(x::Bool, y::Bool)
    # Update counter
    op.counter = max(op.counter + 2 * x - 2 * op.zand, -100)

    # Decide output
    r = rand(op.rng, 0:64)
    z = (op.counter > r)

    # Update zand
    op.zand = z && y

    return z
end

"""
    SSignedDivider

A signed stochastic bitstream divide operator.
"""
@kwdef struct SSignedDivider <: SOperator
    pdiv::SDivider = SDivider()
    ndiv::SDivider = SDivider()
    psub::SSaturatingSubtractor = SSaturatingSubtractor()
    nsub::SSaturatingSubtractor = SSaturatingSubtractor()
end
function (op::SSignedDivider)(x::SBit, y::SBit)
    pp = op.pdiv(pos(x), pos(y))
    np = op.ndiv(neg(x), pos(y))

    z = (op.psub(pp, np), op.nsub(np, pp))

    return z
end

@kwdef mutable struct SFixedGainDivider
    counter::Int = 0
end
function (op::SFixedGainDivider)(x::Bool, y::Real)
    # Update counter
    op.counter += 255 * x

    # Decide output
    z = (op.counter >= round(255 * y))

    # Decrement counter
    op.counter -= z * round(255 * y)

    return z
end

"""
    SSignedFixedGainDivider

A stochastic bitstream fixed gain divide operator.
"""
@kwdef struct SSignedFixedGainDivider <: SOperator
    pdiv::SFixedGainDivider = SFixedGainDivider()
    ndiv::SFixedGainDivider = SFixedGainDivider()
end
(op::SSignedFixedGainDivider)(x::SBit, y::Real) = (op.pdiv(pos(x), y), op.ndiv(neg(x), y))

"""
    SSquareRoot

A stochastic bitstream square root operator.
"""
@kwdef mutable struct SSquareRoot <: SOperator
    counter::Int = 0
    zand::Bool = false
    buffer::CircularBuffer{Bool} = CircularBuffer{Bool}(1)
    rng::MersenneTwister = MersenneTwister(rand(UInt))

    function SSquareRoot(counter, zand, buffer, rng)
        op = new(counter, zand, buffer, rng)
        fill!(op.buffer, false)

        return op
    end
end
function (op::SSquareRoot)(x::SBit)
    # Update counter
    op.counter = max(op.counter + 4 * pos(x) - 4 * op.zand, -100)

    # Decide output
    r = rand(op.rng, 0:511)
    z = (op.counter >= r)

    # Update zand
    op.zand = z && pop!(op.buffer)
    push!(op.buffer, z)

    return (z, false)
end

struct SMatMultiplier
    counters::Array{Int, 2}
end
SMatMultiplier(nrows, ncols) = SMatMultiplier(zeros(nrows, ncols))
function (op::SMatMultiplier)(x::VecOrMat{Bool}, y::VecOrMat{Bool})
    # Increment counters
    op.counters .+= x * y

    # Decide output
    z = (op.counters .>= 1)

    # Decrement counters
    op.counters .-= z

    return z
end
(op::SMatMultiplier)(x::BitArray, y::BitArray) =
    op(convert(Array{Bool, ndims(x)}, x), convert(Array{Bool, ndims(y)}, y))

"""
    SSignedMatMultiplier

A stochastic bitstream matrix multiply operator.
"""
struct SSignedMatMultiplier <: SOperator
    ppmult::SMatMultiplier
    pnmult::SMatMultiplier
    npmult::SMatMultiplier
    nnmult::SMatMultiplier
    sub1::Array{SSaturatingSubtractor, 2}
    sub2::Array{SSaturatingSubtractor, 2}
    sub3::Array{SSaturatingSubtractor, 2}
    sub4::Array{SSaturatingSubtractor, 2}
    padder::Array{SAdder, 2}
    nadder::Array{SAdder, 2}
end
SSignedMatMultiplier(nrows, ncols) =
    SSignedMatMultiplier(
        SMatMultiplier(nrows, ncols),
        SMatMultiplier(nrows, ncols),
        SMatMultiplier(nrows, ncols),
        SMatMultiplier(nrows, ncols),
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SAdder() for i in 1:nrows, j in 1:ncols],
        [SAdder() for i in 1:nrows, j in 1:ncols]
    )
function (op::SSignedMatMultiplier)(x::VecOrMat{SBit}, y::VecOrMat{SBit})
    pp = op.ppmult(pos.(x), pos.(y))
    pn = op.pnmult(pos.(x), neg.(y))
    np = op.npmult(neg.(x), pos.(y))
    nn = op.nnmult(neg.(x), neg.(y))

    s₁ = map.(op.sub1, pp, pn)
    s₂ = map.(op.sub2, pn, pp)
    s₃ = map.(op.sub3, np, nn)
    s₄ = map.(op.sub4, nn, np)

    z = zip(map.(op.padder, s₁, s₄), map.(op.nadder, s₂, s₃)) |> collect

    return (size(z, 2) == 1) ? dropdims(z; dims = 2) : z
end

"""
    SL2Normer

A stochastic bitstream L2-norm operator.
"""
@kwdef struct SL2Normer <: SOperator
    dot::SSignedMatMultiplier = SSignedMatMultiplier(1, 1)
    buffer::CircularBuffer{Matrix{SBit}} = CircularBuffer{Matrix{SBit}}(1)
    root::SSquareRoot = SSquareRoot()
end
function (op::SL2Normer)(x::Vector{SBit})
    # get row vector
    dummybit = map(λ -> SBit((false, false), λ.value, λ.id), permutedims(x))
    xt = isempty(op.buffer) ? dummybit : pop!(op.buffer)
    push!(op.buffer, permutedims(x))

    # compute inner product
    xtx = op.dot(xt, x)[1]

    # compute root
    z = op.root(SBit(xtx, x[1].value, x[1].id))

    return z
end


text:
#=
    Install Package:
        In Julia REPL:
            1. type  ']'  to enter the  'pkg'  interface
            2. type  ' add "https://github.com/Mengzhen-Zhang/Symplectic.jl.git" '
=#
module Symplectic

using LinearAlgebra

const tolerance = 10^-6

const Graph = AbstractMatrix{Bool}
const Colors = Vector
const Vertex = Int
const Colored = Bool
const NofModes = Int

setprecision(BigFloat, 2048)

# Define DataType Symp for Symplectic Matrices
include("symp.jl")

# Functions for Generating Symplectic Matrices
include("symp_generators.jl")

# Functions for Generating Color Sets
include("color_sets.jl")

# Implemention of Interfernce Based Algorithms
include("interference_based_sequence.jl")

# Implemention of Teleportation Based Symplectic Matrix Generator
include("teleportation_based_generator.jl")

# Implemention of (Inverse) Symplectic Cayley Transform
include("cayley_transform.jl")

end


text:
export chain

"""
    @chain head tail

Calls `head` on `tail`

# Examples
```julia
@test vcat(1) == @chain 1 vcat
```
"""
chain(head, tail) = Expr(:call, tail, head)

"""
    @chain head tail::Expr

Reinterprets `\_` in `tail` as `head`. Note that dot vectorization is broken by
this macro. Instead, use [`unweave`](@ref), [`@map`](@ref), or [`@broadcast`](@ref)

# Examples
```julia
@test vcat(2, 1) == @chain 1 vcat(2, _)
```
"""
chain(head, tail::Expr) = Expr(:let, tail, Expr(:(=), :_, head))

chain(head, tail::AnnotatedLine) =
    AnnotatedLine(tail.line, chain(convert(Expr, head), tail.expr) )

"""
    @chain es...

`reduce` `@chain` over `es`

# Examples
```julia
@test ( @chain 1 vcat(_, 2) vcat(_, 3) ) ==
    @chain ( @chain 1 vcat(_, 2) ) vcat(_, 3)
```
"""
chain(es...) = foldl(chain, es)

"""
    @chain e

Separate single begin blocks out into lines and [`chain`](@ref) them;
return single non-blocks.

# Examples
```julia
@test 1 == @chain 1

chain_block = @chain begin
    1
    vcat(_, 2)
end

@test chain_block == @chain 1 vcat(_, 2)
```
"""
chain(e) =
    if MacroTools.isexpr(e, :block)
        convert(Expr, chain(annotate(e.args)...) )
    else
        e
    end

@nonstandard chain
export @chain


text:
# tests to debug addHybridizationUpdateSmart
# and deleteHybridizationUpdate
# Claudia March 2015

# starting topology: Case G
include("../examples/case_g_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output.csv")
d = readDataCF(df)

currloglik,currxmin = optBL!(currT,d)
newT = deepcopy(currT);
count = 0
N = 100
move = whichMove(currT)
move = :CHdir
move = :MVorigin
move = :MVtarget
move = :nni

flag = proposedTop!(move,newT,true,count,N)
printEdges(newT)
printNodes(newT)

newloglik, newxmin = optBL!(newT,d)
isValid(newT)
node=searchHybridNode(newT);
node[1].number
deleteHybridizationUpdate!(newT,node[1])
printEdges(newT)
printNodes(newT)
success=addHybridizationUpdate!(newT);
success[1]
printEdges(newT)
printNodes(newT)

# ---------- debug addHybridizationUpdateSmart

# starting topology: Case G
include("../examples/case_g_example.jl");
currT = deepcopy(net);
printEdges(currT)

# real network: Case H
df = readtable("CaseH_output.csv")
d = readDataCF(df)

currloglik,currxmin = optBL!(currT,d)
newT = deepcopy(currT);
count = 0
N = 100
move = whichMove(currT)
move = :CHdir
move = :MVorigin
move = :MVtarget
move = :nni

flag = proposedTop!(move,newT,true,count,N)
printEdges(newT)
printNodes(newT)

newloglik, newxmin = optBL!(newT,d)
isValid(newT)
node=searchHybridNode(newT);
node[1].number
deleteHybridizationUpdate!(newT,node[1])
printEdges(newT)
printNodes(newT)
success=addHybridizationUpdateSmart!(newT)
printEdges(newT)
printNodes(newT)

afterOptBL!(newT,d)
currT=deepcopy(newT);


text:
const coords = [
(1,1),
(1,2),(2,1),
(3,1),(2,2),(1,3),
(1,4),(2,3),(3,2),(4,1),
(5,1),(4,2),(3,3),(2,4),(1,5),
(1,6),(2,5),(3,4),(4,3),(5,2),(6,1),
(7,1),(6,2),(5,3),(4,4),(3,5),(2,6),(1,7),
(1,8),(2,7),(3,6),(4,5),(5,4),(6,3),(7,2),(8,1),
(8,2),(7,3),(6,4),(5,5),(4,6),(3,7),(2,8),
(3,8),(4,7),(5,6),(6,5),(7,4),(8,3),
(8,4),(7,5),(6,6),(5,7),(4,8),
(5,8),(6,7),(7,6),(8,5),
(8,6),(7,7),(6,8),
(7,8),(8,7),
(8,8),
]

struct DCACBlock
    DC::Int
    AC::Vector{Int}
end

function zigzag(x::Matrix)::Vector
    values = []
    for c in coords
        push!(values, x[c...])
    end
    return values
end

function unzigzag(x::Vector)::Matrix
    retval = zeros(8,8)
    for (i,v) in enumerate(x)
        retval[coords[i]...] = v
    end
    return retval
end

function value_to_dc_code(v::Int)::Tuple{Int,BitArray}
    if v == 0
        return (0,[])
    else
        numbits = Int(floor(log2(abs(v))))+1
        bits = []
        if v > 0
            bits = [parse(Int, i) for i in string(v, base=2)]
        else
            bits = [parse(Int, i) for i in string(abs(v), base=2)]
            bits = 1 .- bits
        end
        return numbits,bits
    end
end

function dc_code_to_value(numbits,bits::BitArray)::Int
    if numbits == 0
        return 0
    else
        if bits[1] == 1
            s = *([i ? "1" : "0" for i in bits]...)
            v = parse(Int,s,base=2)
            return v
        else
            s = *([i ? "0" : "1" for i in bits]...)
            v = parse(Int,s,base=2)
            return -1*v
        end
    end
end

function length_to_huffman_string(b::Int)::BitArray
    d::Dict{Int,BitArray} = Dict(
    0=>[0,0],
    1=>[0,1,0],
    2=>[0,1,1],
    3=>[1,0,0],
    4=>[1,0,1],
    5=>[1,1,0],
    6=>[1,1,1,0],
    7=>[1,1,1,1,0],
    8=>[1,1,1,1,1,0],
    9=>[1,1,1,1,1,1,0],
    10=>[1,1,1,1,1,1,1,0],
    11=>[1,1,1,1,1,1,1,1,0],
    )
    return d[b]
end

function huffman_string_to_length(b::BitArray)::Union{Int,Nothing}
    try
        d::Dict{BitArray,Int} = Dict(
        [0,0]=>0,
        [0,1,0]=>1,
        [0,1,1]=>2,
        [1,0,0]=>3,
        [1,0,1]=>4,
        [1,1,0]=>5,
        [1,1,1,0]=>6,
        [1,1,1,1,0]=>7,
        [1,1,1,1,1,0]=>8,
        [1,1,1,1,1,1,0]=>9,
        [1,1,1,1,1,1,1,0]=>10,
        [1,1,1,1,1,1,1,1,0]=>11,
        )
        return d[b]
    catch
        return nothing
    end
end

function value_to_dc_bits(v::Int)::BitArray
    l, bits = value_to_dc_code(v)
    ls = length_to_huffman_string(l)
    return vcat(ls,bits)
end

function bits_to_dc_value(b::BitArray)::Int
    l = nothing
    ls::BitArray = []
    while isnothing(l)
        push!(ls,popfirst!(b))
        l = huffman_string_to_length(ls)
    end
    dcbits = copy(b[1:l])
    for i=1:l
        popfirst!(b)
    end
    return dc_code_to_value(l,dcbits)
end

function ac_value_to_bits(v::Int)::BitArray
    return [1,1,0,0]
end

function bits_to_ac_value(b::BitArray)::Int
    return 0
end

function encode_zigzagged_block(x::Vector)::BitArray
    dc = x[1]
    dcarray = value_to_dc_bits(dc)
    ac = x[2:end]
    return vcat(dcarray, [1, 1, 0, 0])
end

function decode_zigzagged_block(x::BitArray)::Vector{Int}
    v = bits_to_dc_value(x)
    a = zeros(63)
    idx = 1
    @assert x[1:4] == [1,1,0,0]
    while x[1:4] != [1,1,0,0]
        #parse ac
        acval = bits_to_ac_value(x)
        a[idx] = acval
        idx += 1
    end
    for i=1:4
        popfirst!(x)
    end
    return vcat([v], a)
end


text:
# Define # of binding sites for clusters with n monomers up to imax
function Ns(imax, prefactor)

n  = range(1, stop=imax)
Ns = round.(Int64, prefactor.*n.^(2/3))

return Ns

end


text:
# Difference in NH and SH zenith angles at time x in given year
function zdiff(x, year)
    date = xtomarchdate(x,year)
    theta_s, dist = daily_zenith_angle(date, FT(-45), param_set)
    theta_n, dist = daily_zenith_angle(date, FT(45), param_set)
    return theta_n - theta_s
end

# x is date relative to March 1, with 1.00 representing March 1 00:00
function xtomarchdate(x, year)
    basedate = Dates.DateTime(year, 3, 1)
    deltat = Dates.Second(round((x-1)*Planet.day(param_set)))
    return basedate + deltat
end

days = zeros(length(1900:2100))
for (i,year) in enumerate(1900:2100)
    f = (x -> zdiff(x, year))
    days[i] = find_zeros(f,1.,30)[1]
end

# test mean is about March 21
@test mean(days) ≈ 21 atol=1

# test decreasing
@test mean(days[:100]) > mean(days[100:end])


text:
using JSON
Nproc = 10

rmprocs(workers())
addprocs(Nproc)

println("workers: ",length(workers()))

@everywhere using Base.Test
@everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase
@everywhere using Optim

# local
#@everywhere path="/home/emma/github/fit_behaviour_julia/"

#cluster
@everywhere path="/home/emma/fit_behaviour_julia/"

@everywhere path_functions=path*"functions/"
@everywhere path_scripts=path*"scripts_emma/"

@everywhere include(path_functions*"parameter_functions.jl")
@everywhere include(path_functions*"functions_julia.jl")
@everywhere include(path_functions*"utils.jl")

##############################
Nframes=10
Tframe=0.2
Nstim= 1000

model_num="8"

# Model only c2 and sigma
args=[ ("c2",[1]),("c4",0),("st",0),("x0",0),("bias_d",0),("hbias_d",0),("sigma_a",[1])]
consts=[("c2",[2]),("c4",[1]),("st",[1,2,3,4]),("x0",[1,2]),("bias_d",[1]),("hbias_d",[1]),("sigma_a",0)]

x = [0.6,0.4]
z = [0,1,2,2,0,0,0,0,0,0]
param=make_dict(args,x,consts,z)

ub=[10., 10.]
lb=[-10., -10]
reg=[1,1]

model_f=Dict(:c2_f=>c2_urgency,:c4_f=>c4_const,:c6_f=>c6_zeros,:bias_d_f=>bias_d_const,
:hbias_d_f=>history_bias_d_const,:st_f=>st_logistic_bias2,:x0_f=>initial_condition_bias_hbias,
:ro_f=>read_out_perfect,:lambda=>0.0)

data = Dict()
open(path*"processed_data.json", "r") do f
    global data
    dicttxt = readstring(f)  # file information to string
    data=JSON.parse(dicttxt)  # parse and transform data
end

# print both dictionaries
# println(dict2)
stimulus = data[1]             # Data frames obtained from a gaussian distribution of a determined coherence
coherences = data[2]           # Target coherences used for each trial
rewards = data[3]              # Correct choice side
decisions = data[4]            # Actual decision made by the rat each trial
performances = data[5]         # 0 if the rat chose the correct choice, 1 if the rat chose the incorrect choice
target_sigmas = data[6]        # Target sigmas used for each trial

#one session
stim = zeros(size(stimulus[1][1])[1],size(stimulus[1][1][1])[1])
for s in size(stimulus[1][1])
    for i in size(stimulus[1][1][1])
        stim[s,i] = stimulus[1][1][s][i]
    end
end
println(size(stim))
choices = decisions[1][1]

##############################
# Per 30000 trials fer synthetic data i calcular com de bé calcula els parametres
Nframes=size(stim[1])
Tframe=0.2
Nstim= size(stim)

history_bias=zeros(Nstim)

PR = ProbRight2(stim,history_bias,Tframe,param,model_f)
choices = []
for prob in PR
  if rand() < prob
    append!(choices,1)
  else
    append!(choices,-1)
  end
end

# Call function ComputeLL2 to create a MAtrix of LL parametrically chanching two parameters
# The maximum LL should be for the genereative parameters (the parameters that we use to create the data)

modl = Dict()
for i in 1:100
    println(i)
    x = rand(1:10,2)
    x = x/10
    modl_f = ModelFitting(stim,choices,history_bias, Tframe, args,consts,ub,lb, x,z,model_f)
    modl[string(i)] = modl_f
end

## SAVE THE RESULTS ##
stringdata = JSON.json(modl)

# write the file with the stringdata variable information
open("datafit_rat_2par.json", "w") do f
        write(f, stringdata)
     end
println("done")

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
@testset "Datasets" begin

	mnist = BitsMNIST.Datasets.mnist()
	@test length(mnist["train_x"]) == 60000
	@test length(mnist["train_x"][1]) == 784
	@test length(mnist["test_x"][1]) == 784
	@test length(mnist["train_y"]) == 60000

	noisy = BitsMNIST.Datasets.noisymnist()
	@test noisy["test_y"][1] == BitsMNIST.Datasets.NOISE_LABEL

end