text:
function createStartingPositionsmein2(x0,x1,y0,y1,N,nx,ny,alpha,verschx,verschy)
X = [];
Y = [];
ALPHA = [];
eps = 0.0000001;
unten = 0;
oben = 0;

xmid = (x0+x1)/2-verschx;
ymid = (y0+y1)/2+verschy;

    if (mod(N,2) == 0)
        unten = N/2;
        oben = (N/2)-1;
    else
        unten = (N-1)/2;
        oben = (N-1)/2;
    end

Deltas = -unten:1:oben;
nRaysPerAngle = length(Deltas);
omega = [cos(alpha),sin(alpha)];
omegaPerp = [-sin(alpha),cos(alpha)];
mycolor = rand(1,3);

    for j = 1:nRaysPerAngle

        delta = Deltas[j];

        x = xmid+delta*omegaPerp[1];
        y = ymid+delta*omegaPerp[2];
        vx = -cos(alpha);
        vy = -sin(alpha);
        I,J,L,xresult,yresult = distance2(x,y,nx,ny,x0,x1,y0,y1,vx,vy);

        X = [X;xresult];
        Y = [Y;yresult];
    end

return X,Y,Deltas;

end


text:
using GridGames
using MuZero
using StatsBase

cd(@__DIR__)
include("game.jl")

play(Tac)

config = getconfig(Tac)
play_net = new_test_network(config)
model_path = "~/models/MuZeroModels.jl/tictactoe/1a4cb71d-3dd6-4b6e-ad21-7cc9a5ce5457/model.bson"
load!(play_net, model_path)

struct TacMuZeroPolicy <: AbstractGridPolicy
    net::Network
end

function GridGames.getaction(p::TacMuZeroPolicy, s::AbstractGridState)
    @show s
    output = initial_inference(p.net, s) |> SingleNetworkOutput
    π = output |> getpolicy
    @show π
    argmax(π)
end

policy = TacMuZeroPolicy(
    play_net
)

play(Tac, policy)

ext: jl


text:
@testset "Wester Problems" begin
    (a, b, m, x, ) = @variables a b m x

    #= ::Package:: =#

    #= ::Title:: =#
    #=Michael*Wester=#

    #= Gradshteyn*and*Ryzhik*2.244(8) =#
    @test_int [(-5 + 3*x)^2/(-1 + 2*x)^(7/2), x, 2, -(49/(20*(-1 + 2*x)^(5/2))) + 7/(2*(-1 + 2*x)^(3/2)) - 9/(4*sqrt(-1 + 2*x))]

    #= => 1/(2*m*sqrt*(10))*log*((-5 + e^(m*x)*sqrt*(10))/(-5 - e^(m*x)*sqrt*(10)))
          (Gradshteyn*and*Ryzhik*2.314) =#
    @test_int [1/(-5/ℯ^(m*x) + 2*ℯ^(m*x)), x, 2, -(atanh(sqrt(2/5)*ℯ^(m*x))/(sqrt(10)*m))]

    #= This*example*involves*several*symbolic*parameters
       => 1/sqrt(b^2 - a^2)*log*((sqrt*(b^2 - a^2)*tan*(x/2) + a + b)/
                                (sqrt*(b^2 - a^2)*tan*(x/2) - a - b))   (a^2 < b^2)
          (Gradshteyn*and*Ryzhik*2.553(3)) =#
    #=
    [1/(a + b*cos(x)), x, 0, Assumptions -> a^2 < b^2,
     1/sqrt(b^2 - a^2)*log((sqrt(b^2 - a^2)*tan(x/2) + a + b)/
                           (sqrt(b^2 - a^2)*tan(x/2) - a - b))]
    =#
    @test_int [1/(a + b*cos(x)), x, 2, (2*atan((sqrt(a - b)*tan(x/2))/sqrt(a + b)))/(sqrt(a - b)*sqrt(a + b))]
    #= The*integral*of*1/(a + 3*cos*x + 4*sin*x)*can*have*4*different*forms
       depending*on*the*value*of*a !   (Gradshteyn*and*Ryzhik*2.558(4)) =#
    @test_int [1/(3 + 3*cos(x) + 4*sin(x)), x, 2, (1/4)*log(3 + 4*tan(x/2))]
    @test_int [1/(4 + 3*cos(x) + 4*sin(x)), x, 2, (-(1/3))*log(4 + 3*cot(pi/4 + x/2))]
    @test_int [1/(5 + 3*cos(x) + 4*sin(x)), x, 1, -1/(2 + tan(x/2)), -((4 - 5*sin(x))/(4*(4*cos(x) - 3*sin(x))))]
    #= => (a = 6)*2/sqrt(11)*arctan*((3*tan*(x/2) + 4)/sqrt(11)) =#
    @test_int [1/(6 + 3*cos(x) + 4*sin(x)), x, 3, x/sqrt(11) + (2*atan((4*cos(x) - 3*sin(x))/(6 + sqrt(11) + 3*cos(x) + 4*sin(x))))/sqrt(11)]

    #= => x*log|x^2 - a^2| - 2*x + a*log|(x + a)/(x - a)|
          (Gradshteyn*and*Ryzhik*2.736(1)) =#
    #= [log(Abs(x^2 - a^2)), x, 0, x*log(Abs(x^2 - a^2)) - 2*x + a*log((x + a)/(x - a))] =#
    @test_int [(1/2)*log((-a^2 + x^2)^2), x, 4, -2*x + 2*a*atanh(x/a) + (1/2)*x*log((-a^2 + x^2)^2)]
    end

max_stars_repo_path: test/0 Independent test suites/Wester Problems.jl


text:
module Data

export CMUDict, cmudict

deps(path...) = joinpath(@__DIR__, "..", "..", "deps", path...)

function __init__()
  mkpath(deps())
end

include("cmudict.jl")
using .CMUDict

end


text:
# Sandpile64.jl
module Sandpile64

export sandpile_init64, run_sandpile64!, avalanche64!, add_grain!, is_unstable64

function sandpile_init64(x::Int64, y::Int64, setup::String="random")
    """
            sandpile_init(X, Y, setup)
    Initialize a 2d sandpile grid with size (`X`, `Y`). Can either have a zero or random
    initial condition.
    Arguments:
            x (Int64): length of grid along x
            y (Int64): length of grid along y
            setup (String): Type of inital setup for grid. Can either be `zero` or `random`.
    Returns:
            z (Matrix{Int64}): sandpile grid z with dimensions (x,y).
    """
    if setup == "random"
        return rand(0:3, x, y)
    elseif setup == "zero"
        return zeros(Int64, x, y)
    else
        error("`setup` must have value of `random` or `zero`")
    end
end

function add_grain64!(z::Array{Int64,2})
    """
            add_grain!(z)
    Adds a grain to a random location on the  2d array grid for the Abelian sandpile model, `z`.
    Arguments:
            z (Array{Int64, 2}): 2D sandpile grid.
    """
    return z[rand(1:size(z)[1]), rand(1:size(z)[2])] += 1
end

function is_unstable64(z::Array{Int64,2}, fᵪ::Int64=4)
    """
            is_unstable(z, fᵪ)
    Checks whether or not any location on the grid is unstable (has value above critical value).
    Arguments:
            z (Array{Int64, 2}): 2D sandpile grid.
            fₓ (Int64): Critical value for sandpile model (Must be > 1)
    Returns:
            (bool): Whether grid is stable or unstable.
    """
    if fᵪ <= 1
        error("`fᵪ` must be int greater than  1")
    end
    if length(findall(>=(fᵪ), z)) > 0
        return true
    else
        return false
    end
end

function avalanche64!(z::Array{Int64,2}, fᵪ=4)
    """
            avalanche(z, fᵪ)
    Update grid to follow rules on locations that are unstable.
    Follow Eq 3.2 of Self-organized criticality by Bak, Tang (1988).
    Returns how many avalanches (or slides s) occur for input state after collapsing all
    unstable locations.
    Arguments:
            z (Array{Int64, 2}): 2D sandpile grid.
            fₓ (Int64): Critical value for sandpile model (Must be > 1)
    Returns:
            (Int64): Number of unstable locations collapsed.
    """
    s::Int64 = 0
    unstable_locs = findall(>=(fᵪ), z)
    for crit_loc in unstable_locs
        z[crit_loc] -= 4
        s += 1
        neighbors = [
            CartesianIndex(0, 1),
            CartesianIndex(0, -1),
            CartesianIndex(1, 0),
            CartesianIndex(-1, 0),
        ]
        for neighbor in neighbors
            try
                z[crit_loc + neighbor] += 1
            catch
            end
        end
    end
    return s
end

function run_sandpile64!(z::Array{Int64,2}, N::Int64; N_crit::Bool=false, fᵪ::Int64=4)
    """
            run_sandpile(z, N, fᵪ)
    Run a sandpile model on an input grid `z` by adding `N` grains.
    Arguments:
            z (Array{Int64, 2}): 2D sandpile grid.
            N (Int64): Number of grains to add
            N_crit (boolean): Do you want to to all added `N` grains to result in collapses?
                              If set to `false`, `N` grains will be added to system but
                              returned sizes will only be those that result in avalanches.
                              Warning:  longer runtimes if set to `true`.
            fₓ (Int64): Critical value for sandpile model (Must be > 1)
    Returns
            (Vector{Int64}): Sizes of avalanches for each grain added.
    """
    s_list = zeros(Int64, N)
    if N_crit
            i = 1
            while i <= N
                add_grain64!(z)
                s = 0
                while is_unstable64(z, fᵪ)
                    s += avalanche64!(z)
                end
                if s > 0
                    s_list[i] = s
                    i += 1
                end
            end
    else
        for i in 1:N
            add_grain64!(z)
            s = 0
            while is_unstable64(z, fᵪ)
                s += avalanche64!(z)
            end
            s_list[i] = s
        end
    end
    return s_list[s_list .> 0]
end
end

max_stars_repo_name: romeroraa/sandpile-julia

max_issues_repo_name: romeroraa/sandpile-julia

max_forks_repo_name: romeroraa/sandpile-julia


text:
module Wflow

using Dates
using TOML
using Graphs
using NCDatasets
using StaticArrays
using Statistics
using UnPack
using Random
using BasicModelInterface
using FieldMetadata
using Parameters
using DelimitedFiles
using ProgressLogging
using CFTime
using Base.Threads
using Glob
using NCDatasets: MFDataset

@metadata get_units "mm Δt-1"

const BMI = BasicModelInterface
const Float = Float64
const CFDataset = Union{NCDataset,MFDataset}

mutable struct Clock{T}
    time::T
    iteration::Int
    Δt::Second
end

function Clock(config)
    # this constructor is used by reset_clock!, since if the Clock has already
    # been constructed before, the config is complete
    calendar = get(config, "calendar", "standard")::String
    starttime = cftime(config.starttime, calendar)
    Δt = Second(config.timestepsecs)
    Clock(starttime, 1, Δt)
end

function Clock(config, reader)
    nctimes = reader.dataset["time"][:]
    # if the config file does not have a start or endtime, folow the NetCDF times
    # and add them to the config
    # if the timestep is not given, use the difference between NetCDF time 1 and 2
    starttime = get(config, "starttime", nothing)
    if starttime === nothing
        starttime = first(nctimes)
        config.starttime = starttime
    end
    endtime = get(config, "endtime", nothing)
    if endtime === nothing
        endtime = last(nctimes)
        config.endtime = endtime
    end
    timestepsecs = get(config, "timestepsecs", nothing)
    if timestepsecs === nothing
        timestepsecs = Dates.value(Second(nctimes[2] - nctimes[1]))
        config.timestepsecs = timestepsecs
    end

    calendar = get(config, "calendar", "standard")::String
    starttime = cftime(config.starttime, calendar)
    Δt = Second(timestepsecs)

    fews_run = get(config, "fews_run", false)::Bool
    if fews_run
        starttime = starttime + Δt
    end

    Clock(starttime, 1, Δt)
end

include("io.jl")

"""
    Model{N,L,V,R,W}

Composite type that represents all different aspects of a Wflow Model, such as the
network, parameters, clock, configuration and input and output.
"""
struct Model{N,L,V,R,W}
    config::Config  # all configuration options
    network::N  # connectivity information, directed graph
    lateral::L  # lateral model that holds lateral state, moves along network
    vertical::V  # vertical model that holds vertical state, independent of each other
    clock::Clock  # to keep track of simulation time
    reader::R  # provides the model with dynamic input
    writer::W  # writes model output
end

# prevent a large printout of model components and arrays
Base.show(io::IO, m::Model) = print(io, "model of type ", typeof(m))

include("horizontal_process.jl")
include("hbv.jl")
include("sbm.jl")
include("sediment.jl")
include("reservoir_lake.jl")
include("hbv_model.jl")
include("sbm_model.jl")
include("sediment_model.jl")
include("flow.jl")
include("vertical_process.jl")
include("groundwater/connectivity.jl")
include("groundwater/aquifer.jl")
include("groundwater/boundary_conditions.jl")
include("sbm_gwf_model.jl")
include("utils.jl")
include("bmi.jl")
include("subdomains.jl")

"""
    run(tomlpath::String)
    run(config::Config)
    run(model::Model)
    run()

Run an entire simulation starting either from a path to a TOML settings file,
a prepared `Config` object, or an initialized `Model` object. This allows more flexibility
if you want to for example modify a `Config` before initializing the `Model`.

The 0 argument version expects ARGS to contain a single entry, pointing to the TOML path.
This makes it easier to start a run from the command line without having to escape quotes:

    julia -e "using Wflow; Wflow.run()" "path/to/config.toml"
"""
function run(tomlpath)
    config = Config(tomlpath)
    run(config)
end

function run(config::Config)
    modeltype = config.model.type

    model = if modeltype == "sbm"
        initialize_sbm_model(config)
    elseif modeltype == "sbm_gwf"
        initialize_sbm_gwf_model(config)
    elseif modeltype == "hbv"
        initialize_hbv_model(config)
    elseif modeltype == "sediment"
        initialize_sediment_model(config)
    else
        error("unknown model type")
    end
    load_fixed_forcing(model)
    run(model)
end

function run(model::Model; close_files = true)
    @unpack network, config, writer, clock = model

    # in the case of sbm_gwf it's currently a bit hard to use dispatch
    model_type = config.model.type::String
    update_func = model_type == "sbm_gwf" ? update_sbm_gwf : update

    # determine timesteps to run
    calendar = get(config, "calendar", "standard")::String
    starttime = clock.time
    Δt = clock.Δt
    endtime = cftime(config.endtime, calendar)
    times = range(starttime, endtime, step = Δt)

    @info "Run information" model_type starttime Δt endtime nthreads()
    @progress for (i, time) in enumerate(times)
        @debug "Starting timestep" time timestep = i
        load_dynamic_input!(model)
        model = update_func(model)
    end

    # write output state NetCDF
    # undo the clock advance at the end of the last iteration, since there won't
    # be a next step, and then the output state falls on the correct time
    rewind!(clock)
    write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)

    reset_clock!(model.clock, config)

    # option to support running function twice without re-initializing
    # and thus opening the NetCDF files
    if close_files
        Wflow.close_files(model, delete_output = false)
    end
    return model
end

function run()
    usage = "Usage: julia -e 'using Wflow; Wflow.run()' 'path/to/config.toml'"
    n = length(ARGS)
    if n != 1
        throw(ArgumentError(usage))
    end
    toml_path = only(ARGS)
    if !isfile(toml_path)
        throw(ArgumentError("File not found: $(toml_path)\n" * usage))
    end
    Wflow.run(toml_path)
end

end # module


text:
function QueryOperators.query(source::Observable)
    return source
end

function QueryOperators.map(source::Observable, f::Function, f_expr::Expr)
    return MapObservable(f, source)
end

function DataFrames.DataFrame(source::Observable)
    df = DataFrame()
    subscribe(source, df)
    return df
end

function Base.Array(source::Observable)
    ar = Array{Any}(0)
    subscribe(source, ar)
    return ar
end

function onnext(obs::Array, value)
    push!(obs, value)
end

function oncompleted(obs::Array)
end


text:
module SwagUI

using JSON
using Genie, Genie.Router

include("Options.jl")
include("Builders.jl")

const JS_PATH = joinpath(dirname(dirname(@__FILE__)), "dist")
export Options, render_swagger, serve_assets, build_js_string, build_html_string

"""
    render_swagger(swagger_doc::Union{Dict{String, Any}, Nothing})::String

Render Swagger UI's HTML as `String`, based on the options passed in. If a vaid `swagger_doc` is passed,
the `url` and `urls` in the `Options.swagger_options` are ignored by the UI. If `swagger_doc` is `nothing`,
`url` or `urls` has to be included in the `Options.swagger_options`.

# Arguments
Required:
- `swagger_doc::Union{Dict{String, Any}, Nothing}` : The swagger specification file `swagger.json` as a `Dict{String, Any}`.
"""
function render_swagger(swagger_doc::Union{Dict{String, Any}, Nothing}; options::Options=Options())::String
    serve_assets(JS_PATH, excludes=["index.html"])

    js_string = build_js_string(options, swagger_doc)
    html_string = build_html_string(options, js_string)

    return html_string
end

"""
    serve_assets(path::String)

Serve static files from a given folder `path`. `excludes` contains a list of file names that need to be skipped.

# Arguments
Required:
- `path::String` : The path of the folder that contains the files to serve.
"""
function serve_assets(path::String; excludes::Array{String, 1}=Array{String, 1}())
    for (root, _, files) in walkdir(path)
        for file in files
            if !(file in excludes)
                route(file) do
                    serve_static_file(file, root=root)
                end
            end
        end
    end
end

end # module


text:
@testset "group" begin
    @test group(identity, 11:20)::HashDictionary == HashDictionary(11:20, (x->[x]).(11:20))
    @test group(iseven, 1:10)::HashDictionary == dictionary(true => [2,4,6,8,10], false => [1,3,5,7,9])

    @test group(iseven, x -> x*2, 1:10)::HashDictionary == dictionary(true => [4,8,12,16,20], false => [2,6,10,14,18])

    @test group((x,y) -> iseven(x+y), (x,y) -> x, 1:10, [1,3,4,2,5,6,4,2,3,9])::HashDictionary == dictionary(true => [1,4,5,6,8,9], false => [2,3,7,10])
end

@testset "grouponly" begin
    @test grouponly(identity, 11:20)::HashDictionary == HashDictionary(11:20, 11:20)
    @test_throws IndexError grouponly(iseven, 11:20)
end

@testset "groupunique" begin
    @test groupunique(identity, 11:20)::HashDictionary{Int, HashIndices{Int}} == HashDictionary(11:20, map(x -> HashIndices([x]), 11:20))
    @test groupunique(iseven, 11:20)::HashDictionary{Bool, HashIndices{Int}} == dictionary(true => HashIndices([12,14,16,18,20]), false => HashIndices([11,13,15,17,19]))
end

@testset "groupfind" begin
    @test groupfind(identity, 11:20) == Dictionary(11:20, (x->[x]).(1:10))
    @test groupfind(iseven, 11:20) == dictionary(true => [2,4,6,8,10], false => [1,3,5,7,9])
end

@testset "groupview" begin
    @test groupview(identity, 11:20)::GroupDictionary == group(identity, 11:20)::HashDictionary
    @test groupview(iseven, 11:20)::GroupDictionary == group(iseven, 11:20)::HashDictionary
end

@testset "groupreduce" begin
    @test groupreduce(identity, +, 1:10) == HashDictionary(1:10, 1:10)
    @test groupreduce(iseven, +, 1:10) == dictionary(true => 30, false => 25)

    @test groupreduce(iseven, x -> x*2, +, 1:10) == dictionary(true => 60, false => 50)

    @test groupreduce(iseven, x -> x*2, +, 1:10; init=10) == dictionary(true => 70, false => 60)

    @test groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, 1:10; init=10) == dictionary(true => 120)
    @test groupreduce((x,y) -> iseven(x+y), (x,y) -> x+y, +, 1:10, [1,3,4,2,5,6,4,2,3,9]; init=10) == dictionary(true => 62, false => 52)

    @test groupcount(iseven, 1:10) == dictionary(true => 5, false => 5)
    @test groupsum(iseven, 1:10) == dictionary(true => 30, false => 25)
    @test groupprod(iseven, 1:10) == dictionary(true => 2*4*6*8*10, false => 1*3*5*7*9)
    @test groupfirst(iseven, 1:10) == dictionary(true => 2, false => 1)
    @test grouplast(iseven, 1:10) == dictionary(true => 10, false => 9)
end


text:
struct BlockArrayCoo{T,N,A,X} <: AbstractBlockArray{T,N}
  blocks::Vector{A}
  blockids::Vector{NTuple{N,Int}}
  axes::X
  ptrs::Array{Int,N}
  zero_blocks::Vector{A}
  function BlockArrayCoo(
    blocks::Vector{A},
    blockids::Vector{NTuple{N,Int}},
    axes::NTuple{N},
    ptrs::Array{Int,N}=_compute_ptrs(blockids,axes),
    zero_blocks::Vector{A}=_compute_zero_blocks(A,ptrs,axes)) where {T,N,A<:AbstractArray{T,N}}

    X = typeof(axes)
    new{T,N,A,X}(blocks,blockids,axes,ptrs,zero_blocks)
  end
end

function BlockArrayCoo(
  blocks::Vector{A},ax::NTuple{N},ptrs::Array{Int,N},zero_blocks::Vector{A}) where {A,N}
  cis = CartesianIndices(ptrs)
  blockids = NTuple{N,Int}[]
  for ci in cis
    p = ptrs[ci]
    if p>0
      push!(blockids,Tuple(ci))
    end
  end
  BlockArrayCoo(blocks,blockids,ax,ptrs,zero_blocks)
end

function BlockArrayCoo(allblocks::AbstractArray{A,N},ax::NTuple{N},mask::AbstractArray{Bool,N}=fill(true,size(allblocks))) where {A,N}
  blocks = A[]
  zero_blocks = A[]
  ptrs = zeros(Int,size(allblocks))
  pp = 1
  pn = 1
  cis = CartesianIndices(allblocks)
  for ci in cis
    if mask[ci]
      push!(blocks,allblocks[ci])
      ptrs[ci] = pp
      pp += 1
    else
      push!(zero_blocks,allblocks[ci])
      ptrs[ci] = -pn
      pn += 1
    end
  end
  BlockArrayCoo(blocks,ax,ptrs,zero_blocks)
end

#  axs = _compute_axes_from_allblocks(allblocks)
#
#function _compute_axes_from_allblocks(allblocks)
#  N = ndims(allblocks)
#  s = size(allblocks)
#  axtmp = [ zeros(Int,s[n]) for n in 1:N]
#  cis = CartesianIndices(allblocks)
#  for ci in cis
#    block = allblocks[ci]
#    sb = size(block)
#    for (n,cin) in enumerate(Tuple(ci))
#      axtmp[n][cin] = sb[n]
#    end
#  end
#  axs = map(blockedrange,Tuple(axtmp))
#  axs
#end

const BlockMatrixCoo = BlockArrayCoo{T,2} where T
const BlockVectorCoo = BlockArrayCoo{T,1} where T

zeros_like(a::AbstractArray) = zeros_like(a,axes(a))

function zeros_like(a::BlockArrayCoo,axs::Tuple)
  blocks = eltype(a.blocks)[]
  blockids = eltype(a.blockids)[]
  BlockArrayCoo(blocks,blockids,axs)
end

function zeros_like(a::AbstractArray,axs::Tuple)
  T = eltype(a)
  b = similar(a,T,axs)
  fill!(b,zero(T))
  b
end

function zeros_like(::Type{<:BlockArrayCoo{T,N,A}},axs::Tuple) where {T,N,A}
  blocks = A[]
  blockids = NTuple{N,Int}[]
  BlockArrayCoo(blocks,blockids,axs)
end

function zeros_like(::Type{A},axs::Tuple) where A <: AbstractArray
  a = similar(A,axs)
  fill!(a,zero(eltype(a)))
  a
end

Base.similar(a::BlockArrayCoo) = similar(a,eltype(a),axes(a))

Base.similar(a::BlockArrayCoo,axs::Tuple) = similar(a,eltype(a),axs)

function Base.similar(a::BlockArrayCoo{S,N} where S,::Type{T}, axs::NTuple{N,<:BlockedUnitRange}) where {T,N}
  _similar_block_array_coo(a,T,axs)
end

function _similar_block_array_coo(a,::Type{T},axs) where T
  if length(a.blocks) != 0
    A = typeof(similar(first(a.blocks),T))
  else
    A = typeof(similar(first(a.zero_blocks),T))
  end
  blocks = A[]
  for p in 1:length(a.blocks)
    I = a.blockids[p]
    laxs = map( local_range, axs, I)
    block = similar(a.blocks[p],T,laxs)
    push!(blocks,block)
  end
  BlockArrayCoo(blocks,a.blockids,axs,a.ptrs)
end

function local_range(a::BlockedUnitRange,k::Integer)
  n = length(a[Block(k)])
  Base.OneTo(n)
end

function _compute_ptrs(blockids,axes)
  s = map(i->first(blocksize(i)),axes)
  ptrs = zeros(Int,s)
  for (i,c) in enumerate(blockids)
    ptrs[c...] = i
  end
  m = 1
  for (k,p) in enumerate(ptrs)
    if p ==0
      ptrs[k] = -m
      m += 1
    end
  end
  ptrs
end

function _compute_zero_blocks(::Type{A},ptrs,axs) where A
  cis = CartesianIndices(ptrs)
  zero_blocks = A[]
  for ci in cis
    p = ptrs[ci]
    if p<0
      I = Tuple(ci)
      laxs = map( local_range, axs, I)
      block = zeros_like(A,laxs)
      push!(zero_blocks,block)
    end
  end
  zero_blocks
end

function BlockArrays.getblock(a::BlockArrayCoo,i::Integer...)
  p = a.ptrs[i...]
  if p>0
    a.blocks[p]
  else
    a.zero_blocks[-p]
  end
end

@inline is_nonzero_block(a,b...) = ! is_zero_block(a,b...)

function is_zero_block(a,b::Block)
  i = convert(Tuple,b)
  is_zero_block(a,i...)
end

function is_zero_block(a,b::Block...)
  i = map(Int,b)
  is_zero_block(a,i...)
end

function is_zero_block(a,b::CartesianIndex)
  i = Tuple(b)
  is_zero_block(a,i...)
end

function is_zero_block(a::BlockArrayCoo,i::Integer...)
  p = a.ptrs[i...]
  @assert p != 0
  p < 0
end

function BlockArrays.getblock!(block,a::BlockArrayCoo,i::Integer...)
  p = a.ptrs[i...]
  if p>0
    block .= a.blocks[p]
  else
    block .= a.zero_blocks[-p]
  end
end

function BlockArrays.eachblock(a::BlockArrayCoo)
  cis = CartesianIndices(blocksize(a))
  blocks = map(ci->Block(Tuple(ci)),cis)
  ( a[block] for block in blocks  )
end

function enumerateblocks(a)
  blocks = eachblockindex(a)
  zip(blocks,eachblock(a))
end

function eachblockindex(a)
  cis = CartesianIndices(blocksize(a))
  map(ci->Block(Tuple(ci)),cis)
end

Base.size(a::BlockArrayCoo) = map(length,Base.axes(a))

Base.axes(a::BlockArrayCoo) = a.axes

function Base.getindex(a::BlockArrayCoo{T,N} where T,i::Vararg{Integer,N}) where N
  s = map(findblockindex,a.axes,i)
  ai = a[s...]
  ai
end

function Base.getindex(a::BlockVectorCoo,i::Integer,j::Integer...)
  @assert all( j .== 1)
  s = map(findblockindex,a.axes,(i,))
  ai = a[s...]
  ai
end

function Base.setindex!(a::BlockArrayCoo{T,N} where T,v,i::Vararg{Integer,N}) where N
  s = map(findblockindex,a.axes,i)
  s = map(findblockindex,a.axes,(i,))
  I = Block(map(i->i.I[1],s)...)
  α = CartesianIndex(map(BlockArrays.blockindex,s))
  a[I][α] = v
end

function Base.setindex!(a::BlockVectorCoo,v,i::Integer,j::Integer...)
  @assert all( j .== 1)
  s = map(findblockindex,a.axes,(i,))
  I = Block(map(i->i.I[1],s)...)
  α = CartesianIndex(map(BlockArrays.blockindex,s))
  a[I][α] = v
end

function Base.:+(a::BlockArrayCoo{Ta,N},b::BlockArrayCoo{Tb,N}) where {Ta,Tb,N}
  @assert axes(a) == axes(b)
  @assert blockaxes(a) == blockaxes(b)
  I1 = first(eachblockindex(a))
  A = typeof(a[I1]+b[I1])
  blocks = A[]
  blockids = NTuple{N,Int}[]
  for (I,aI) in enumerateblocks(a)
    bI = b[I]
    if is_nonzero_block(a,I) || is_nonzero_block(b,I)
      block = aI + bI
      push!(blocks,block)
      push!(blockids,I.n)
    end
  end
  BlockArrayCoo(blocks,blockids,a.axes)
end

function Base.:-(a::BlockArrayCoo{Ta,N},b::BlockArrayCoo{Tb,N}) where {Ta,Tb,N}
  @assert axes(a) == axes(b)
  @assert blockaxes(a) == blockaxes(b)
  I1 = first(eachblockindex(a))
  A = typeof(a[I1]-b[I1])
  blocks = A[]
  blockids = NTuple{N,Int}[]
  for (I,aI) in enumerateblocks(a)
    bI = b[I]
    if is_nonzero_block(a,I) || is_nonzero_block(b,I)
      block = aI - bI
      push!(blocks,block)
      push!(blockids,I.n)
    end
  end
  BlockArrayCoo(blocks,blockids,a.axes)
end

function Base.:*(a::BlockArrayCoo,b::Number)
  *(b,a)
end

function Base.:*(a::Number,b::BlockArrayCoo)
  f(block)= a*block
  blocks = f.(b.blocks)
  BlockArrayCoo(blocks,b.blockids,b.axes,b.ptrs,b.zero_blocks)
end

function Base.:*(a::BlockMatrixCoo,b::BlockVectorCoo)
  @assert blocksize(a,2) == blocksize(b,1)
  A = typeof(a[Block(1,1)]*b[Block(1)])
  blocks = A[]
  blockids = Tuple{Int}[]
  for i in 1:blocksize(a,1)
    block = zeros_like(a[Block(i,1)]*b[Block(1)])
    for j in 1:blocksize(a,2)
      if is_nonzero_block(a,Block(i,j)) && is_nonzero_block(b,Block(j))
        block += a[Block(i,j)]*b[Block(j)]
      end
    end
    push!(blocks,block)
    push!(blockids,(i,))
  end
  axs = (axes(a)[1],)
  BlockArrayCoo(blocks,blockids,axs)
end

function Base.:*(a::BlockMatrixCoo,b::BlockMatrixCoo)
  @assert blocksize(a,2) == blocksize(b,1)
  A = typeof(a[Block(1,1)]*b[Block(1,1)])
  blocks = A[]
  blockids = Tuple{Int,Int}[]
  for i in 1:blocksize(a,1)
    for j in 1:blocksize(b,2)
      block = zeros_like(a[Block(i,1)]*b[Block(1,j)])
      for k in 1:blocksize(a,2)
        if is_nonzero_block(a,Block(i,k)) && is_nonzero_block(b,Block(k,j))
          block += a[Block(i,k)]*b[Block(k,j)]
        end
      end
      push!(blocks,block)
      push!(blockids,(i,j))
    end
  end
  axs = (axes(a)[1],axes(b)[2])
  BlockArrayCoo(blocks,blockids,axs)
end

function LinearAlgebra.mul!(c::BlockVectorCoo,a::BlockMatrixCoo,b::BlockVectorCoo)
  fill!(c,zero(eltype(c)))
  mul!(c,a,b,1,0)
end

function LinearAlgebra.mul!(c::BlockMatrixCoo,a::BlockMatrixCoo,b::BlockMatrixCoo)
  fill!(c,zero(eltype(c)))
  mul!(c,a,b,1,0)
end

function scale_entries!(c::BlockArrayCoo,β)
  for block in c.blocks
    scale_entries!(block,β)
  end
  c
end

function LinearAlgebra.mul!(c::BlockVectorCoo,a::BlockMatrixCoo,b::BlockVectorCoo,α::Number,β::Number)
  for I in 1:blocksize(a,1)
    cI = c[Block(I)]
    if is_nonzero_block(c,Block(I))
      scale_entries!(cI,β)
    end
    for J in 1:blocksize(a,2)
      if is_nonzero_block(a,Block(I,J)) && is_nonzero_block(b,Block(J))
        @assert is_nonzero_block(c,Block(I))
        aIJ = a[Block(I,J)]
        bJ = b[Block(J)]
        mul!(cI,aIJ,bJ,α,1)
      end
    end
  end
  c
end

function LinearAlgebra.mul!(c::BlockMatrixCoo,a::BlockMatrixCoo,b::BlockMatrixCoo,α::Number,β::Number)
  for I in 1:blocksize(a,1)
    for J in 1:blocksize(b,2)
      cIJ = c[Block(I,J)]
      if is_nonzero_block(c,Block(I,J))
        scale_entries!(cIJ,β)
      end
      for K in 1:blocksize(a,2)
        if is_nonzero_block(a,Block(I,K)) && is_nonzero_block(b,Block(K,J))
          @assert is_nonzero_block(c,Block(I,J))
          cIJ = c[Block(I,J)]
          aIK = a[Block(I,K)]
          bKJ = b[Block(K,J)]
          mymul!(cIJ,aIK,bKJ,α) # Hack to avoid allocations
          #mul!(cIJ,aIK,bKJ,α,1) # Why this leads to memory allocations?? Specially with Ints
        end
      end
    end
  end
  c
end

# Hack to avoid allocations
@inline function mymul!(cIJ,aIK,bKJ,α)
  @boundscheck begin
    @assert size(cIJ,1) == size(aIK,1)
    @assert size(cIJ,2) == size(bKJ,2)
    @assert size(aIK,2) == size(bKJ,1)
  end
  for i in 1:size(cIJ,1)
    for j in 1:size(bKJ,2)
      for k in 1:size(bKJ,1)
        @inbounds cIJ[i,j] += α*aIK[i,k]*bKJ[k,j]
      end
    end
  end
end

# Hack to avoid allocations
@inline function mymul!(cIJ::BlockMatrixCoo,aIK::BlockMatrixCoo,bKJ::BlockMatrixCoo,α)
  mul!(cIJ,aIK,bKJ,α,1)
end

function Base.fill!(a::BlockArrayCoo,v)
  @notimplementedif v != zero(v)
  for b in a.blocks
    fill!(b,v)
  end
  a
end

function LinearAlgebra.transpose(a::BlockMatrixCoo)
  blocks = transpose.(a.blocks)
  zero_blocks = transpose.(a.zero_blocks)
  blockids = [ (j,i) for (i,j) in a.blockids ]
  ax,ay = axes(a)
  axs = (ay,ax)
  ptrs = collect(transpose(a.ptrs))
  BlockArrayCoo(blocks,blockids,axs,ptrs,zero_blocks)
end

function Base.copy(a::BlockArrayCoo)
  blocks = copy.(a.blocks)
  blockids = copy(a.blockids)
  axes = copy.(a.axes)
  ptrs = copy(a.ptrs)
  zero_blocks = copy.(a.zero_blocks)
  BlockArrayCoo(blocks,blockids,axes,ptrs,zero_blocks)
end

function Base.copy!(a::BlockArrayCoo,b::BlockArrayCoo)
  copyto!(a,b)
end

function Base.copyto!(a::BlockArrayCoo,b::BlockArrayCoo)
  if a.ptrs === b.ptrs || a.ptrs == b.ptrs
    for p in 1:length(a.blocks)
      copyto!(a.blocks[p],b.blocks[p])
    end
  else
    for I in eachblockindex(a)
      if is_nonzero_block(a,I)
        copyto!(a[I],b[I])
      end
    end
  end
 a
end

# For blocks of blocks

struct TwoLevelBlockedUnitRange{CS} <: AbstractUnitRange{Int}
  global_range::BlockedUnitRange{CS}
  local_ranges::Vector{BlockedUnitRange{CS}}
end

function TwoLevelBlockedUnitRange(local_ranges::Vector{<:BlockedUnitRange})
  global_range = blockedrange(length.(local_ranges))
  TwoLevelBlockedUnitRange(global_range,local_ranges)
end

function BlockArrays.blockedrange(local_ranges::Vector{<:BlockedUnitRange})
  TwoLevelBlockedUnitRange(local_ranges)
end

function Base.similar(a::BlockArrayCoo{S,N} where S,::Type{T}, axs::NTuple{N,<:TwoLevelBlockedUnitRange}) where {T,N}
  _similar_block_array_coo(a,T,axs)
end

function Base.similar(::Type{T},axs::NTuple{N,<:TwoLevelBlockedUnitRange}) where {N,T<:Array}
  similar(T,map(length,axs))
end

local_range(a::TwoLevelBlockedUnitRange,k::Integer) = a.local_ranges[k]

Base.first(a::TwoLevelBlockedUnitRange) = first(a.global_range)
Base.last(a::TwoLevelBlockedUnitRange) = last(a.global_range)
BlockArrays.blockaxes(a::TwoLevelBlockedUnitRange) = blockaxes(a.global_range)
BlockArrays.blocklasts(a::TwoLevelBlockedUnitRange) = blocklasts(a.global_range)
BlockArrays.findblock(a::TwoLevelBlockedUnitRange,k::Integer) = findblock(a.global_range,k)
Base.getindex(a::TwoLevelBlockedUnitRange,i::Integer) = a.global_range[i]
Base.getindex(a::TwoLevelBlockedUnitRange,i::Block{1}) = a.global_range[i]
Base.axes(a::TwoLevelBlockedUnitRange) = axes(a.global_range)
Base.Broadcast.axistype(a::T, b::T) where T<:TwoLevelBlockedUnitRange =  Base.Broadcast.axistype(a.global_range,b.global_range)
Base.Broadcast.axistype(a::TwoLevelBlockedUnitRange, b::TwoLevelBlockedUnitRange) = Base.Broadcast.axistype(a.global_range,b.global_range)
Base.Broadcast.axistype(a::TwoLevelBlockedUnitRange, b) = Base.Broadcast.axistype(a.global_range,b)
Base.Broadcast.axistype(a, b::TwoLevelBlockedUnitRange) =Base.Broadcast.axistype(a,b.global_range)
Base.print_matrix_row(
  io::IO,
  X::TwoLevelBlockedUnitRange,
  A::Vector,
  i::Integer,
  cols::AbstractVector,
  sep::AbstractString) = print_matrix_row(io,X.global_range,A,i,cols,sep)
Base.show(io::IO, mimetype::MIME"text/plain", a::TwoLevelBlockedUnitRange) = show(io,mimetype,a.global_range)


text:
function parse_disjuncts(head)
    ret = MacroTools.postwalk(head) do ex
        if @capture(ex, num_::call_(args__))
            quote $call($(args...), p = $num) end
        else
            ex
        end
    end
    ret
end

function parse_call(call, names)
    names = map(names) do name
        if name isa QuoteNode
            quote Constant(String($name)) end
        else
            name
        end
    end
    return quote sp.add_fact($call($(names...))) end
end

function parse_body(body)
    if @capture(body, head_ :- tail_)
        disjuncts = parse_disjuncts(head)
        return quote sp.add_clause(AnnotatedDisjunction($disjuncts, $tail)) end

    elseif @capture(body, head_ << tail_)
        return quote sp.add_clause($head << $tail) end

    elseif @capture(body, call_(names__))
        return parse_call(call, names)

    else
        body
    end
end

function _problox(expr)
    if @capture(expr, function name_(args__) body__ end)
        trans = map(body) do ex
            if @capture(ex, begin body__; end)
                length(body) == 1 || error("Parsing error at $body.")
                parse_body(body[1])
            else
                error("Parsing error at $body.")
            end
        end
        expr = quote
            function $name($(args...))
                sp = SimpleProgram()
                $(trans...)
                sp
            end
        end

    elseif @capture(expr, begin body__ end)
        trans = map(body) do ex
            if @capture(ex, begin body__; end)
                length(body) == 1 || error("Parsing error at $body.")
                parse_body(body[1])
            else
                error("Parsing error at $body.")
            end
        end
        expr = quote
            sp = SimpleProgram()
            $(trans...)
            sp
        end
    end
    expr
end

macro problox(expr)
    new = _problox(expr)
    new = MacroTools.postwalk(rmlines ∘ unblock, new)
    esc(new)
end


text:
using Rebugger, Test

include("edit.jl")
include("interpret.jl")
if Sys.isunix() && VERSION >= v"1.1.0"
    include("interpret_ui.jl")
else
    @warn "Skipping UI tests"
end
println("done")  # there is so much terminal manipulation, best to let the user know


text:
using Flase, Test

@testset "Dog" begin
    dog = Flase.Dog( [1.,2.], zeros(2) )
    @test typeof(dog) == Flase.Dog{Float64}
    @test Flase.Dog( dog; state = Flase.OCCUPIED ).state == Flase.OCCUPIED
end # testset


text:
using Omega: withkernel, kseα

"Replica Exchange (Parallel Tempering)"
struct ReplicaAlg <: SamplingAlgorithm end

"Single Site Metropolis Hastings"
const Replica = ReplicaAlg()

softhard(::Type{ReplicaAlg}) = IsSoft{ReplicaAlg}()

defΩ(::ReplicaAlg) = Omega.LinearΩ{ID, UnitRange{Int64}, Vector{Real}}
defΩ(x, ::ReplicaAlg; inneralg...) = defΩ(inneralg)

isapproximate(::ReplicaAlg) = true

function swap!(v, i, j)
  temp = v[i]
  v[i] = v[j]
  v[j] = temp
end

"Swap adjacent chains"
function exchange!(rng, logdensity, ωs, temps, kernel)
  for i in length(ωs):-1:2
    j = i - 1
    E_i_x = withkernel(kernel(temps[i])) do
      logdensity(ωs[i])
    end
    E_j_x = withkernel(kernel(temps[j])) do
      logdensity(ωs[i])
    end
    E_i_y = withkernel(kernel(temps[i])) do
      logdensity(ωs[j])
    end
    E_j_y = withkernel(kernel(temps[j])) do
      logdensity(ωs[j])
    end
    k = (E_i_y + E_j_x) - (E_i_x + E_j_y)
    doswap = log(rand(rng)) < k
    if doswap
      swap!(ωs, i, j)
    end
  end
end

"Logarithmically spaced temperatures"
logtemps(n, k = 10) = exp.(k * range(-2.0, stop = 1.0, length = n))

struct PreSwap end
struct PostSwap end
struct InLoop end

"""Sample from `density` using Replica Exchange

$(SIGNATURES)

Replica exchange (aka parallel tempemring) runs `nreplicas` independent mcmc
chains in parallel.
Returns samples from lowest temperature chain

# Arguments
- `ΩT`:
- `logdensity`: Real-valued `RandVar`
- `n`: Number of samples
- `swapevery` : performs swap every swapevery iterations
- `nreplicas` : number of replica chains to run
- `temps` : temperatures of different chains
- `inneralg` : Algorithm used for each chain
- `algargs::NamedTuple` : keyword arguments to be passed to `inneralg` in:
   `rand(ΩT, density, swapevery, inneralg; algargs...)`
- `kernel`: Kernel to use for soft constraints (DEPRECATE ME)

# Returns

"""
function Base.rand(rng,
                   ΩT::Type{OT},
                   logdensity::RandVar,
                   n::Integer,
                   alg::ReplicaAlg;
                   inneralg = SSMH,
                   algargs = NamedTuple(),
                   swapevery = 1,
                   nreplicas = 4,
                   temps = logtemps(nreplicas),
                   kernel = Omega.kseα) where {OT <: Ω}
  @pre issorted(temps)
  @pre n % swapevery == 0
  @pre nreplicas == length(temps)
  @show temps

  ωsamples = OT[]
  ωs = [ΩT() for i = 1:nreplicas]

  # Do swapevery steps for each chain, then swap ωs
  lowesti = nreplicas
  for j = 1:div(n, swapevery)
    for i = 1:nreplicas
      withkernel(kernel(temps[i])) do
        try
          ωst = rand(rng, ΩT, logdensity, swapevery, inneralg;
                     ωinit = ωs[i],
                     offset = (j - 1) * swapevery,
                     algargs...)
          if i == lowesti # keep lowest temperatre
            append!(ωsamples, ωst)
          end
          lens(InLoop, (ωs = ωst, j = j, i = i, temp = temps[i]))
          ωs[i] = ωst[end]
        catch e
          # rethrow(e)
          println("Chain at temp $(temps[i]) Failed due to:", e)
        end
      end
    end
    lens(PreSwap, (ωs = ωs, temps = temps))
    exchange!(rng, logdensity, ωs, temps, kernel)
    lens(PostSwap, (ωs = ωs, temps = temps))
  end
  ωsamples
end


text:
using Omega
using Spec
using Test
using Pkg

# FIXME: Make this a package
# include("TestLib.jl")

# Add TestModels as submodule
# Pkg.develop(PackageSpec(url=joinpath(dirname(pathof(Omega)), "..", "test", "OmegaTestModels")))

# walktests(Omega, exclude = ["typestable.jl", "simple.jl"])

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(Tuple{Core.kwftype(typeof(ITensors.eigen)),NamedTuple{(:ishermitian, :which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho),Tuple{Bool,Nothing,TagSet,Int64,Int64,Float64,Nothing,String}},typeof(eigen),ITensor{2},IndexSet{1,Index{Int64},Tuple{Index{Int64}}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}})
    precompile(Tuple{Core.kwftype(typeof(ITensors.factorize)),NamedTuple{(:which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho),Tuple{Nothing,TagSet,Int64,Int64,Float64,Nothing,String}},typeof(factorize),ITensor{3},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}})
    precompile(Tuple{Core.kwftype(typeof(ITensors.replacebond!)),NamedTuple{(:maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho, :which_decomp),Tuple{Int64,Int64,Float64,Nothing,String,Nothing}},typeof(replacebond!),MPS,Int64,ITensor{3}})
    precompile(Tuple{Core.kwftype(typeof(KrylovKit.eigsolve)),NamedTuple{(:ishermitian, :tol, :krylovdim, :maxiter),Tuple{Bool,Float64,Int64,Int64}},typeof(KrylovKit.eigsolve),ProjMPO,ITensor{3},Int64,Symbol})
    precompile(Tuple{Core.kwftype(typeof(KrylovKit.expand!)),NamedTuple{(:verbosity,),Tuple{Int64}},typeof(KrylovKit.expand!),KrylovKit.LanczosIterator{ProjMPO,ITensor{3},KrylovKit.ModifiedGramSchmidt2},KrylovKit.LanczosFactorization{ITensor{3},Float64}})
    precompile(Tuple{Core.kwftype(typeof(NDTensors.eigen)),NamedTuple{(:ishermitian, :which_decomp, :tags, :maxdim, :mindim, :cutoff, :eigen_perturbation, :ortho),Tuple{Bool,Nothing,TagSet,Int64,Int64,Float64,Nothing,String}},typeof(eigen),LinearAlgebra.Hermitian{Float64,NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}}}})
    precompile(Tuple{Core.kwftype(typeof(NDTensors.svd)),NamedTuple{(:alg,),Tuple{String}},typeof(svd),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}}})
    precompile(Tuple{Core.kwftype(typeof(NDTensors.svd)),NamedTuple{(:alg,),Tuple{String}},typeof(svd),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}}})
    precompile(Tuple{ProjMPO,ITensor{3}})
    precompile(Tuple{Type{MPO},AutoMPO,Array{Index{Int64},1}})
    precompile(Tuple{Type{TagSet},String})
    precompile(Tuple{typeof(*),Float64,NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}}})
    precompile(Tuple{typeof(*),ITensor{2},ITensor{3}})
    precompile(Tuple{typeof(*),ITensor{3},ITensor{1}})
    precompile(Tuple{typeof(*),ITensor{3},ITensor{2}})
    precompile(Tuple{typeof(*),ITensor{3},ITensor{3}})
    precompile(Tuple{typeof(*),ITensor{4},ITensor{1}})
    precompile(Tuple{typeof(*),ITensor{4},ITensor{4}})
    precompile(Tuple{typeof(==),Array{Index{Int64},1},Array{Index{Int64},1}})
    precompile(Tuple{typeof(==),IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}})
    precompile(Tuple{typeof(==),TagSet,TagSet})
    precompile(Tuple{typeof(Base.Broadcast.materialize!),ITensor{4},Base.Broadcast.Broadcasted{ITensors.ITensorStyle,Nothing,typeof(+),Tuple{ITensor{4},ITensor{4}}}})
    precompile(Tuple{typeof(Base.afoldl),typeof(*),ITensor{4},ITensor{3}})
    precompile(Tuple{typeof(ITensors.setstore!),ITensor{2},NDTensors.Dense{Float64,Array{Float64,1}}})
    precompile(Tuple{typeof(KrylovKit.orthogonalize!),ITensor{3},ITensor{3},KrylovKit.ModifiedGramSchmidt2})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Number,2,NDTensors.Combiner,IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,0,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{0,Union{},Tuple{}}},Tuple{},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Number,2,NDTensors.Combiner,IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},Tuple{Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},Tuple{Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(NDTensors.contract!!),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64}})
    precompile(Tuple{typeof(NDTensors.contraction_output),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},IndexSet{0,Union{},Tuple{}}})
    precompile(Tuple{typeof(NDTensors.contraction_output),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}})
    precompile(Tuple{typeof(NDTensors.getperm),IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}})
    precompile(Tuple{typeof(NDTensors.permutedims!!),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},Function})
    precompile(Tuple{typeof(NDTensors.permutedims!!),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},Function})
    precompile(Tuple{typeof(Random.randn!),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}}})
    precompile(Tuple{typeof(axpy!),Float64,ITensor{3},ITensor{3}})
    precompile(Tuple{typeof(combiner),Index{Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Number,2,NDTensors.Combiner,IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NTuple{4,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Number,2,NDTensors.Combiner,IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Number,2,NDTensors.Combiner,IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,2,NDTensors.Diag{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},Tuple{Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},Tuple{Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NTuple{4,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},Tuple{Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},Tuple{}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},Tuple{Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,1,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}},Tuple{Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,3,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}}},Tuple{Int64,Int64,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NTuple{4,Int64}})
    precompile(Tuple{typeof(contract),NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64},NDTensors.Tensor{Float64,4,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{4,Index{Int64},NTuple{4,Index{Int64}}}},NTuple{4,Int64}})
    precompile(Tuple{typeof(hassameinds),ITensor{3},Tuple{Index{Int64},Index{Int64},Index{Int64}}})
    precompile(Tuple{typeof(hastags),TagSet,String})
    precompile(Tuple{typeof(itensor),Array{Float64,2},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(itensor),Array{Float64,3},Index{Int64},Index{Int64},Vararg{Index,N} where N})
    precompile(Tuple{typeof(mul!),ITensor{3},ITensor{3},Float64})
    precompile(Tuple{typeof(mul!),ITensor{3},KrylovKit.OrthonormalBasis{ITensor{3}},SubArray{Float64,1,Array{Float64,2},Tuple{UnitRange{Int64},Int64},true}})
    precompile(Tuple{typeof(prime),ITensor{3},ITensors.Not{TagSet}})
    precompile(Tuple{typeof(prime),IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},Int64,ITensors.Not{IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}}})
    precompile(Tuple{typeof(randomITensor),Index{Int64},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(randomMPS),Array{Index{Int64},1},Int64})
    precompile(Tuple{typeof(replaceind!),ITensor{2},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(replaceind!),ITensor{3},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(replaceind),IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},Index{Int64},Index{Int64}})
    precompile(Tuple{typeof(replaceinds),IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}},Tuple{Index{Int64},Index{Int64}},Tuple{Index{Int64},Index{Int64}}})
    precompile(Tuple{typeof(replaceinds),IndexSet{3,Index{Int64},Tuple{Index{Int64},Index{Int64},Index{Int64}}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}},IndexSet{1,Index{Int64},Tuple{Index{Int64}}}})
    precompile(Tuple{typeof(setindex!),Dict{Array{ITensors.SiteOp,1},Array{Float64,2}},Array{Float64,2},Array{ITensors.SiteOp,1}})
    precompile(Tuple{typeof(svd),NDTensors.Tensor{Complex{Float64},2,NDTensors.Dense{Complex{Float64},Array{Complex{Float64},1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}}})
    precompile(Tuple{typeof(svd),NDTensors.Tensor{Float64,2,NDTensors.Dense{Float64,Array{Float64,1}},IndexSet{2,Index{Int64},Tuple{Index{Int64},Index{Int64}}}}})
    precompile(Tuple{typeof(zip),IndexSet{1,Index{Int64},Tuple{Index{Int64}}},Vararg{IndexSet{1,Index{Int64},Tuple{Index{Int64}}},N} where N})
end


text:
using Plots
using LaTeXStrings

function plot_example_abs(f, u, x, data, û_FDM, û)
    x_opt = range(x[1], x[end]; length=1000)
    f_opt = f.(x_opt)
    u_opt = u.(x_opt)

    append!(û_FDM, NaN)

    pf = plot(x_opt, f_opt; c=:grey, ylabel=L"f", label=L"f")
    pf = plot!(x, data; c=:2, ylabel=L"f", label=L"f_{noisy}", legend=:bottomright)
    pu = plot(x_opt, u_opt; c=:grey, ylabel=L"u", label=L"u")
    pu = plot!(
        x,
        [û_FDM, û];
        c=[:2 :1],
        ylabel=L"u",
        label=[L"\hat{u}_{FDM}" L"\hat{u}_{tvdiff}"],
        legend=:bottomright,
    )

    return plot(pf, pu; layout=(2, 1), show=true)
end

function plot_FDM(f)
    n = length(f)
    dx = 1 / (n - 1)

    û_FDM = diff(f) / dx  # FDM
    return plot_FDM(f, û_FDM)
end

function plot_FDM(f, û)
    append!(û, NaN)

    n = length(f)
    x = range(0, 1; length=n)

    pf = plot(x, f; c=:1, ylabel=L"f")
    pu = plot(x, û; c=:1, ylabel=L"\hat{u}_{FDM}")

    return plot(pf, pu; layout=(2, 1), legend=false, show=true)
end

function plot_demo_large_diff(f)
    dx = 1
    û = diff(f) / dx  # FDM
    append!(û, NaN) # hide

    pf = plot(f; c=:1, ylabel=L"f (L/min)")
    pu = plot(û; c=:1, ylabel=L"\hat{u}_{FDM} (L/min/s)", xlabel=L"t (s)")

    return plot(pf, pu; layout=(2, 1), legend=false, show=true, dpi=300)
end

function plot_demo_large_tvdiff(f, û)
    append!(û, NaN) # hide

    pf = plot(f; c=:1, ylabel=latexstring("f (L/min)"))
    pu = plot(û; c=:1, ylabel=L"\hat{u}_{tvdiff} (L/min/s)", xlabel=L"t (s)")

    return plot(pf, pu; layout=(2, 1), legend=false, show=true, dpi=300)
end

function plot_tvdiff(f, û)
    n = length(f)
    dx = 1 / (n - 1)

    return plot_tvdiff(f, û, dx)
end

"""
Plot data f and total variance regularized
numerical differences `u`.
"""
function plot_tvdiff(f, û, dx::Real)
    n = length(f)
    x = range(0; step=dx, length=n)

    pf = plot(x, f; ylabel=L"f")
    pu = plot(x, û; ylabel=L"\hat{u}_{tvdiff}")

    return plot(pf, pu; layout=(2, 1), legend=false, show=true)
end

function plot_tvdiff_all(f, û_TVR)
    n = length(f)
    dx = 1 / (n - 1)
    û_FDM = diff(f) / dx  # FDM

    return plot_tvdiff_all(f, û_FDM, û_TVR)
end

function plot_tvdiff_all(f, û_FDM, û_TVR)
    pf = plot(f; c=:1, ylabel=L"f")
    pu1 = plot(û_FDM; c=:1, ylabel=L"\hat{u}_{FDM}")
    pu2 = plot(û_TVR; c=:1, ylabel=L"\hat{u}_{tvdiff}")

    return plot(pf, pu1, pu2; layout=(3, 1), legend=false, show=true, dpi=300)
end


text:
### Creating destination popularity: INPUT is datamatrix, e, T, number_zones and OUTPUT is p_dest

function createpdestin(datamatrix, number_zones)
    p_dest = zeros(number_zones, number_zones, T);
    max_x = zeros(number_zones, number_zones);
    min_x = zeros(number_zones, number_zones);
    normalization_factor = zeros(number_zones, T);

    ## Calculating the minima and maxima of mean travel time
    for i=1:number_zones
        # Calculating the minima and maxima
        for j=1:number_zones # destination zone ID
            max_x[i,j] = maximum(datamatrix[i, j, :, 1])
            min_x[i,j] = minimum(datamatrix[i, j, :, 1])
        end
    end

    ## Calculating the unnormalized p_dest
    for i=1:number_zones
        for j=1:number_zones
            for t=1:T
                if (max_x[i, j] > 0)
                    mean = datamatrix[i, j, t, 1];
                    p_dest[i, j, t] = ((mean - min_x[i, j])/ (max_x[i, j] - min_x[i, j]))^e_dest;
                end
            end
        end
    end

    ## Calculating normalization factors
    for i=1:number_zones
        for t=1:T
            normalization_factor[i, t] = sum(p_dest[i, :, t]);
        end
    end

    ## Normalizing the probability distributions of p_dest
    for i=1:number_zones
        for t=1:T
            if(normalization_factor[i, t] > 0) # A zero would lead to infinity given the formula below and create NaN entries.
                for j=1:number_zones
                    p_dest[i, j, t] = p_dest[i, j, t] / normalization_factor[i, t];
                end
            end
        end
    end

    return p_dest

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
# surrogate_ground_truth_nodal_mala-logistic
#
# This script runs a long chain based on the
# Metropolis-adjusted Langevin algorithm (MALA). This is used as a gold
# standard as this procedure yields an asymptotically unbiased chain
# and thus is used to compare the efficiency of the biased chains produced
# in the experiment bias_variance-arwmh-logistic.
#
# The target distribution is a Bayesian logistic regression posterior
# under independent Gaussian priors based on a prostate cancer dataset with
# binary outcomes indicating whether cancer has spread to surrounding lymph
# nodes.

using Distributions: Normal
using RDatasets
using StatsFuns: logistic
using SteinDistributions: SteinLogisticRegressionGaussianPrior
using SteinSamplers: runmala

include("experiment_utils.jl")

# Set random number generator seed
@parseintcli seed "s" "seed" 7
# setup the data
data_set = "nodal"
X = y = nothing
nodal = dataset("boot", "nodal")
y = 2 * convert(Array{Float64},nodal[:R]) - 1
# X includes an intercept term
X = convert(Array{Float64, 2}, nodal[[:M, :Aged, :Stage, :Grade, :Xray, :Acid]])

# N is number of samples, p is dimension of parameter
N, d = size(X)
# the max number of samples to run in the solver
n = 10000000
# setup distribution
target = SteinLogisticRegressionGaussianPrior(X, y)
# name of prior + posterior
distname = "logisticgaussianprior"
# intiial guess for beta
beta0 = zeros(d)
# the epsilons to run the experiemnt
epsilons = [1e-2]

for epsilon = epsilons
    # Draw MALA samples
    (betas_mala, numgrad_mala) =
        @setseed runmala(target, beta0;
                         epsilonfunc=t -> epsilon,
                         numiter=n)
    # estimate the final beta moments
    betahat = vec(mean(betas_mala, 1))
    betavcov = cov(betas_mala)
    betasd = sqrt(diag(betavcov))
    # compute the running estimated mu (i.e. P(Y_i = 1 | X_i))
    # when using the first k values from the chain
    ks = round(Int, linspace(0, n, 1001))[2:end]
    meanpredprobs = Array(Float64, N, length(ks))
    for (i,k) in enumerate(ks)
        meanbeta = mean(betas_mala[1:k,:], 1)
        predprobs = logistic(X * meanbeta')
        meanpredprobs[:,i] = vec(predprobs)
    end
    # scale the mu's by the l_infty norm of each X_i
    scaledpredprobs = meanpredprobs ./ vec(maximum(X, 2))
    # fetch constants for computing confidence intervals of beta
    phi = Normal()
    zscore = quantile(phi, 1 - 0.05/2)
    # Package and save results
    summary_data = Dict{Any, Any}(
        "n" => n,
        "seed" => seed,
        "distname" => distname,
        "epsilon" => epsilon,
        "beta.hat" => betahat,
        "beta.vcov" => betavcov,
        "beta.sd" => betasd,
        "beta.hat.low" => betahat - zscore .* betasd,
        "beta.hat.high" => betahat + zscore .* betasd,
        "scaled.pred.probs" => scaledpredprobs,
        "scaled.pred.probs.i" => ks,
        "numgrad" => numgrad_mala,
        "sampler" => "MALA-long-nodal"
    )
    save_json(
        summary_data;
        dir="bias_variance-arwmh-logistic-MALA",
        distname=distname,
        dataset=data_set,
        n=n,
        epsilon=epsilon,
        d=d,
        seed=seed
    )
end


text:
factor_roots_comb(p::AbstractPolynomialLike) = factor_main(p, factor_sqfr_roots)
factor_roots_comb(eq) = wrap(factor_roots_comb, eq)

factor_roots_SSP(p::AbstractPolynomialLike) = factor_main(p, factor_sqfr_SSP)
factor_roots_SSP(eq) = wrap(factor_roots_SSP, eq)

###################### Common Factoring ************************************

function factor_main(p::AbstractPolynomialLike, fs...; resolution=1e-5, abstol=1e-10)
    c₀, p = integer_poly_coef(p)
    x = var(p)

    P = FactoredPoly()
    lc = leading(p)

    for w in factors((decompose(p)))
        v, k = first(w), last(w)
        F, G, lc = combine_factorization_algs(v, lc, fs...; abstol, resolution)

        for q in [F; G]
            add_factor!(P, prim(q), k)
        end
    end

    add_coefficient!(P, p, c₀)

    P
end

function add_coefficient!(P, p::AbstractPolynomialLike, c₀)
    c, ρ = integer_poly_coef(p ÷ poly(P))
    !isone(prim(ρ)) && add_factor!(P, prim(ρ), 1)
    c *= cont(ρ) * c₀
    c = denominator(c) == 1 ? numerator(c) : c
    !isone(c) && add_factor!(P, c, 1)
end

function combine_factorization_algs(p::AbstractPolynomialLike, lc, fs...; resolution=0, abstol=1e-10)
    F₁, G₁, lc = fs[1](p, lc; resolution, abstol)

    if isempty(G₁) || length(fs) == 1
        return F₁, G₁, lc
    end

    # fs = (length(G₁) == 1 ? fs[2:end] : fs)
    fs = fs[2:end]
    G = []
    F = F₁

    for g in G₁
        F₂, G₂, lc = combine_factorization_algs(g, lc, fs...; resolution, abstol)
        append!(F, F₂)
        append!(G, G₂)
    end

    return F, G, lc
end

###############################################################################

function factor_sqfr_roots(p::AbstractPolynomialLike, lc; max_deg=100000, max_misses=10000, resolution=1e-5, abstol=1e-10)
    x = var(p)
    v = polynomial(Float64.(coefficients(p)), terms(p))
    r, s = find_roots(v, x; abstol)
    s = s[1:2:end]
    σ = [r; 2*real.(s)]
    μ = [r; abs2.(s)]
    n = length(μ)

    d = 1
    misses = 0
    F = []
    G = []

    test_fun = function(pat)
        l = [i for i=1:n if testbit(pat, i)]
        S = sum(σ[l]; init=0.0)
        M = prod(μ[l]; init=1.0)

        if near_integer(S*lc) && near_integer(M*lc)
            q = simple_factor(x, σ, μ, l)
            q = integer_poly(q*lc)
            g = gcd(lc, coefficients(q)...)
            q = q ÷ g

            if iszero(p % q)
                push!(F, q)
                p = remove_factor(p, q)
                lc = g
                d *= g
                return true
            else
                misses += 1
                if misses == max_misses
                    error("too many misses!")
                end
            end
        end
        false
    end

    m = length(r) + 2*length(s)
    mask = 0
    cmask = 2^length(μ) - 2^length(r)

    try
        for k = 1:min(m÷2, max_deg)
            n₁ = n - count_ones(mask)
            if n₁*log(n₁) - k*log(k) - (n₁-k)*log(n₁-k) < 15
                mask = traverse_patterns(0, mask, n, k, cmask, test_fun)
            end
        end
    catch e
    end

    deg(p) > 0 && push!(G, p)
    F, G, lc
end

###############################################################################

near_integer(x::Float64; abstol=1e-5) = abs(x - round(x)) < abstol
near_integer(x::BigFloat; abstol=1e-10) = abs(x - round(x)) < abstol
exclude(v, l) = [v[i] for i=1:length(v) if i ∉ l]
testbit(x, n) = isodd(x >> (n-1))
expand_binary(x::Int) = [i for i=1:64-leading_zeros(x) if testbit(x,i)]
expand_binary(x::BigInt) = [i for i=1:128-leading_zeros(x) if testbit(x,i)]
expand_binary(x, n) = [i for i=1:n if testbit(x,i)]

# generates all n-bit numbers with d 1 bits
# excluding mask bits
# test each pattern by calling fun
function traverse_patterns(pat, mask, n, k, cmask, fun)
    bit = 1

    for i = 1:min(trailing_zeros(pat), n)
        pat₁ = pat | bit
        if (pat != pat₁) && (mask & pat₁ == 0)
            nb = count_ones(pat₁) + count_ones(pat₁ & cmask)
            if nb == k
                if fun(pat₁)
                    mask |= pat₁
                end
            elseif nb < k
                mask = traverse_patterns(pat₁, mask, n, k, cmask, fun)
            end
        end
        bit <<= 1
    end
    return mask
end

function simple_factor(x, σ, μ, l)
    prod(μ[i]==σ[i] ? (x - σ[i]) : x^2 - σ[i]*x + μ[i] for i in l; init=one(x))
end

simple_factor(x, σ, μ, i::Integer) = simple_factor(x, σ, μ, [i])

##############################################################################

frac(x) = x - floor(x)
index(x, L) = 1 + floor(Int, L * x)

struct Context{T, S}
    failures::BitMatrix
    σ::Vector{T}
    L::Int
    ls::Vector{S}
    abstol::T
    cmask::S
    cutoff::Int
end

Base.size(ctx::Context) = size(ctx.failures)
index(ctx::Context, t) = index(t, ctx.L)
add_path!(ctx::Context, path) = push!(ctx.ls, path)
@inbounds is_failure(ctx::Context, k, i) = ctx.failures[k,i]
@inbounds fail(ctx::Context, k, i) = (ctx.failures[k,i] = true)
can_grow(ctx::Context, path) = count_ones(path) + count_ones(path & ctx.cmask) < ctx.cutoff

"""
    solve_SSP is the main SSP-solving search function
    solve_SSP implements a pseudo-polynomial time dynamic programming algorithm
    the actual search phase is done by traverse_SSP
"""
function solve_SSP(σ::Array{T}, μ::Array{T}; resolution=1e-5) where T<:Real
    L = round(Int, 1 / resolution)
    σ = frac.(σ)
    m = L * round(Int, max(sum(σ),one(T)))
    n = length(σ)

    failures = falses(n, m)
    S = (n <= 64 ? Int : Int128)
    path = zero(S)
    ls = S[]

    cmask = sum(one(S)<<(i-1) for i=1:n if μ[i] != σ[i]; init=zero(S))
    cutoff = (n + count_ones(cmask)) ÷ 2

    ctx = Context{T,S}(failures, σ, L, ls, resolution, cmask, cutoff)
    traverse_SSP(ctx, 1, zero(T), path)
    return ctx.ls
end

"""
    traverse_SSP is the inner loop of the SSP-solving algorithm
    this function consumes the bulk of time of factorization and needs to be fine-tuned!
"""
function traverse_SSP(ctx::Context, k, t, path)
    n, m = size(ctx)

    if k == n + 1
        if 0 < count_ones(path) < n && near_integer(t; abstol=ctx.abstol)
            add_path!(ctx, path)
            return true
        end
        return false
    end

    i = index(ctx, t)

    if i > m return false end

    if is_failure(ctx, k, i) return false end

    if can_grow(ctx, path)
        success = traverse_SSP(ctx, k+1, t+ctx.σ[k], path | (1 << (k-1))) |
                  traverse_SSP(ctx, k+1, t, path)
    else
        success = traverse_SSP(ctx, k+1, t, path)
    end

    if !success
        fail(ctx, k, i)
    end

    success
end

##############################################################################

"""
    `factor_sqfr_SSP` factors a square-free polynomial `p` using `solve_SSP`
    `lc` is the leading coefficient
"""
function factor_sqfr_SSP(p::AbstractPolynomialLike, lc; abstol=1e-10, resolution=0)
    factor_sqfr_SSP(BigFloat, p, lc; abstol, resolution)
end

function factor_sqfr_SSP(T, p::AbstractPolynomialLike, lc; abstol=1e-10, resolution=0)
    x = var(p)
    p = rationalize(p) / leading(p)
    lc = T(1.0 / cont(p))
    p = prim(p)

    G = find_factors(T, x, p, lc; resolution)
    sort!(G; by=deg)
    G, [], lc
end

"""
    `find_factors` find candidate factors of `p`
"""
function find_factors(T, x, p, lc; resolution=1e-5)
    F = find_factor(T, x, p, lc; resolution)
    length(F) == 1 && return F

    G = []
    for q in F
        Q = find_factors(T, x, q, lc; resolution)
        append!(G, Q)
    end
    return G
end

"""
    `find_factor` is the core of the SSP factorization algorithm
"""
function find_factor(T, x, p, lc; abstol=1e-10, resolution=0)
    if resolution == 0
        resolution = 1e-3 / sqrt(deg(p))
    end

    v = polynomial(T.(coefficients(p)), terms(p))
    r, s = find_roots(T, v, x; abstol)
    σ = T.([r; 2*real.(s[1:2:end])])
    μ = T.([r; abs2.(s[1:2:end])])
    n = length(σ)

    w = sortperm(frac.(σ); rev=true)
    σ = σ[w]
    μ = μ[w]

    ls = solve_SSP(σ*lc, μ*lc; resolution)
    ls = sort(ls; by=count_ones)

    d = zero(eltype(ls))
    F = []
    x₀ = find_test_x(p)
    α = round(BigInt, lc*p(x₀))
    # α = lc*p(im)

    for c in ls
        c = c ⊻ d
        if c > 0
            l = expand_binary(c, n)
            S = sum(σ[l]; init=0.0)

            if near_integer(S*lc)
                β = round(BigInt, lc*simple_factor(x₀, σ, μ, l))
                if !iszero(β) && !iszero(α % β) continue end
                # β = lc*simple_factor(im, σ, μ, l)
                # if !iszero(β)
                #     αβ = α * conj(β) / abs2(β)
                #     if !near_integer(real(αβ)) || !near_integer(imag(αβ))
                #         continue
                #     end
                # end

                q = simple_factor(x, σ, μ, l)
                if maximum(coefficients(q)) < typemax(Int)/10
                    q = prim(integer_poly(q*lc))
                    g = q #gcd(p, q)

                    if !isone(g) && iszero(p % g)
                        d |= c
                        p = prim(remove_factor(p, g))
                        push!(F, g)
                    end
                end
            end
        end
    end

    !isone(p) && push!(F, p)
    return F
end

"""
    `find_test_x` finds the integer with the smallest absolute value as a test number
    note that 1 is used implicitely by the core algorithm
    `find_test_x` searches integers in -1, 2, -2, 3, -3,... list and returns the
    first one that results in a nonzero `p`
"""
function find_test_x(p)
    i = 1
    while true
        x₀ = BigInt(1 + i÷2)
        x₀ = isodd(i) ? -x₀ : x₀
        if p(x₀) != 0 return x₀ end
        i += 1
    end
end


text:
# Base type for all experiments
import Base: run
abstract type BanditExpBase end

function run( experiment::BanditExpBase,
              noOfTimeSteps::Integer,
              noOfRounds::Integer )
    error( "[BanditExpBase]: No Implementation of run() for ", typeof(experiment) )
end


text:
module MainPaths
using Graphs
using SparseArrays

export mainpath
export SPCEdge, SPCVertex, GKP
export ForwardLocal, BackwardLocal, ForwardBackwardLocal

include("utils.jl")
include("weights_spc.jl")
include("weights_gkp.jl")
include("traversal.jl")
include("mainpath.jl")
include("segments.jl")

end


text:
using TimedSound
using Base.Test

include("test_sound_construction.jl")
include("test_stream_construction.jl")


text:
@benchgroup "connectivity" begin
  @benchgroup "digraphs" begin
    for (name, g) in DIGRAPHS
      @bench "$(name): strongly_connected_components" LightGraphs.strongly_connected_components($g)
    end
  end # digraphs
  @benchgroup "graphs" begin
    for (name, g) in GRAPHS
        @bench "$(name): connected_components" LightGraphs.connected_components($g)
    end
  end # graphs
end # connectivity


text:
using GraphRecipes
using Plots
using Test

include("readme.jl")

@testset "utils.jl" begin

    @test GraphRecipes.directed_curve(0., 1., 0., 1.) == GraphRecipes.directed_curve(0, 1, 0, 1)

    @testset "Functions from Plots.jl" begin

        @test GraphRecipes.isnothing(nothing) == Plots.isnothing(nothing)
        @test GraphRecipes.isnothing(missing) == Plots.isnothing(missing)
        @test GraphRecipes.isnothing(NaN) == Plots.isnothing(NaN)
        @test GraphRecipes.isnothing(0) == Plots.isnothing(0)
        @test GraphRecipes.isnothing(1) == Plots.isnothing(1)
        @test GraphRecipes.isnothing(0.0) == Plots.isnothing(0.0)
        @test GraphRecipes.isnothing(1.0) == Plots.isnothing(1.0)

        for (s, e) in [ (rand(), rand()) for i in 1:100 ]
            @test GraphRecipes.partialcircle(s, e) == Plots.partialcircle(s, e)
        end

        @testset "nearest_intersection" begin
            @test GraphRecipes.nearest_intersection(0, 0, 3, 3, [(1,0), (0,1)]) == (0,0,0.5,0.5)
            @test GraphRecipes.nearest_intersection(1, 2, 1, 2, []) == (1, 2, 1, 2)
        end

        @testset "unoccupied_angle" begin
            @test GraphRecipes.unoccupied_angle(1, 1, [1, 1, 1, 1], [2, 0, 3, -1]) == 2pi
        end

        @testset "islabel" begin
            @test GraphRecipes.islabel("hi") == true
            @test GraphRecipes.islabel(1) == true
            @test GraphRecipes.islabel(missing) == false
            @test GraphRecipes.islabel(NaN) == false
            @test GraphRecipes.islabel(false) == false
            @test GraphRecipes.islabel("") == false
        end

        @testset "control_point" begin
            @test GraphRecipes.control_point(0, 0, 6, 0, 4) == (4,3)
        end

    end

end

# -----------------------------------------
# marginalhist

# using Distributions
# n = 1000
# x = rand(Gamma(2), n)
# y = -0.5x + randn(n)
# marginalhist(x, y)

# -----------------------------------------
# portfolio composition map

# # fake data
# tickers = ["IBM", "Google", "Apple", "Intel"]
# N = 10
# D = length(tickers)
# weights = rand(N,D)
# weights ./= sum(weights, 2)
# returns = sort!((1:N) + D*randn(N))

# # plot it
# portfoliocomposition(weights, returns, labels = tickers')

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


text:
"""

In this example, the centers of mass of the protein and of the POPC residues
are computed, and the difference between the z-coordinates of them is ploted.

"""

using M3GTools

println(" Loading Plots... ")
using Plots
ENV["GKSwstype"]="nul" # This supresses the need of a display while plotting

# Initialize simulation data and path to VMD executable

println(" Reading simulation data ... ")

mysim = M3GTools.init(psf="../simulation_files/structure.psf",
                      dcd="../simulation_files/structure.dcd",
                      vmd="vmd")

println(" Defining selections... ")
popc = M3GTools.select(mysim,"resname POPC")
prot = M3GTools.select(mysim,"protein")

println(" Computing center of masses... ")
popc_cm = Array{Float32}(undef,mysim.nframes,3)
prot_cm = Array{Float32}(undef,mysim.nframes,3)
cm = Array{Float32}(undef,3)

for i in 1:mysim.nframes
  sides, x, y, z = M3GTools.nextframe(mysim)

  cm = M3GTools.cm(popc,mysim,x,y,z)
  for j in 1:3
    popc_cm[i,j] = cm[j]
  end
  cm = M3GTools.cm(prot,mysim,x,y,z)
  for j in 1:3
    prot_cm[i,j] = cm[j]
  end

end
M3GTools.closedcd(mysim)

# Printing the difference in the Z-coordinate of the center of
# masses of the selections

println(" Plotting... ")
diff = Vector{Float32}(undef,mysim.nframes)
for i in 1:mysim.nframes
  diff[i] = prot_cm[i,3] - popc_cm[i,3]
end
x = [ i for i in 1:mysim.nframes ]
plot(x,diff,xlabel="Frame",ylabel="Z-distance")
savefig("center_of_mass.pdf")
println(" Created figure: center_of_mass.pdf")


text:
module SampleSimpleFail

using Test

test() = @test false

end # module


text:
#!/usr/bin/env julia
import TOML

toml = Dict()
toml["interval"] = 60
toml["output_directory"] = "/www-root"
toml["title"] = "status.julialang.org"
toml["desc"] = "Status page for services hosted by JuliaLang"
toml["influxdb"] = Dict(
    "server" => "http://influxdb:8086",
    "db" => "telegraf",
    "measurement" => "http_response",
)

services = []
for pkgserver in [
 "https://eu-central.pkg.julialang.org",
 "https://us-west.pkg.julialang.org",
 "https://us-east.pkg.julialang.org",
 "https://us-east2.pkg.julialang.org",
 "https://us-east-ci.pkg.julialang.org",
 "https://kr.pkg.julialang.org",
 "https://sg.pkg.julialang.org",
 "https://in.pkg.julialang.org",
 "https://au.pkg.julialang.org",
 "https://cn-southeast.pkg.julialang.org",
 "https://cn-east.pkg.julialang.org",
 "https://cn-northeast.pkg.julialang.org",
 ]
    name = String(match(r"https://(.*)", pkgserver)[1])
    # checks = [pkgserver * "/meta", pkgserver * "/registries"]
    checks = [pkgserver * "/meta"]
    group = "Package servers"
    push!(services, Dict("name"=>name, "checks"=>checks, "group"=>group))
end

for sserver in ["https://us-east.storage.juliahub.com", "https://kr.storage.juliahub.com"]
    name = String(match(r"https://(.*)", sserver)[1])
    checks = [sserver * "/meta"]
    group = "Storage servers"
    push!(services, Dict("name"=>name, "checks"=>checks, "group"=>group))
end
toml["services"] = services

open(joinpath(@__DIR__, "config.toml"), "w") do io
    TOML.print(io, toml)
end


text:
module TestModule

using Dates
using Mongoc: BSONObjectId

struct ChildType
    c1::String
    c2::Int
    c3::UInt8
    c4::Int8
    c5::UInt16
    c6::Int16
    c7::DateTime
    c8::Date
    c9::Vector{Date}
    c10::UInt32
end

function Base.:(==)(c1::ChildType, c2::ChildType)
    return (c1.c1 == c2.c1
            && c1.c2 == c2.c2
            && c1.c3 == c2.c3
            && c1.c4 == c2.c4
            && c1.c5 == c2.c5
            && c1.c6 == c2.c6
            && c1.c7 == c2.c7
            && c1.c8 == c2.c8
            && c1.c9 == c2.c9
            && c1.c10 == c2.c10)
end

struct FatherType
    f1::String
    f2::ChildType
    f3::Vector{Int}
    f4::BSONObjectId
    f5::Symbol
    f6::Char
end

function Base.:(==)(f1::FatherType, f2::FatherType)
    return (f1.f1 == f2.f1
            && f1.f2 == f2.f2
            && f1.f3 == f2.f3
            && f1.f4 == f2.f4
            && f1.f5 == f2.f5
            && f1.f6 == f2.f6)
end

struct Periods
    y::Year
    m::Month
    d::Day
    hh::Hour
    mm::Minute
end

struct ManyDicts
    d1::Dict{Symbol, Date}
    d2::Dict{Int, Year}
    d3::Dict{Int, Periods}
end

function Base.:(==)(m1::ManyDicts, m2::ManyDicts)
    return m1.d1 == m2.d1 && m1.d2 == m2.d2 && m1.d3 == m2.d3
end

struct SingletonStruct
end

struct OptionNothing
    val::Union{Nothing, Int}
end

struct OptionMissing
    val::Union{Missing, String}
end

struct DateEncodedAsString
    date::Date
end

module Submodule
    struct SubStruct
        val::Int
    end
end

struct StructFloat
    val::Float64
end

struct StructUInt64
    val::UInt64
end

struct StructAbsTypes
    val::Any
end

struct UnionPeriods
    val::Union{Dates.Day, Dates.Month, Dates.Year}
end

function Base.:(==)(s1::StructAbsTypes, s2::StructAbsTypes)
    return s1.val == s2.val
end

struct Lift2
    op::Function
    arg1::Number
    arg2::Number
end

exec_op(l::Lift2) = l.op(l.arg1, l.arg2)

end # TestModule


text:
abstract type AbstractSolver end

function optimize(solver::AbstractSolver, cache, x0, lx, ux, lg, ug, rows, cols)
    xopt = [0.0, 0.0]
    fopt = 0.0
    info = "example"
    out = "other data"
    return xopt, fopt, info, out
end

"""
    Options(;sparsity=DensePattern(), derivatives=ForwardFD(), solver=IPOPT())

Options for SNOW.  Default is dense, forward finite differencing, and IPOPT.

# Arguments
- `sparsity::AbstractSparsityPattern`: specify sparsity pattern
- `derivatives::AbstractDiffMethod`: specific differentiation methods to use
    or `derivatives::Vector{AbstractDiffMethod}`: vector of length two,
    first for gradient differentiation method, second for jacobian differentiation method
- `solver::AbstractSolver`: specificy which optimizer to use
"""
struct Options{T1,T2,T3}
    sparsity::T1  # AbstractSparsityPattern
    derivatives::T2  # AbstractDiffMethod
    solver::T3  # AbstractSolver
end

# defaults
Options(; sparsity=DensePattern(), derivatives=ForwardFD(), solver=IPOPT()
    ) = Options(sparsity, derivatives, solver)

# private convenience method to size bounds of length 1 to required length
function resizebounds(lx, nx)
    if length(lx) == 1 && nx > 1
        lx = lx*ones(nx)
    end
    return lx
end

"""
    minimize(func!, x0, ng, lx=-Inf, ux=Inf, lg=-Inf, ug=0.0, options=Options())

solve the optimization problem

``\\text{minimize} \\quad    f(x) \\\\``
``\\text{subject to} \\quad  l_g \\le g(x) \\le u_g \\\\``
``\\quad\\quad\\quad\\quad   l_x \\le x \\le u_x``

`f = func!(g, x)`, unless user-supplied derivatives then: `f = func!(g, df, dg, x)`

equality constraint for the ith constraint: `lg[i] = ug[i]`
"""
function minimize(func!, x0, ng, lx=-Inf, ux=Inf, lg=-Inf, ug=0.0, options=Options())

    # initialize
    nx = length(x0)

    # parse bounds
    lx = resizebounds(lx, nx)
    ux = resizebounds(ux, nx)
    lg = resizebounds(lg, ng)
    ug = resizebounds(ug, ng)

    # create cache
    cache = createcache(options.sparsity, options.derivatives, func!, nx, ng)

    # determine sparsity pattern
    rows, cols = getsparsity(options.sparsity, nx, ng)

    # returns xopt, fopt, info, out
    return optimize(options.solver, cache, x0, lx, ux, lg, ug, rows, cols)
end

# -------- Alternate Interface ----------

# mutable struct Problem{TFUN, TF, TS, TI}
#     func::TFUN
#     name::TS
#     x0::Vector{TF}
#     lx::Vector{TF}
#     ux::Vector{TF}
#     xlen::Vector{TI}
#     xnames::Vector{TS}
#     lg::Vector{TF}
#     ug::Vector{TF}
#     glen::Vector{TI}
#     fnames::Vector{TS}
# end

# function createproblem(func!, name="OptProb")
#     return Problem(func!, name, Float64[], Float64[], Float64[], Int[], String[], Float64[], Float64[], Int[], String["obj"])
# end

# function add_dv!(prob, x, lx=-Inf, ux=Inf, names=[])
#     nx = length(x)
#     lx = resizebounds(lx, nx)
#     ux = resizebounds(ux, nx)

#     append!(prob.x0, x)
#     append!(prob.lx, lx)
#     append!(prob.ux, ux)
#     append!(prob.xlen, nx)
#     if length(names) == nx
#         append!(prob.xnames, names)
#     end
# end

# function add_con!(prob, ng, lg=-Inf, ug=0.0, names=[])
#     lg = resizebounds(lg, ng)
#     ug = resizebounds(ug, ng)

#     append!(prob.lg, lg)
#     append!(prob.ug, ug)
#     append!(prob.glen, ng)
#     if length(names) == ng
#         append!(prob.fnames, names)
#     end
# end

# function minimize(prob, options)

#     ng = length(prob.lg)

#     # rewrite names if applicable
#     if options.solver isa Snopt && length(prob.xnames) == nx && length(prob.fnames) == 1+ng
#         names = Snopt.Names(prob.name, prob.xnames, prob.fnames)
#         solver = Snopt(options.solver.options, names, options.solver.warmstart)
#         options = Options(options.sparsity, options.derivatives, solver)
#     end
#     return minimize(prob.func, prob.x0, ng, lx=prob.lx, ux=prob.ux, lg=prob.lg, ug=prob.ug, options=options)
# end


text:
using Elefridge
using ZfpCompression
using PyPlot
using PyCall
using ColorSchemes
using StatsBase, Statistics
xr = pyimport("xarray")
ccrs = pyimport("cartopy.crs")

path = "/Users/milan/cams/gridded/"
filelist = filter(x->endswith(x,"_go3.grib"),readdir(path))
grib = xr.open_dataarray(joinpath(path,filelist[end]),engine="cfgrib")
X = grib.data
lat = grib.latitude.data
lon = grib.longitude.data

level = 85
X = permutedims(X,[3,2,1])      # for lossless longitude first
o3 = X[:,:,level]

## compression round
o3r = round(o3,5)

## compression zfp
X = permutedims(X,[3,2,1])      # for zfp vert x lat x lon
o3z = zfp_decompress(zfp_compress(X,precision=10))[level,:,:]'

## PLOT
lat_div(n::Integer) = Array(-90:180/(n-1):90)
lon_div(n::Integer) = Array(0:360/n:360)
lond = lon_div(6)

# pygui(true)
# ion()
ioff()
cmap = ColorMap(ColorSchemes.turku.colors)
vmin,vmax = extrema(o3)
fig,(ax0,ax1,ax2) = plt.subplots(3,1,figsize=(7,10),sharex=true,sharey=true,
                    subplot_kw=Dict("projection"=>ccrs.Robinson(central_longitude=180)))

ax2.pcolormesh(lon,lat,o3',transform=ccrs.PlateCarree();vmin,vmax,cmap)
ax1.pcolormesh(lon,lat,o3',transform=ccrs.PlateCarree();vmin,vmax,cmap)
ax0.pcolormesh(lon,lat,o3',transform=ccrs.PlateCarree();vmin,vmax,cmap)

# ax0.coastlines(color="k",lw=0.5)
# ax1.coastlines(color="k",lw=0.5)
#
# for longi in lond[2:end-1]
#     ax0.plot(repeat([longi,],31),lat_div(31),"w",transform=ccrs.PlateCarree(),lw=0.25)
#     ax1.plot(repeat([longi,],31),lat_div(31),"w",transform=ccrs.PlateCarree(),lw=0.25)
# end

ax0.text(0.9,0.95,"a",fontweight="bold",transform=ax0.transAxes)
ax1.text(0.9,0.95,"b",fontweight="bold",transform=ax1.transAxes)
ax2.text(0.9,0.95,"c",fontweight="bold",transform=ax2.transAxes)

tight_layout(rect=[0,0,1,1])
savefig("/Users/milan/git/Elefridge.jl/maps/turing2_o3_$level.png",dpi=200)
close(fig)