text:
using Gen
using Statistics

import Base.+
import Base.*
import Base.==

#### Type defs

abstract type UnrealizedValue end
abstract type UnrealizedNormalValue <: UnrealizedValue end

abstract type UnrealizedNormalLeafValue <: UnrealizedNormalValue end

mutable struct UnrealizedNormalIndyLeafValue <: UnrealizedNormalLeafValue
    #Use only when all instances are posterior independent by construction. This means committing a priori to which observations will be made.
    baseMean::Float64
    currentCondSD::Float64
    currentMean::Float64
    resetCount::Int64
end

UnrealizedNormalIndyLeafValue(μ, σ) = UnrealizedNormalIndyLeafValue(μ, σ, μ)

UnrealizedNormalIndyLeafValue(μ, σ, μt) = UnrealizedNormalIndyLeafValue(μ, σ, μt, 0)

function copyTo!(source::UnrealizedNormalIndyLeafValue, dest::UnrealizedNormalIndyLeafValue)
    dest.baseMean = source.baseMean
    dest.currentCondSD = source.currentCondSD
    dest.currentMean = source.currentMean
    dest.resetCount = dest.resetCount + 1
end

mutable struct UnrealizedNormalNodeValue <: UnrealizedNormalValue
    offset::Float64
    components::Vector{Tuple{Float64,UnrealizedNormalLeafValue}} #should be unique
end

mutable struct ObservedNormal
    unrealized::UnrealizedNormalNodeValue
    SD::Float64
end

struct UnrealizedDistribution{T} <: Gen.Distribution{T}
    func
end

##### commute definitions (so we only have to define once)

Base.:+(x::Int64,y::UnrealizedNormalValue) = y+x #definitions below should put the "more abstract" thing 1st.
Base.:+(x::UnrealizedNormalLeafValue,y::UnrealizedNormalNodeValue) = y+x #definitions below should put the "more abstract" thing 1st.
Base.:*(x::Int64,y::UnrealizedNormalValue) = y*x #definitions below should put the "more abstract" thing 1st.
Base.:-(x::UnrealizedValue,y::UnrealizedValue) = x + y * -1. #subtraction as addition

##### define ==

function ==(x::UnrealizedNormalIndyLeafValue,y::UnrealizedNormalIndyLeafValue)
    x.baseMean == y.baseMean #FIXME: Incomplete.
end

##### define +

Base.:+(x::UnrealizedNormalLeafValue,y::Float64) = UnrealizedNormalNodeValue(y,[(1.,x)])

Base.:+(x::UnrealizedNormalNodeValue,y::Float64) = UnrealizedNormalNodeValue(x.offset + y,x.components)

Base.:+(x::UnrealizedNormalLeafValue,y::UnrealizedNormalLeafValue) = UnrealizedNormalNodeValue(0.,[(1.,x),(1.,y)])

function Base.:+(anode::UnrealizedNormalNodeValue,aleaf::UnrealizedNormalLeafValue, amul = 1.)
    c = copy(anode.components)
    ret = UnrealizedNormalNodeValue(anode.offset,c)
    for (i, (nmul, nleaf)) in enumerate(c)
        if aleaf === nleaf
            #TODO: remove if nmul + amul == 0.
            c[i] = (nmul + amul, aleaf)
            return ret
        end
    end
    push!(c, (amul,aleaf))
    ret
end

function Base.:+(anode::UnrealizedNormalNodeValue,bnode::UnrealizedNormalNodeValue)
    ret = anode
    for (bmul, bleaf) in bnode.components
        ret = +(ret, bleaf, bmul)
    end
    ret.offset = anode.offset + bnode.offset
    ret
end

##### define *

function Base.:*(x::UnrealizedNormalLeafValue,y::Float64)
    if y == 0.
        UnrealizedNormalNodeValue(0.,[])
    end
    UnrealizedNormalNodeValue(0.,[(y,x)])
end

function Base.:*(anode::UnrealizedNormalNodeValue,y::Float64)
    if y == 0.
        @debug "Why multiply by zero?"
    end
    ret = UnrealizedNormalNodeValue(anode.offset * y,copy(anode.components))
    for (i,(nmul,nleaf)) in enumerate(ret.components)
        ret.components[i] = (nmul * y, nleaf)
        if nmul == 0.
            @debug "Why nmul == zero?"
        end
    end
    ret
end

##### utility functions

function meanvarmean(node::UnrealizedNormalNodeValue)
    (μ0, σ2, μt) = node.offset, 0., node.offset
    for (mul,l) in node.components
        μ0 += l.baseMean * mul
        μt += l.currentMean * mul
        σ2 += (l.currentCondSD * mul)^2
        if isnan(μt)
            @debug "bad leaf?" l node
        end
    end
    (μ0, σ2, μt)
end

function meanvarmean(o::ObservedNormal)
    (μ0, σ2, μt) = meanvarmean(o.unrealized)
    σ2 += o.SD^2
    (μ0, σ2, μt)
end

##### define as if distributions. There must be a better way.

#makeUnrealizedIndyNormal(μ::Float64, σ::Float64) = UnrealizedNormalIndyLeafValue(μ, σ)
unrealizedIndyNormal = UnrealizedDistribution{UnrealizedNormalIndyLeafValue}(UnrealizedNormalIndyLeafValue)

ObservedNormal(μ::UnrealizedNormalLeafValue,σ::Float64) = ObservedNormal(μ + 0., σ)
observedNormal = UnrealizedDistribution{Float64}(ObservedNormal)

##### implement distribution behavior

function sampleFrom!(unilv::UnrealizedNormalIndyLeafValue)
    #@debug "sampleFrom 1"
    if unilv.currentCondSD > 0.
        val = normal(unilv.currentMean, unilv.currentCondSD)
        unilv.currentMean = val
        unilv.currentCondSD = 0.
    else
        val = unilv.currentMean
    end
    val
end

function sampleFrom!(unnv::UnrealizedNormalNodeValue)
    #@debug "sampleFrom 2"
    val = unnv.offset
    for (mul, leaf) in unnv.components
        val += mul * sampleFrom!(leaf)
    end
    val
end

function sampleFrom!(obsn::ObservedNormal)
    #@debug "sampleFrom 3"
    sampleFrom!(obsn.unrealized) + normal(0.,obsn.SD)
end

Gen.random(d::UnrealizedDistribution, args...) = d.func(args...) #This is the whole point - instead of sampling, just record

function Gen.logpdf(d::UnrealizedDistribution{UnrealizedNormalIndyLeafValue}, v::UnrealizedValue, args...)
    #@debug "logpdf UnrealizedNormalIndyLeafValue; resetting" v.resetCount+1
    baseVal = d.func(args...)
    copyTo!(baseVal, v)

    return 0.
    # if (d.func(args...) == v)
    #     return 0.
    # else
    #     @debug "unrealized value must equal itself" d.func(args...) v
    # end
end

function Gen.logpdf_grad(d::UnrealizedDistribution{UnrealizedNormalIndyLeafValue}, v::UnrealizedValue, args...)

    return 0.
    # if (d.func(args...) == v)
    #     return 0.
    # else
    #     @debug "logpdf_grad: unrealized value must equal itself" d.func(args...) v
    # end
end
(d::UnrealizedValue)(args...) = Gen.random(d, args...)
Gen.has_output_grad(d::UnrealizedValue) = true
Gen.has_argument_grads(d::DistributionsBacked) = false

function Gen.random(wrapped::UnrealizedDistribution{Float64}, args...)
    #@debug "random outer"
    d = wrapped.func(args...)
    sampleFrom!(d)
end

(d::UnrealizedDistribution{Float64})(args...) = Gen.random(d, args...)

Gen.logpdf_grad(d::UnrealizedDistribution{Float64}, v::UnrealizedValue, args...) = throw("unimplemented")
Gen.has_output_grad(d::UnrealizedDistribution{Float64}) = false
Gen.has_argument_grads(d::UnrealizedDistribution{Float64}) = false

function Gen.logpdf(wrapped::UnrealizedDistribution{Float64}, v::Float64, args...) #Observe the value. Side effects! This is a hack, for now!
    #@debug "logpdf outer"
    d = wrapped.func(args...)
    Gen.logpdf(d, v)
end

function Gen.logpdf(d::ObservedNormal, v::Float64)
    μ0, σ2, μt = meanvarmean(d)
    lpdf = Gen.logpdf(normal, v, μt, sqrt(σ2))

    #@debug "logpdf" d μ0 σ2 μt v lpdf
    rawDelta = v - μ0
    un = d.unrealized
    com = un.components

    #update all independent leafs
    for (i, (mul, leaf)) in enumerate(com)
        if leaf.currentCondSD > 0. && mul != 0.
            delta = rawDelta + mul * (leaf.baseMean - leaf.currentMean) #use unconditional mean for other components, as they're independent
            likePrec = 1/d.SD^2 #likelihood precision
            totalPrec = likePrec + 1/(mul * leaf.currentCondSD)^2 #ignore variance from other (inpependent) components; that is, condition on point values
            cm, cσ = leaf.currentMean, leaf.currentCondSD
            leaf.currentMean += delta * likePrec / totalPrec / mul
            leaf.currentCondSD = sqrt(1 / totalPrec / mul^2)
            if isnan(leaf.currentMean) || leaf.currentCondSD == 0
                #@debug "updating..." μ0 σ2 μt v delta cm cσ delta/d.SD d.SD likePrec mul totalPrec leaf.currentMean leaf.currentCondSD
            end
        end
    end
    lpdf
end

function compareModels(integratedModel, explicitModel, args, constraints, locus, meansteps = meansteps)
    tr, w = generate(integratedModel, args, constraints)
    tr2s, logws, w2 = importance_sampling(explicitModel, args, constraints, meansteps)
    ws = AnalyticWeights(exp.(logws[1:meansteps]))
    locusvals = [tr2s[i][locus] for i in 1:meansteps]
    m = mean(locusvals, ws)
    vv = var(locusvals, ws)
    @warn "Comparing models:" w w2 m tr[locus].currentMean vv tr[locus].currentCondSD^2
end

max_stars_repo_name: jamesonquinn/Gen

max_issues_repo_name: jamesonquinn/Gen

max_forks_repo_name: jamesonquinn/Gen


text:
import Distributions: logpdf, pdf

struct SDT{T1,T2} <: ContinuousUnivariateDistribution
    d::T1
    c::T2
end

logpdf(d::SDT, data::Vector{Int64}) = logpdf(d, data...)

logpdf(d::SDT, data::Tuple{Vararg{Int64}}) = logpdf(d, data...)

function logpdf(d::SDT, hits, fas, Nd)
    @unpack d,c = d
    θhit = cdf(Normal(0, 1), d/2-c)
    θfa = cdf(Normal(0, 1), -d/2-c)
    loghits = logpdf(Binomial(Nd, θhit), hits)
    logfas = logpdf(Binomial(Nd, θfa), fas)
    return loghits + logfas
end

pdf(d::SDT, data::Vector{Int64}) = exp(logpdf(d, data...))

ext: jl


text:
using Compat, ClassicalLaminateTheory

forcesAndMoments =
 @compat Dict(
    :Nx => 1000,
    :Ny => 300,
    :Nxy => 150,
    :Mx => 2.5,
    :My => 5,
    :Mxy => 10)

materialProperties =@compat Dict(
  :E1 => [20.0e6, 13.0*10^6],
  :E2 => [1.5*10^6, 1.0*10^6],
  :G12 => [1.0*10^6, 0.75*10^6],
  :ν12 => [0.34, 0.34],
  :alpha1 => [0.2e-6, 0.2e-5],
  :alpha2 => [0.2e-4, 0.2e-3],
  :beta1 => [0.2e-6, 0.2e-5],
  :beta2 => [0.2e-4, 0.2e-3],
  :thickness => [0.0075, 0.0055])

laminateProperties =@compat Dict(
  :nplies => 4,
  :orientation => [0, 90, 0, 90],
  :materials => [2, 2, 1, 1],
  :repeats => 2,
  :symmetric => true,
  :symmetricrepeats => 1,
  :deltaTemp => -280.0,
  :deltaMoisture => 2
)

materialProperties |> display
println()
createLaminate!(laminateProperties, materialProperties, forcesAndMoments)
show(laminateProperties)

println("ABD matrix:")
round(laminateProperties[:ABD], 3) |> display

println()
println("Total forces and moments:")
round(laminateProperties[:fTotal], 7) |> display


text:
export transform

"""
    transform(problem::InitialValueProblem, options::Options)

Interface function that calls the respective transformation function.

### Input

- `problem` -- discrete or continuous initial-value problem
- `option`  -- problem options

### Output

A tuple containing the transformed problem and the transformed options.

### Notes

The functions that are called in the background should return a the transformed
system components `A`, `X0`, and `U`, and also an inverse transformation matrix `M`.
If the system has an invariant, it is transformed as well.
"""
function transform(problem::InitialValueProblem, options::Options)
    method = options[:coordinate_transformation]
    if method == ""
        nothing # no-op
    elseif method == "schur"
        options[:transformation_matrix] = Tm
        problem = schur_transform(problem)
    else
        error("the transformation method $method is undefined")
    end
    return (problem, options)
end

"""
    schur_transform(problem::InitialValueProblem)

Applies a Schur transformation to a discrete or continuous initial-value problem.

### Input

- `S` -- discrete or continuous initial-value problem

### Output

Transformed problem.

### Algorithm

We use Julia's default `schurfact` function to compute a
[Schur decomposition](https://en.wikipedia.org/wiki/Schur_decomposition)
of the coefficients matrix ``A``.
"""
function schur_transform(problem::InitialValueProblem{PT, ST}
                         ) where {PT <: Union{ConstrainedLinearControlDiscreteSystem, ConstrainedLinearControlContinuousSystem}, ST<:LazySet}

    n = size(problem.s.A, 1)

    # full (dense) matrix is required by schur
    A_new, T_new = schur(Matrix(problem.s.A))

    # recall that for Schur matrices, inv(T) == T'
    Z_inverse = copy(transpose(T_new))

    # apply transformation to the initial states
    X0_new = Z_inverse * problem.x0

    # apply transformation to the inputs
    U_new = Z_inverse * problem.s.U

    # obtain the transformation matrix for reverting the transformation again
    T_inverse = F[:vectors]

    # apply transformation to the invariant
    if hasmethod(stateset, Tuple{typeof(problem.s)})
        invariant_new = T_inverse * problem.s.X
    else
        invariant_new = Universe(n)
    end

    system_new = (A_new, I(n), invariant_new, U_new)
    return InitialValueProblem(PT(system_new), X0_new)
end


text:
struct FDGZProblem{TBCL, TBCR}
    a # (a,b) define the domain
    b
    Δx # mesh spacing
    n # number of points in the interior, excluding a and b
    x # interior mesh
    ϵ # scalar parameter
    F # holds the residual values
    left_bc::TBCL
    right_bc::TBCR
end

struct DirichletBC
    u_bc
end

struct NeumannBC
    v_bc
end

struct PeriodicBC
end


text:
gear = [KPDSupply("map", 9, 150),
        KPDSupply("compass", 13, 35),
        KPDSupply("water", 153, 200),
        KPDSupply("sandwich", 50, 160),
        KPDSupply("glucose", 15, 60),
        KPDSupply("tin", 68, 45),
        KPDSupply("banana", 27, 60),
        KPDSupply("apple", 39, 40),
        KPDSupply("cheese", 23, 30),
        KPDSupply("beer", 52, 10),
        KPDSupply("suntan cream", 11, 70),
        KPDSupply("camera", 32, 30),
        KPDSupply("T-shirt", 24, 15),
        KPDSupply("trousers", 48, 10),
        KPDSupply("umbrella", 73, 40),
        KPDSupply("waterproof trousers", 42, 70),
        KPDSupply("waterproof overclothes", 43, 75),
        KPDSupply("note-case", 22, 80),
        KPDSupply("sunglasses", 7, 20),
        KPDSupply("towel", 18, 12),
        KPDSupply("socks", 4, 50),
        KPDSupply("book", 30, 10)]

pack = solve(gear, 400)
println("The hicker should pack: \n - ", join(pack, "\n - "))
println("\nPacked weight: ", mapreduce(x -> x.weight, +, pack), " kg")
println("Packed value: ", mapreduce(x -> x.value, +, pack), " €")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module GeometryBasics

using StaticArrays, Tables, StructArrays, IterTools, LinearAlgebra
using EarCut_jll

using Base: @propagate_inbounds

include("fixed_arrays.jl")
include("offsetintegers.jl")
include("basic_types.jl")

include("primitives/rectangles.jl")
include("primitives/spheres.jl")
include("primitives/cylinders.jl")
include("primitives/pyramids.jl")
include("primitives/particles.jl")

include("interfaces.jl")
include("metadata.jl")
include("viewtypes.jl")
include("geometry_primitives.jl")
include("meshes.jl")
include("triangulation.jl")
include("lines.jl")
include("boundingboxes.jl")

export AbstractGeometry, GeometryPrimitive
export Mat, Point, Vec
export LineFace, Polytope, Line, NgonFace, convert_simplex
export LineString, AbstractPolygon, Polygon, MultiPoint, MultiLineString, MultiPolygon
export Simplex, connect, Triangle, NSimplex, Tetrahedron
export QuadFace, metafree, coordinates, TetrahedronFace
export TupleView, SimplexFace, Mesh, meta
export Triangle, TriangleP
export AbstractFace, TriangleFace, QuadFace, GLTriangleFace
export OffsetInteger, ZeroIndex, OneIndex, GLIndex
export FaceView, SimpleFaceView
export AbstractPoint, PointMeta, PointWithUV
export PolygonMeta, MultiPointMeta, MultiLineStringMeta, MeshMeta, LineStringMeta,
       MultiPolygonMeta
export decompose, coordinates, faces, normals, decompose_uv, decompose_normals,
       texturecoordinates
export Tesselation, pointmeta, Normal, UV, UVW
export GLTriangleFace, GLNormalMesh3D, GLPlainTriangleMesh, GLUVMesh3D, GLUVNormalMesh3D
export AbstractMesh, Mesh, TriangleMesh
export GLNormalMesh2D, PlainTriangleMesh
export MetaT, meta_table

# all the different predefined mesh types
# Note: meshes can contain arbitrary meta information,
export AbstractMesh, TriangleMesh, PlainMesh, GLPlainMesh, GLPlainMesh2D, GLPlainMesh3D
export UVMesh, GLUVMesh, GLUVMesh2D, GLUVMesh3D
export NormalMesh, GLNormalMesh, GLNormalMesh2D, GLNormalMesh3D
export NormalUVMesh, GLNormalUVMesh, GLNormalUVMesh2D, GLNormalUVMesh3D
export NormalUVWMesh, GLNormalUVWMesh, GLNormalUVWMesh2D, GLNormalUVWMesh3D

# mesh creation functions
export triangle_mesh, triangle_mesh, uv_mesh
export uv_mesh, normal_mesh, uv_normal_mesh

export height, origin, radius, width, widths, xwidth, yheight
export HyperSphere, Circle, Sphere
export Cylinder, Cylinder2, Cylinder3, Pyramid, extremity
export Rect, Rect2D, Rect3D, IRect, IRect2D, IRect3D, FRect, FRect2D, FRect3D
export before, during, isinside, isoutside, meets, overlaps, intersects, finishes
export centered, direction, area, update
export max_dist_dim, max_euclidean, max_euclideansq, min_dist_dim, min_euclidean
export min_euclideansq, minmax_dist_dim, minmax_euclidean, minmax_euclideansq
export self_intersections, split_intersections

if Base.VERSION >= v"1.4.2"
    include("precompile.jl")
    _precompile_()
end

end # module


text:
import PWS
using PWS: SMCEstimate, DirectMCEstimate, TIEstimate, marginal_configuration, MarginalEnsemble, gene_expression_system, generate_configuration, logpdf

system_fn = () -> PWS.gene_expression_system(dtimes=0:0.1:10)
smc = SMCEstimate(256)
dmc = DirectMCEstimate(8*256)
ti  = TIEstimate(0, 8, 256)

system = system_fn()

conf = generate_configuration(system)

sresult = PWS.mutual_information(system, smc, num_responses=200)
tresult = PWS.mutual_information(system, ti, num_responses=5)
dresult= PWS.mutual_information(system, dmc, num_responses=5)

using Plots
plot(conf.s_traj, xlim=(0,2), seriescolor=:green, label="S")
plot!(conf.x_traj, seriescolor=:cornflowerblue, label="X")
savefig("~/Downloads/example_traj.pdf")

plot(system.dtimes, result.MutualInformation, color=:gray, legend=false)
plot!(system.dtimes, dresult.MutualInformation, color=:pink, label="")

using Statistics
plot(system.dtimes, mean(result.MutualInformation), ribbon=sqrt.(var(result.MutualInformation)./size(result, 1)), label="SMC", ylabel="Path mutual information", xlabel="Trajectory length")
plot!(system.dtimes, mean(dresult.MutualInformation), ribbon=sqrt.(var(dresult.MutualInformation)./size(dresult, 1)), label="DMC")
plot!(system.dtimes, mean(tresult.MutualInformation), ribbon=sqrt.(var(tresult.MutualInformation)./size(tresult, 1)), label="TI")

using CSVFiles, DrWatson, DataFrames

zechner_res = DataFrame(load(datadir("zechner/gene_exp.csv")))
plot!(zechner_res.Duration[1:20], zechner_res.PMI[1:20])

savefig(plot!(dpi=100, size=(6*72,3.5*72)), "~/Downloads/gene_expression.pdf")

ens = MarginalEnsemble(system)
r = map(1:50) do i
    smc = SMCEstimate(100)
    a = PWS.log_marginal(PWS.simulate(smc, conf, ens))[end]

    dmc = DirectMCEstimate(100)
    b = PWS.log_marginal(PWS.simulate(dmc, conf, ens))[end]
    (a, b)
end

r

scatter(randn(50)*0.01, cond .- getindex.(r, 1))
scatter!(randn(50)*0.01 .+ 1, cond .- getindex.(r, 2))

using Statistics
mean(getindex.(r, 1))
mean(getindex.(r, 2))

smc = SMCEstimate(10000)
PWS.log_marginal(PWS.simulate(smc, conf, ens))[end]

dmc = DirectMCEstimate(10000)
PWS.log_marginal(PWS.simulate(dmc, conf, ens))[end]

cond = -PWS.energy_difference(conf, ens)

histogram(PWS.simulate(dmc, conf, ens).samples[end,:])

r = map(1:1000) do i
    result = PWS.sample(conf, ens)
    PWS.energy_difference(result, ens) => result
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Implementation of HdrImage for the manipulation and elaboration of and HDR image

"""
    HdrImage

Wrapper of a `Matrix` of elements of type `RGB{Float32}`, used to represent an image in hdr format.
"""
struct HdrImage
    pixel_matrix::Matrix{RGB{Float32}}
end

@doc """
    HdrImage(pixel_matrix::Matrix{RGB{Float32}})

Constructor for a [`HdrImage`](@ref) instance.
""" HdrImage(::Matrix{RGB{Float32}})

"""
    HdrImage(img_width::Integer, img_height::Integer)

Construct an [`HdrImage`](@ref) wrapping a zero-initialized matrix of size `(img_width, img_height)`.

# Examples

```jldoctest
julia> HdrImage(3, 2)
3x2 HdrImage:
 (0.0 0.0 0.0)  (0.0 0.0 0.0)
 (0.0 0.0 0.0)  (0.0 0.0 0.0)
 (0.0 0.0 0.0)  (0.0 0.0 0.0)
```
"""
HdrImage(img_width::Integer, img_height::Integer) = HdrImage(zeros(RGB{Float32}, img_width, img_height))

"""
    HdrImage(arr::AbstractArray{<:Any, 1}, im_width::Integer, im_height::Integer)
    HdrImage(arr::AbstractArray{<:Any, 1}, shape)

Construct an [`HdrImage`](@ref) wrapping a matrix obtained from `reshape`.

# Examples

```jldoctest
julia> arr = [RGB( 1.,  2.,  3.), RGB( 4.,  5.,  6.), RGB( 7.,  8.,  9.),
              RGB(10., 11., 12.), RGB(13., 14., 15.), RGB(16., 17., 18.)];

julia> HdrImage(arr, 3, 2)
3x2 HdrImage:
 (1.0 2.0 3.0)  (10.0 11.0 12.0)
 (4.0 5.0 6.0)  (13.0 14.0 15.0)
 (7.0 8.0 9.0)  (16.0 17.0 18.0)
```
"""
HdrImage(arr::AbstractArray{<:Any, 1}, im_width::Integer, im_height::Integer) = HdrImage(reshape(arr, im_width, im_height))
HdrImage(arr::AbstractArray{<:Any, 1}, shape) = HdrImage(reshape(arr, shape))

#############
# Iterations

length(image::HdrImage) = length(image.pixel_matrix)

firstindex(image::HdrImage) = firstindex(image.pixel_matrix)
firstindex(image::HdrImage, d) = firstindex(image.pixel_matrix, d)

lastindex(image::HdrImage) = lastindex(image.pixel_matrix)
lastindex(image::HdrImage, d) = lastindex(image.pixel_matrix, d)

getindex(image::HdrImage, inds...) = getindex(image.pixel_matrix, inds...)

setindex!(image::HdrImage, value, key) = setindex!(image.pixel_matrix, value, key)
setindex!(image::HdrImage, X, inds...) = setindex!(image.pixel_matrix, X, inds...)

iterate(image::HdrImage, state = 1) = state > lastindex(image) ? nothing : (image[state], state + 1)

###############
# Broadcasting

axes(image::HdrImage) = axes(image.pixel_matrix)
axes(image::HdrImage, d) = axes(image.pixel_matrix, d)

broadcastable(image::HdrImage) = image
broadcastable(::Type{<:HdrImage}) = HdrImage

BroadcastStyle(::Type{<:HdrImage}) = Style{HdrImage}()
BroadcastStyle(::Style{HdrImage}, ::BroadcastStyle) = Style{HdrImage}()

copy(bc::Broadcasted{Style{HdrImage}}) = copy(convert(Broadcasted{Broadcast.DefaultArrayStyle{RGB{Float32}}}, bc))

###############
# Tone mapping

@doc raw"""
    luminosity(image::HdrImage; δ::Float32 = eps(Float32))

Return the average luminosity an [`HdrImage`](@ref) as the logaritmic mean of the [`luminosity(::RGB)`](@ref)
``l_i`` of each pixel:

```math
\left< l \right> = 10^{\frac{\sum_i \log_{10}(\delta + l_i)}{N}}
```

The parameter `δ` avoid singularities for ``l_i = 0`` (black pixels).

# Examples

```jldoctest
julia> arr = [RGB( 1.,  2.,  3.), RGB( 4.,  5.,  6.), RGB( 7.,  8.,  9.),
              RGB(10., 11., 12.), RGB(13., 14., 15.), RGB(16., 17., 18.)];

julia> luminosity(HdrImage(arr, 3, 2))
7.706255f0
```
"""
function luminosity(image::HdrImage; δ::Float32 = eps(Float32))
    10^(sum(map(x -> log10(δ + luminosity(x)), image))/length(image))
end

"""
    normalize(image::HdrImage, α::Float32
              ; luminosity::Float32 = average_luminosity(image))

Normalize a [`HdrImage`](@ref) for a given luminosity.

If the `luminosity` parameter is not specified, the image will be normalized according to
the result of [`luminosity(::HdrImage; ::Float32)`](@ref).

# Examples

```julia
julia> arr = [RGB( 1.,  2.,  3.), RGB( 4.,  5.,  6.), RGB( 7.,  8.,  9.),
              RGB(10., 11., 12.), RGB(13., 14., 15.), RGB(16., 17., 18.)];

julia> normalize(HdrImage(arr, 3, 2), 1f0)
3x2 HdrImage:
 (0.12976472 0.25952944 0.38929415)  (1.2976472 1.4274119 1.5571766)
 (0.5190589 0.6488236 0.7785883)     (1.6869414 1.8167061 1.9464709)
 (0.90835303 1.0381178 1.1678824)    (2.0762355 2.2060003 2.335765)
```
"""
function normalize(image::HdrImage, α::Float32; luminosity::Float32 = luminosity(image))
    HdrImage([α / luminosity * pix for pix ∈ image], size(image))
end

"""
    clamp(image::HdrImage)

Adjust the color levels of the brightest pixels of a [`HdrImage`](@ref), by applying the [`clamp(::RGB)`](@ref)
function to each pixel.

# Examples

```julia
julia> arr = [RGB( 1.,  2.,  3.), RGB( 4.,  5.,  6.), RGB( 7.,  8.,  9.),
              RGB(10., 11., 12.), RGB(13., 14., 15.), RGB(16., 17., 18.)];

julia> clamp(HdrImage(arr, 3, 2))
3x2 HdrImage:
 (0.5 0.6666667 0.75)        (0.90909094 0.9166667 0.9230769)
 (0.8 0.8333333 0.85714287)  (0.9285714 0.93333334 0.9375)
 (0.875 0.8888889 0.9)       (0.9411765 0.9444444 0.94736844)
```
"""
clamp(image::HdrImage) = HdrImage(clamp.(image))

"""
    γ_correction(image::HdrImage, γ::Float32)

Compute the γ correction of a [`HdrImage`](@ref), by applying the [`γ_correction(::RGB, ::Float32)`](@ref)
function to each pixel.

Before calling this function, you should apply a tone-mapping algorithm to the image and be sure that
the R, G, and B values of the colors in the image are all in the range ``[0, 1]``. Use
[`normalize(::HdrImage, ::Float32; ::Float32)`](@ref) and [`clamp(image::HdrImage)`](@ref) to do this.

# Examples

```julia
julia> arr = [RGB( 1.,  2.,  3.), RGB( 4.,  5.,  6.), RGB( 7.,  8.,  9.),
              RGB(10., 11., 12.), RGB(13., 14., 15.), RGB(16., 17., 18.)];

julia> image = normalize(HdrImage(arr, 3, 2), 1f0) |> clamp;

julia> γ_correction(image, 1f0)
3x2 HdrImage:
 (0.11485995 0.20605269 0.28021002)  (0.5647722 0.58803856 0.6089437)
 (0.34169766 0.393507 0.43775633)    (0.6278296 0.64497536 0.660611)
 (0.47598794 0.5093512 0.53872037)   (0.67492735 0.6880849 0.7002187)

julia> γ_correction(image, 0.8f0)
3x2 HdrImage:
 (0.06686684 0.13882665 0.20387058)  (0.48960024 0.51494074 0.53792465)
 (0.26124772 0.31166682 0.35607424)  (0.558859 0.57800144 0.5955692)
 (0.39536202 0.43030033 0.4615346)   (0.6117462 0.6266897 0.64053386)

julia> γ_correction(image, 2.4f0)
3x2 HdrImage:
 (0.40588558 0.5177947 0.58855206)  (0.7881591 0.8015287 0.8132807)
 (0.63926977 0.6780008 0.7087834)   (0.82369685 0.83299613 0.8413514)
 (0.73394746 0.7549599 0.77280176)  (0.8489011 0.8557578 0.86201346)
```
"""
γ_correction(image::HdrImage, γ::Float32) = HdrImage(γ_correction.(image, γ))

################
# Miscellaneous

show(io::IO, image::HdrImage) = print_matrix(io, image.pixel_matrix)

function show(io::IO, ::MIME"text/plain", image::HdrImage)
    println(io, "$(join(map(string, size(image)), "x")) $(typeof(image)):")
    print_matrix(io, image.pixel_matrix)
end

eltype(::HdrImage) = RGB{Float32}

fill!(image::HdrImage, x) = HdrImage(fill!(image.pixel_matrix, x))

size(image::HdrImage) = size(image.pixel_matrix)


text:
using GlyphTools
using Base.Test

if VERSION.minor < 5
  macro testset(name, block)
    eval(block)
  end
end

@testset "GlyphTools" begin
  @testset "FreeType" begin
    include("freetype.jl")
  end
end


text:
using Plots
using LaTeXStrings

function arch(θ₁,θ₂;radius=1.,Δθ=1.)
  θ₁ = π*θ₁/180
  θ₂ = π*θ₂/180
  Δθ = π*Δθ/180
  l = round(Int,(θ₂-θ₁)/Δθ)
  x = zeros(l)
  y = zeros(l)
  for i in 1:l
    θ = θ₁ + i*Δθ
    x[i] = radius*cos(θ)
    y[i] = radius*sin(θ)
  end
  return x, y
end

plot()

x, y = arch(0,360)
plot(x,y,seriestype=:shape,label="",alpha=0.5)

x, y = arch(0,360,radius=0.95)
plot!(x,y,seriestype=:shape,label="",fillcolor=:white)

x, y = arch(0,360,radius=0.7)
plot!(x,y,seriestype=:shape,label="",alpha=0.5,fillcolor=:red)

x, y = arch(0,360,radius=0.65)
plot!(x,y,seriestype=:shape,label="",fillcolor=:white)

plot!([0,0],[0,1.1],arrow=true,color=:black,linewidth=2,label="")
plot!([0,1.1],[0,0],arrow=true,color=:black,linewidth=2,label="")

x, y = arch(15,16,radius=0.65)
plot!([0,x[1]],[0,y[1]],arrow=true,color=:black,linewidth=1,label="")

x, y = arch(35,36,radius=0.95)
plot!([0,x[1]],[0,y[1]],arrow=true,color=:black,linewidth=1,label="")

plot!(draw_arrow=true)
plot!(showaxis=:no,ticks=nothing,xlim=[-0.1,1.1],ylim=[-0.1,1.1],)
plot!(xlabel="x",ylabel="y",size=(400,400))

annotate!(0.58,-0.07,text(L"\Delta v_1",10))
annotate!(0.88,-0.07,text(L"\Delta v_2",10))

savefig("./velocities.pdf")


text:
"""

Assuming that we are running this module directly in a google cloud machine
with the correct permissions, there is no need to authenticate
"""
module GCSCLIProvider

using ...SimpleTasks.Types

import SimpleTasks.Services.Bucket
import SimpleTasks.Services.CLIBucket

export Details

const GCS_PREFIX = `gs:/`
const GCS_COMMAND = `gsutil`

Details() = CLIBucket.Provider(GCS_PREFIX, GCS_COMMAND)

end # module GCSCLIProvider

max_forks_repo_name: seung-lab/SimpleTasks


text:
export
   QWEvolution
"""
    QWEvolution(model, parameters, check)

Type describing standard quantum walk evolution. `check_qwdynamics` is executed if and only iff `check` is true.

 Needs implementation of
* `evolve(::QWEvolution{<:QWModelDiscr}, state)` or `evolve(::QWEvolution{<:QWModelCont}, state, time::Real)`
* `measure(::QWEvolution, state)`
* `check_qwdynamics(QWEvolution, ::QWModel, parameters::Dict{Symbol})`
* proper constructors.

Offers functions
* `execute`
* `execute_single`
* `execute_single_measured`
* `execute_all`
* `execute_all_measured`
"""
struct QWEvolution{T} <: QWDynamics{T}
  model::T
  parameters::Dict{Symbol}

  function QWEvolution(model::T,
                       parameters::Dict{Symbol};
                       check::Bool=false) where T<:QWModel
    if check
       check_qwdynamics(QWEvolution, model, parameters)
    end
    new{T}(model, parameters)
  end
end

function check_qwdynamics(qwe::QWEvolution)
  check_qwdynamics(QWEvolution, model(qwe), parameters(qwe))
end


text:
# using Revise
using ShaleDrillingLikelihood
using ShaleDrillingLikelihood.SDLParameters

# detect if using SLURM
if "SLURM_JOBID" in keys(ENV)
    SLURM_JOBID = ENV["SLURM_JOBID"]
    using ClusterManagers
else
    SLURM_JOBID = ""
end

using CountPlus, Distributed, JLD2
using Optim: minimizer, Options, BFGS, NelderMead
using ShaleDrillingLikelihood: value_function, EVobj, cost,
    theta_royalty_κ, update_kappa_level_to_cumsum!, theta_royalty
using SparseArrays: nonzeros
using Formatting: generate_formatter

num_nonzeros(x) = length(nonzeros(x))
intstr = generate_formatter("%'d")

# ------------------- number of simulations ----------------------

pargs = parse_commandline()
# using ArgParse: parse_args
# pargs = parse_args(["--noPar"], ShaleDrillingLikelihood.SDLParameters.arg_settings())
print_parsed_args(pargs)

DATAPATH = pargs["dataset"]

# estimation
M_cnstr = pargs["Mcnstr"]
M_full  = pargs["Mfull"]
DO_CNSTR = pargs["doCnstr"]
DO_FULL  = !pargs["noFull"]
MAXTIME_CNSTR = pargs["maxtimeCnstr"]
MAXTIME_FULL  = pargs["maxtimeFull"]
DO_PAR = !pargs["noPar"]

# parms
COMPUTE_INITIAL_VALUES = pargs["computeStarting"]
THETA0_FULL_OVERRIDE = pargs["theta"]
CONVERT_KAPPA = pargs["convertKappa"]

#payoffs
COST = pargs["cost"]
EXT = pargs["extension"]
REV = pargs["revenue"]
SCRAP = pargs["scrap"]
# DO_FULL = false
# REV = DrillingRevenue(Unconstrained(), TimeTrend(), GathProcess() )
# COST = DrillingCost_TimeFE(2008, 2012)
# SCRAP = ScrapValue_Constant()
# THETA0_FULL_OVERRIDE = [-12.489, -8.965, -7.696, -7.131, -6.782, 1.576, -1.495, 0.1, -2.709, 0.597, 0.34, 0.022, 0.681, 0.12, 0.594, 1.183, -1.706, 0.14, 3.798, 0.819, 1.307, 1.342, 1.073, -14.926, 0x1.8c72626bfc357p-4, 0x1.59eda12744aacp-2, ]
REWARD = DrillReward(REV, COST, EXT, SCRAP)
ANTICIPATE = pargs["anticipateT1EV"]
DISCOUNT = pargs["discount"]

# approximation
PSI = PsiSpace(pargs["numPsi"])
NUM_P = pargs["numP"]
EXTEND_GRID = pargs["extendPriceGrid"]
MINP = pargs["minTransProb"]

# print some things
println("Model is\n\t$REWARD\n")
println("Firms anticipate T1EV shocks? $ANTICIPATE")
println("Psi-space is\n\t$PSI with length $(length(PSI))")

if cost(REWARD) isa DrillingCost_TimeFE
    HASRIGS = "no-rigs"
elseif cost(REWARD) isa DrillingCost_TimeFE_rigrate
    HASRIGS = "WITH-rigs"
elseif cost(REWARD) isa DrillingCost_TimeFE_costdiffs
    HASRIGS = "COSTDIFF"
elseif cost(REWARD) isa DrillingCost_DoubleTimeFE
    HASRIGS = "no-rigs-dbl-timefe"
else
    throw(error("don't have values for this cost fct"))
end

if "SLURM_JOBID" in keys(ENV)
    DATADIR = "/home/magerton/haynesville/intermediate_data"
else
    DATADIR = "E:/projects/haynesville/intermediate_data"
end

# --------------- create data ---------------

# load in data from disk
rdatapath = joinpath(DATADIR, DATAPATH)
println_time_flush("loading $rdatapath")
data_royalty    = DataRoyalty(REWARD, rdatapath)
data_produce    = DataProduce(REWARD, rdatapath)
data_drill_prim = DataDrillPrimitive(REWARD, rdatapath)
println_time_flush("Data imported")

thetarho0 = ThetaRho()
theta0_royalty = Theta(data_royalty; θρ = thetarho0)
update_kappa_level_to_cumsum!(theta_royalty_κ(data_royalty, theta0_royalty))
theta0_produce = Theta(data_produce)
theta0_drill   = Theta(REWARD; θρ=thetarho0)

# transition matrices
zrng, ztrans = GridTransition(data_drill_prim, EXTEND_GRID, NUM_P; minp=MINP)
println_time_flush("Z range is\n\t$zrng with lengths $(length.(zrng))")
println_time_flush("Transition matrix drops probs < $MINP and has $(intstr(num_nonzeros(ztrans))) elements")

# create drilling dataset
wp = statespace(data_drill_prim)
ddm = DynamicDrillModel(REWARD, DISCOUNT, wp, zrng, ztrans, PSI, ANTICIPATE)
data_drill = DataDrill(data_drill_prim, ddm)

evdims = size(EVobj(value_function(ddm)))
println_time_flush("EV dimension is $evdims, implying $(intstr(prod(evdims))) states")

# full dataset
dataset_full = DataSetofSets(data_drill, data_royalty, data_produce, CoefLinks(data_drill))

# constrained model
rwrd_cnstr = ConstrainedProblem(REWARD, theta0_drill)
theta0_cnstr = ThetaConstrained(REWARD, theta0_drill)
data_cnstr = DataDrill(data_drill, DynamicDrillModel(ddm, rwrd_cnstr))
dataset_cnstr = DataSetofSets(data_cnstr, EmptyDataSet(), EmptyDataSet())
println_time_flush("Data created")

# ------------------- get initial values -----------------------

if COMPUTE_INITIAL_VALUES
    roy_short = DataRoyalty(RoyaltyModelNoHet(), data_royalty)
    # thet_roy = theta0_royalty[3:end]
    # kap = theta_royalty_κ(roy_short, thet_roy)
    # update_kappa_level_to_cumsum!(kap)
    # res_royalty = solve_model(roy_short, thet_roy)
    res_royalty = solve_model(roy_short, theta0_royalty[3:end])
    theta0_royalty[3:end] .= minimizer(res_royalty)

    theta0_produce .= ThetaProduceStarting(REWARD, rdatapath)
    println_time_flush("Pdxn / Royalty starting values done")
end

# ------------------- set up workers -----------------------

if DO_PAR
    pids = start_up_workers(ENV; nprocs=8)
    @everywhere using ShaleDrillingLikelihood
    println_time_flush("Library loaded on workers")
end

# Solve constrained simpler model
if DO_CNSTR
    res_c, ew_c = solve_model(dataset_cnstr, theta0_cnstr, M_cnstr, MAXTIME_CNSTR)
    theta1_cnstr = minimizer(res_c)
    updateThetaUnconstrained!(REWARD, theta0_drill, theta1_cnstr)
end

# Solve unconstrained full model
theta0s = (theta0_drill, theta0_royalty, theta0_produce)
theta0_full = merge_thetas(theta0s, dataset_full)

if length(THETA0_FULL_OVERRIDE) == length(theta0_full)
    theta0_full .= THETA0_FULL_OVERRIDE
elseif length(THETA0_FULL_OVERRIDE) != 0
    @warn "Not using theta supplied - incorrect length"
end

if false # CONVERT_KAPPA
    let r = theta_royalty(dataset_full, theta0_full),
        kap = theta_royalty_κ(data_royalty, r)
        update_kappa_level_to_cumsum!(kap)
    end
end

let thts = split_thetas(dataset_full, theta0_full)
    println_time_flush("Parameters are")
    println("theta0_drill = $(round.(thts[1]; digits=3))")
    println("theta0_roy   = $(round.(thts[2]; digits=3))")
    println("theta0_pdxn  = $(round.(thts[3]; digits=3))")
end

if DO_FULL
    println_time_flush("Starting full model solution")
    res_u, ew_u = solve_model(dataset_full, theta0_full, M_full, MAXTIME_FULL)
else
    println_time_flush("Evaluating likelihood only")
    ew_u = evaluate_likelihood(dataset_full, theta0_full, M_full)
end

# ------------------- save stuff -----------------------

datapathout = replace(DATAPATH, r"\.RData" => ".jld2")
fn = "estimation-results-" * SLURM_JOBID * "-" * HASRIGS * "-" * datapathout

leo = LocalEstObj(ew_u)
reo = RemoteEstObj(ew_u)

leodata = ShaleDrillingLikelihood.data(leo)
leodatadrill = ShaleDrillingLikelihood.drill(leodata)
ddmnovf = DDM_NoVF(_model(leodatadrill))

println_time_flush("saving data to\n\t$fn")
# jldopen(fn, "w"; compress=true) do file
jldopen(fn, true, true, true, IOStream; compress=true) do file
    file["DATAPATH"] = DATAPATH
    file["M"]        = ShaleDrillingLikelihood._num_sim(ShaleDrillingLikelihood.sim(reo))
    file["ddm_novf"] = ddmnovf
    file["LL"]       = ShaleDrillingLikelihood.LL(reo)
    file["grad"]     = grad(leo)
    file["hess"]     = hess(leo)
    file["theta1"]   = theta1(leo)
    file["invhess"]  = invhess(leo)
end
println_time_flush("Done!")


text:
module ReferenceImplementations

export @refimpl

using MacroTools
using MacroTools: postwalk
using Cassette
using Test

struct RefImpl end
struct RefImplAll end  # used to dispatch when invoking all reference implementations
Cassette.@context RefImplCtx

# used as ReferenceImplementations.overdub, etc. in macros to shorten the expressions a bit
const overdub = Cassette.overdub
const recurse = Cassette.recurse

# generate a standard RefImplCtx context (no hooks!)
refimplctx(valT) = Cassette.disablehooks(RefImplCtx(metadata=valT))

# define a reference implementation
function refimpl_def(func)
    funcdef = splitdef(func)
    funcargs = copy(funcdef[:args])
    whereparams = funcdef[:whereparams]
    pushfirst!(funcdef[:args], :(::ReferenceImplementations.RefImpl))

    funcname = funcdef[:name]
    newfuncdef = MacroTools.combinedef(funcdef)

    if length(funcdef[:kwargs]) > 0  # we have kwargs
        overdub_block = quote
            ReferenceImplementations.overdub(ctx::ReferenceImplementations.RefImplCtx{Val{TREFIMPL}},
                kwf::Core.kwftype(typeof($funcname)),
                kwargs::Any, func::typeof($funcname), $(funcargs...)) where {
                        TREFIMPL <: ReferenceImplementations.RefImplAll, $(whereparams...)} =
                    ReferenceImplementations.recurse(ctx, kwf, kwargs, func, ReferenceImplementations.RefImpl(), $(funcargs...))
            ReferenceImplementations.overdub(ctx::ReferenceImplementations.RefImplCtx{Val{TREFIMPL}},
                kwf::Core.kwftype(typeof($funcname)),
                kwargs::Any, func::TREFIMPL, $(funcargs...)) where {
                        TREFIMPL <: typeof($funcname), $(whereparams...)} =
                    ReferenceImplementations.recurse(ctx, kwf, kwargs, func, ReferenceImplementations.RefImpl(), $(funcargs...))
        end
    else  # just args
        overdub_block = quote
            ReferenceImplementations.overdub(ctx::ReferenceImplementations.RefImplCtx{Val{TREFIMPL}}, func::typeof($funcname),
                $(funcargs...)) where {
                    TREFIMPL <: ReferenceImplementations.RefImplAll, $(whereparams...)} =
                ReferenceImplementations.recurse(ctx, func, ReferenceImplementations.RefImpl(), $(funcargs...))
            ReferenceImplementations.overdub(ctx::ReferenceImplementations.RefImplCtx{Val{TREFIMPL}}, func::TREFIMPL,
                $(funcargs...)) where {
                    TREFIMPL <: typeof($funcname), $(whereparams...)} =
                ReferenceImplementations.recurse(ctx, func, ReferenceImplementations.RefImpl(), $(funcargs...))
        end
    end

    return esc(Expr(:block, newfuncdef, overdub_block))
end

# used to convert kwarg iterator to NamedTuple
extractkwargs(; kwargs...) = values(kwargs)

# functions for overdubbing
refimplall(f, args...) = overdub(refimplctx(Val(RefImplAll)), f, args...)
refimplone(f, ref_func, args...) = overdub(refimplctx(Val(typeof(ref_func))), f, args...)

# evaluate expression with all reference implementations inside
function refimpl_call(ex)
    newex = postwalk(ex) do x

        # throw an error if there is a function definition in the expression
        if @capture(MacroTools.longdef1(x), function (fcall_ | fcall_) body_ end)
            throw(ArgumentError("To define a new reference implementation, the argument of @refimpl must"*
                " contain only a function definition."))
        elseif @capture(x, f_(args__; kwargs__))
            return quote
                ReferenceImplementations.refimplall(Core.kwfunc($(f)),
                    ReferenceImplementations.extractkwargs(;$(kwargs...)), $((f)), $(args...))
            end
        elseif @capture(x, f_(args__))
            return quote
                ReferenceImplementations.refimplall($(f), $(args...))
            end
        else
            return x
        end
    end
    return esc(newex)
end

# evaluate expression with where only one function's reference implementations are invoked
function refimpl_call(ref_func, ex)
    newex = postwalk(ex) do x
        if @capture(MacroTools.longdef1(x), function (fcall_ | fcall_) body_ end)
            throw(ArgumentError("To define a new reference implementation, the argument of @refimpl must"*
                " contain only a function definition."))
        elseif @capture(x, f_(args__; kwargs__))
            return quote
                ReferenceImplementations.refimplone(Core.kwfunc($(f)), $(ref_func),
                    ReferenceImplementations.extractkwargs(;$(kwargs...)), $(f), $(args...))
            end
        elseif @capture(x, f_(args__))
            return quote
                ReferenceImplementations.refimplone($(f), $(ref_func), $(args...))
            end
        else
            return x
        end
    end
    return esc(newex)
end

"""
    @refimpl (method definition)
    @refimpl (expression)
    @refimpl (function name) (expression)

* If prefacing a function definition, defines a reference implementation for that function.
* If prefacing an expression that isn't a function definition, invokes reference implementations of methods defined with [`@refimpl`](@ref), even if nested.

Preface an expression performs a Cassette pass on every top-level function
in the expression, recursively looking for methods with @refimpl implementations.
If a function is passed before the expression (separated by a space), only that
method is switched with its reference implementation.

# Examples

Calling `@refimpl` on an expression calls every method with a reference implementation
in the nested sequence of calls for that expression.

```julia
using ReferenceImplementations
@refimpl mysin(x) = begin println("ref mysin"); return sin(x) end
mysin(x) = begin println("mysin"); return sin(x) end

# call the reference implementation
@refimpl mysin(0.)  # prints "ref mysin"
mysin(0.)           # prints "mysin"
```

This works for `@refimpl` functions that are nested inside other functions in the expression.

```julia
@refimpl f(x) = begin println("ref f"); return mysin(x)^2 end
f(x) = begin println("f"); return mysin(x)^2 end

# call the reference implementation
@refimpl f(0.)  # prints "ref f", "ref mysin"
f(0.)           # prints "f", "mysin"
```

You can target individual functions to be replaced with their reference implementation by passing that function after `@refimpl`.

```julia
@refimpl mysin f(0.)  # prints "f", "ref mysin"
@refimpl f f(0.)      # prints "ref f", "mysin"
```
"""
macro refimpl(ex)
    if @capture(MacroTools.longdef1(ex), function (fcall_ | fcall_) body_ end)
        return refimpl_def(ex)
    else
        return refimpl_call(ex)
    end
end

# invoke a a specific function's reference implementation
macro refimpl(ref_func, ex)
    return refimpl_call(ref_func, ex)
end

end  # module


text:
function buildindex{S, R <: Integer}(invindex::Dict{S, R})
    index = Array(S, length(invindex))
    for (v, i) in invindex
        index[i] = v
    end
    return index
end

function buildinvindex{T}(index::Vector{T}, R=DefaultRefType)
    invindex = Dict{T, R}()
    for (i, v) in enumerate(index)
        invindex[v] = i
    end
    return invindex
end

function buildvalues!{T, R, V}(pool::CategoricalPool{T, R, V})
    n = length(levels(pool))
    resize!(pool.valindex, n)
    for i in 1:n
        pool.valindex[i] = V(i, pool)
    end
    return pool.valindex
end

function buildorder!{S, R <: Integer}(order::Array{R},
                                      invindex::Dict{S, R},
                                      levels::Vector{S})
    for (i, v) in enumerate(levels)
        order[invindex[convert(S, v)]] = i
    end
    return order
end

function buildorder{S, R <: Integer}(invindex::Dict{S, R}, levels::Vector)
    order = Array(R, length(invindex))
    return buildorder!(order, invindex, levels)
end


text:
# The following variables need to be defined before `include`-ing this file:
# `testset_name_suffix`
# `n_values`
# `k_values`
# `m_values`
for T ∈ (ComplexF32, ComplexF64, Complex{Int}, Complex{Int32})
    @time @testset "Matrix Multiply $T $(testset_name_suffix)" begin
        for n ∈ n_values
            for k ∈ k_values
                for m ∈ m_values
                    A = rand(T, m, k)
                    B = rand(T, k, n)

                    Are = real.(A)
                    Bre = real.(B)

                    A′ = permutedims(A)'
                    B′ = permutedims(B)'
                    AB = A * B;
                    A′B = A′*B
                    AB′ = A*B′
                    A′B′= A′*B′

                    AreB = Are*B
                    ABre = A*Bre

                    @info "" T n k m
                    @test @time(Octavian.matmul(A, B)) ≈ AB
                    @test @time(Octavian.matmul(A, Bre)) ≈ ABre
                    @test @time(Octavian.matmul(Are, B)) ≈ AreB
                    @test @time(Octavian.matmul(A′, B)) ≈ A′B
                    @test @time(Octavian.matmul(A, B′)) ≈ AB′
                    @test @time(Octavian.matmul(A′, B′)) ≈ A′B′

                    @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                    @test @time(Octavian.matmul_serial(A, Bre)) ≈ ABre
                    @test @time(Octavian.matmul_serial(Are, B)) ≈ AreB
                    @test @time(Octavian.matmul_serial(A′, B)) ≈ A′B
                    @test @time(Octavian.matmul_serial(A, B′)) ≈ AB′
                    @test @time(Octavian.matmul_serial(A′, B′)) ≈ A′B′

                    C = Matrix{T}(undef, n, m)'
                    @test @time(Octavian.matmul!(C, A, B)) ≈ AB

                    C1 = rand(T, m, n)
                    C2 = copy(C1)
                    α, β = T(1 - 2im), T(3 + 4im)
                    @test @time(Octavian.matmul!(C1, A, B, α, β)) ≈ Octavian.matmul!(C2, A, B, α, β)
                end
            end
        end
    end
end

@time @testset "Matrix Multiply Float64 $(testset_name_suffix)" begin
    T = Float64
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "Matrix Multiply Float32 $(testset_name_suffix)" begin
    T = Float32
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "Matrix Multiply Int32 $(testset_name_suffix)" begin
    T = Int32
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "Matrix Multiply Int64 $(testset_name_suffix)" begin
    T = Int64
    for n ∈ n_values
        for k ∈ k_values
            for m ∈ m_values
                A = rand(T, m, k)
                B = rand(T, k, n)
                A′ = permutedims(A)'
                B′ = permutedims(B)'
                AB = A * B;
                @info "" T n k m
                @test @time(Octavian.matmul(A, B)) ≈ AB
                @test @time(Octavian.matmul(A′, B)) ≈ AB
                @test @time(Octavian.matmul(A, B′)) ≈ AB
                @test @time(Octavian.matmul(A′, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B)) ≈ AB
                @test @time(Octavian.matmul_serial(A, B′)) ≈ AB
                @test @time(Octavian.matmul_serial(A′, B′)) ≈ AB
            end
        end
    end
    m = k = n = max(8Octavian.OCTAVIAN_NUM_TASKS[], 400)
    A = rand(T, m, k);
    B = rand(T, k, n);
    A′ = permutedims(A)';
    B′ = permutedims(B)';
    AB = A * B;
    @test matmul_pack_ab!(similar(AB), A, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A, B′) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B) ≈ AB
    @test matmul_pack_ab!(similar(AB), A′, B′) ≈ AB
end

@time @testset "zero-sized-matrices" begin
    @test Octavian.matmul_serial(randn(0,0), randn(0,0)) == zeros(0, 0)
    @test Octavian.matmul_serial(randn(2,3), randn(3,0)) == zeros(2, 0)
    @test Octavian.matmul_serial(randn(2,0), randn(0,2)) == zeros(2, 2)
    @test Octavian.matmul_serial!(ones(2,2),randn(2,0), randn(0,2), 1.0, 2.0) == ones(2, 2) .* 2
    @test Octavian.matmul(randn(0,0), randn(0,0)) == zeros(0, 0)
    @test Octavian.matmul(randn(2,3), randn(3,0)) == zeros(2, 0)
    @test Octavian.matmul(randn(2,0), randn(0,2)) == zeros(2, 2)
    @test Octavian.matmul!(ones(2,2),randn(2,0), randn(0,2), 1.0, 2.0) == ones(2, 2) .* 2
end


text:
function debug_log_api_call(logger::AbstractLogger, user_model::AVKUserModel, message::String)
    with_logger(logger) do

        # get user email -
        alphavantage_api_email = user_model.alphavantage_api_email

        # current timestamp -
        now_stamp = now()

        # formulate log message -
        log_mesage = "$(now_stamp)::$(alphavantage_api_email)::$(message)"

        # log -
        @debug(log_mesage)
    end
end

max_stars_repo_name: teamalphavantage/AlphaVantageAPI

max_issues_repo_name: teamalphavantage/AlphaVantageAPI

max_forks_repo_name: teamalphavantage/AlphaVantageAPI


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

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

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

     convert_hashes_to_columns!(matrix, symbol_ht)

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

     # exact_sparse_linear_algebra!(matrix, basis)
     linear_algebra!(matrix, basis, Val(linalg), rng)

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

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

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

function standardize_basis!(basis::Basis, ht, ord)
    for i in 1:basis.nlead
        idx = basis.nonred[i]
        # basis.lead[i] = basis.lead[idx]
        basis.nonred[i] = i
        basis.isred[i] = 0
        basis.coeffs[i] = basis.coeffs[idx]
        basis.gens[i] = basis.gens[idx]
    end
    basis.size = basis.ndone = basis.ntotal = basis.nlead
    resize!(basis.coeffs, basis.ndone)
    resize!(basis.gens, basis.ndone)
    resize!(basis.lead, basis.ndone)
    resize!(basis.nonred, basis.ndone)
    resize!(basis.isred, basis.ndone)

    sort_gens_by_lead_increasing_in_standardize!(basis, ht, ord)
    normalize_basis!(basis)
end

function export_basis_data(basis::Basis{C}, ht::MonomialHashtable) where {C<:Coeff}
    exps   = Vector{Vector{ExponentVector}}(undef, basis.nlead)
    coeffs = Vector{Vector{C}}(undef, basis.nlead)

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

    exps, coeffs
end

function hash_to_exponents(basis::Basis, ht::MonomialHashtable)
    exps   = Vector{Vector{ExponentVector}}(undef, basis.nlead)
    for i in 1:basis.nlead
        idx = basis.nonred[i]
        poly = basis.gens[idx]
        exps[i] = Vector{ExponentVector}(undef, length(poly))
        @inbounds for j in 1:length(poly)
            exps[i][j] = ht.exponents[poly[j]]
        end
    end
    exps
end

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

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

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

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

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

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

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

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

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

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

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

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

    basis, basis_ht
end

# Initializes Basis and MonomialHashtable structures,
# fills input data from exponents and coeffs
#
# Hashtable initial size is set to tablesize
function initialize_structures_no_normalize(
            ring::PolyRing,
            exponents::Vector{Vector{ExponentVector}},
            coeffs_qq,
            coeffs_ff::Vector{Vector{C}},
            rng::Random.AbstractRNG,
            tablesize::Int) where {C<:Coeff}

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

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

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

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

    basis, basis_ht
end

# Initializes Basis and MonomialHashtable structures,
# fills input data from exponents and coeffs
#
# Hashtable initial size is set to tablesize
function initialize_structures_ff(
            ring::PolyRing,
            exponents::Vector{Vector{ExponentVector}},
            coeffs::Vector{Vector{C}},
            rng::Random.AbstractRNG,
            tablesize::Int) where {C<:Coeff}

    coeffs_ff = [Vector{CoeffFF}(undef, length(c)) for c in coeffs]
    initialize_structures_no_normalize(ring, exponents, coeffs, coeffs_ff, rng, tablesize)
end

# Initializes Basis and MonomialHashtable structures,
# fills input data from exponents and coeffs
function initialize_structures(
            ring::PolyRing,
            exponents::Vector{Vector{ExponentVector}},
            coeffs_qq::Vector{Vector{CoeffQQ}},
            coeffs_zz::Vector{Vector{CoeffZZ}},
            coeffs::Vector{Vector{C}},
            rng::Random.AbstractRNG,
            tablesize::Int) where {C<:Coeff}

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

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

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

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

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

    basis, basis_ht
end

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

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

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

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

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

    basis, present_ht
end

# Initializes Basis with the given hashed exponents and coefficients
function initialize_structures(
            ring::PolyRing,
            hashedexps::Vector{Vector{Int}},
            coeffs::Vector{Vector{C}},
            present_ht::MonomialHashtable) where {C<:Coeff}

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

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

    basis, present_ht
end

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

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

    for i in 1:basis.nlead
        etmp  = ht.exponents[1]
        etmp  = zero(etmp)
        # etmp is now set to zero, and has a zero hash

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

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

        # add all non redundant elements from basis
        # as matrix upper rows

        matrix.nrows += 1
        uprows[matrix.nrows] = multiplied_poly_to_matrix_row!(
                                    symbol_ht, ht, UInt32(0), etmp,
                                    basis.gens[basis.nonred[i]])

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

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

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

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

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

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

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

        sort_matrix_rows_decreasing!(matrix)

        println(matrix.ncols)

        println(basis.ndone)

        interreduce_matrix_rows!(matrix, basis)

        #=
        @warn "after interreduce"
        dump(matrix, maxdepth=5)
        =#
        # TODO
        # convert_matrix_rows_to_basis_elements_use_symbol!(matrix, basis)
        convert_matrix_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
        # no longer need in two hashtables
        # TODO
        # ht = symbol_ht

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

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

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

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

    end

    # TODO
    # sort_gens_by_lead_increasing_in_reduce!(basis, ht)
end

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

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

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

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

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

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

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

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

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

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

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

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

    sort_matrix_rows_decreasing!(matrix)

    interreduce_matrix_rows!(matrix, basis)

    #=
    @warn "after interreduce"
    dump(matrix, maxdepth=5)
    =#
    # TODO
    # convert_matrix_rows_to_basis_elements_use_symbol!(matrix, basis)
    convert_matrix_rows_to_basis_elements!(matrix, basis, ht, symbol_ht)
    # no longer need in two hashtables
    # TODO
    # ht = symbol_ht

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

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

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

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

    # TODO
    # sort_gens_by_lead_increasing_in_reduce!(basis, ht)
end

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

#=
    Input ivariants:
        - ring is set, and ring.ch == basis.ch, and ring ~ ht
        - divmasks in ht are set
        - basis is filled so that
            basis.ntotal = actual number of elements
            basis.ndone  = 0
            basis.nlead  = 0

    Output invariants:
        - basis.ndone == basis.ntotal == basis.nlead
        - basis.gens and basis.coeffs of size basis.ndone
        - basis elements are sorted increasingly wrt ordering on lead elements
        - divmasks in basis are filled and coincide to divmasks in hashtable

=#
function f4!(ring::PolyRing,
             basis::Basis{Coefftype},
             ht,
             reduced,
             linalg,
             rng) where {Coefftype<:Coeff}

    # print("input: $(basis.ntotal) gens, $(ring.nvars) vars. ")

    @assert ring.ch == basis.ch
    @assert ring.ord == ht.ord && ring.nvars == ht.nvars && ring.explen == ht.explen
    # @error "hashtable divmasks"
    # println(ht.exponents[2:ht.load])
    # println("###########")
    # println(map(x->bitstring(x.divmask), ht.hashdata[2:ht.load]))
    # println("###########")
    @assert basis.ndone == 0

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

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

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

    # @warn "ht initialily" ht.load ht.size

    # makes basis fields valid,
    # does not copy,
    # checks for redundancy of new elems
    plcm = Vector{Int}(undef, 0)
    update!(pairset, basis, ht, update_ht, plcm)

    # @warn "ht update" ht.load ht.size

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

        # selects pairs for reduction from pairset following normal strategy
        # (minimal lcm degrees are selected),
        # and puts these into the matrix rows
        select_normal!(pairset, basis, matrix, ht, symbol_ht)
        # TODO: learn, and select\discard S-polynomials

        # @warn "ht select" ht.load ht.size

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

        # @warn "ht symbolic" ht.load ht.size

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

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

        # @warn "ht update" ht.load ht.size

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

        if d > 10000
            @error "Something is probably wrong in f4. Please submit an issue."
            break
        end
    end

    # remove redundant elements
    filter_redundant!(basis)

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

    standardize_basis!(basis, ht, ht.ord)

    # assertion
    #=
    println("#####")
    println(basis.lead)
    println(map(x->bitstring(x.divmask), ht.hashdata[2:ht.load]))
    println("#####")
    =#
    #=
    for i in 1:basis.nlead
        #=
        println(basis.lead[i])
        println(basis.gens[i][1])
        println(ht.exponents[basis.gens[i][1]])
        println(UInt32(ht.hashdata[basis.gens[i][1]].divmask))
        =#
        @assert basis.lead[i] == ht.hashdata[basis.gens[i][1]].divmask
    end
    =#
    # println("HT: $(ht.load)/$(ht.size)")

    nothing
end


text:
################################################################################
#
#  Print functions
#
################################################################################
function _print_to_file(list::Vector{FieldsTower}, f::String)
  c = open(f, "a+")
  _print_to_file(list, c)
  close(c)
end

function _print_to_file(list::Vector{FieldsTower}, f::IOStream)
  for i = 1:length(list)
    K = list[i].field
    auts = list[i].generators_of_automorphisms
    embs = list[i].subfields
    p = K.pol
    p1 = fmpz[numerator(coeff(p, i)) for i = 0:degree(p)]
    Qx = parent(p)
    imgs = fmpq_poly[Qx(x.prim_img) for x in auts]
    imgs_coeff = Vector{fmpq}[fmpq[coeff(x1, i) for i = 0:degree(x1)] for x1 in imgs]
    embs_imgs = Tuple{fmpq_poly, fmpq_poly}[(codomain(x2).pol, parent(codomain(x2).pol)(x2.prim_img)) for x2 in embs]
    embs_imgs_coeff = Vector{Tuple{Vector{fmpq}, Vector{fmpq}}}(undef, length(embs_imgs))
    for i = 1:length(embs_imgs)
      y, z = embs_imgs[i]
      embs_imgs_coeff[i] = ([coeff(y, j) for j = 0:degree(y)], [coeff(z, k) for k = 0:degree(z)])
    end
    write(f, "$(p1),$(imgs_coeff),$(embs_imgs_coeff)\n")
  end
  return nothing
end

Base.convert(::Type{fmpq}, a::Rational{<: Integer}) = fmpq(a)

function _read_from_file(f::String)
  c = open(f, "r")
  v = _read_from_file(c)
  close(c)
  return v
end

function _read_from_file(f::IOStream)
  Qx, x = PolynomialRing(FlintQQ, "x")
  list = Vector{FieldsTower}()
  for l in eachline(f)
    a = eval(Meta.parse(l))
    K = number_field(Qx(a[1]), "a", cached = false, check = false)[1]
    auts = Vector{NfToNfMor}(undef, length(a[2]))
    for i = 1:length(auts)
      img = K(Qx(a[2][i]))
      auts[i] = NfToNfMor(K, K, img)
    end
    embs = Vector{NfToNfMor}(undef, length(a[3]))
    Kdomain = number_field(x-1)[1]
    Kcodomain = number_field(Qx(a[3][1][1]), cached = false, check = false)[1]
    embs[1] = NfToNfMor(Kdomain, Kcodomain, Kcodomain(Qx(a[3][1][2])))
    for i = 2:(length(a[3]))
      Kdomain = Kcodomain
      Kcodomain = number_field(Qx(a[3][i][1]), cached = false, check = false)[1]
      embs[i] = NfToNfMor(Kdomain, Kcodomain, Kcodomain(Qx(a[3][i][2])))
    end
    push!(list, FieldsTower(K, auts, embs))
  end
  return list
end


text:
### Case control: correcting new parameters in case they cross any feasible bounds: INPUT p_min[N+1], p_max[N+1], p_min[N], p_max[N] OUTPUT p_min_next, p_max_next

function correctparameters(p_min_next, p_max_next, p_min, p_max)

    if(p_min_next < 0)
        p_min_next = 0;
    else
        if(p_min_next > p_max)
            p_min_next = p_max;
        end
    end

    if(p_max_next > 1)
        p_max_next = 1;
    else
        if(p_max_next < p_min)
            p_max_next = p_min;
        end
    end

    return p_min_next, p_max_next

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
"""
    IArray
"""
struct IArray{T,N,P<:AbstractArray{T,N},I<:NTuple{N,AbstractIndex}} <: IndicesArray{T,N,P,I}
    _parent::P
    _indices::I

    function IArray{T,N,P,I}(parent, indices, uc, lc) where {T,N,P,I}
        check_index_params(parent, indices, uc, lc)
        new{T,N,P,I}(parent, indices)
    end
end

function IArray(
    p::AbstractArray{T,N},
    inds::Tuple,
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked,
   ) where {T,N}
    return IArray{T,N,typeof(p)}(
        p,
        map((a, i) -> _process_index(a, i, uc, lc), axes(p), inds),
        AllUnique,
        LengthChecked
       )
end

function IArray(p::AbstractArray{T,N}, axs::Tuple{Vararg{<:Union{Symbol,Nothing},N}}) where {T,N}
    return IArray{T,N,typeof(p)}(p, indices(p, axs), AllUnique, LengthChecked)
end

function IArray(
    p::AbstractVector,
    inds::AbstractVector,
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked,
   )
    return IArray(p, (inds,), uc, lc)
end

function IArray{T,N,P}(
    p::P,
    inds::I,
    uc::Uniqueness,
    lc::AbstractLengthCheck
   ) where {T,N,P<:AbstractArray{T,N},I}
    return IArray{T,N,P,I}(p, inds, uc, lc)
end

#IArray(p::AbstractArray) = IArray(p, indices(p), AllUnique, LengthChecked)

function IArray(x::AbstractArray{T,N}; kwargs...) where {T,N}
    if isempty(kwargs)
        return IArray(x, indices(x), AllUnique, LengthChecked)
    else
        return IArray(x, Tuple([Index{k}(v) for (k,v) in kwargs]))
    end
end

function _process_index(a, i, uc, lc)
    check_index_length(a, i, lc)
    return Index(i, uc)
end

_process_index(a, ::Nothing, uc, lc) = Index(a, uc)
function _process_index(a, i::AbstractIndex, uc, lc)
    check_index_length(a, i, lc)
    return i
end

"""
    IMatrix{T,P,I1,I2}

Alias for two-dimensional `IArray`  with elements of type `T`, parent of
type `P`, and index of type `I1` and `I2`. Alias for [`IArray{T,1}`](@ref).
"""
const IMatrix{T,P<:AbstractMatrix{T},I1,I2} = IArray{T,2,P,Tuple{I1,I2}}

"""
    IVector{T,P,I}

Alias for one-dimensional `IArray`  with elements of type `T`, parent of
type `P`, and index of type `I`. Alias for [`IArray{T,1}`](@ref).
"""
const IVector{T,P<:AbstractVector{T},I1} = IArray{T,1,P,Tuple{I1}}

"""
    IVecOrMat{T}

Union type of [`IVector{T}`](@ref) and [`IMatrix{T}`](@ref).
"""
const IVecOrMat{T} = Union{IMatrix{T},IVector{T}}

const IAdjoint{T,P,I1,I2} = IMatrix{T,Adjoint{T,P},I1,I2}

const ITranspose{T,P,I1,I2} = IMatrix{T,Transpose{T,P},I1,I2}

const IDiagonal{T,V,I1,I2} = IMatrix{T,Diagonal{T,V},I1,I2}

const IQRCompactWY{T,P,I1,I2} = LinearAlgebra.QRCompactWY{T,IMatrix{T,P,I1,I2}}

const IQRPivoted{T,P,I1,I2} = QRPivoted{T,IMatrix{T,P,I1,I2}}

const IQR{T,P,I1,I2} = QR{T,IMatrix{T,P,I1,I2}}

const IQRUnion{T,P,I1,I2} = Union{IQRCompactWY{T,P,I1,I2},
                                  IQRPivoted{T,P,I1,I2},
                                  IQR{T,P,I1,I2}}

Base.parent(x::IArray) = getfield(x, :_parent)

Base.axes(x::IArray) = getfield(x, :_indices)


text:
function sudoku(n, start, ::Val{:raw})
    N = n^2
    d = domain(1:N)

    m = model(;kind=:sudoku)

    # Add variables
    if isnothing(start)
        foreach(_ -> variable!(m, d), 1:(N^2))
    else
        foreach(((x, v),) -> variable!(m, 1 ≤ v ≤ N ? domain(v) : d), pairs(start))
    end

    e = (x; param=nothing, dom_size=N) -> error_f(
        usual_constraints[:all_different])(x; param=param, dom_size=dom_size
    )

    # Add constraints: line, columns; blocks
    foreach(i -> constraint!(m, e, (i * N + 1):((i + 1) * N)), 0:(N - 1))
    foreach(i -> constraint!(m, e, [j * N + i for j in 0:(N - 1)]), 1:N)

    for i in 0:(n - 1)
        for j in 0:(n - 1)
            vars = Vector{Int}()
            for k in 1:n
                for l in 0:(n - 1)
                    push!(vars, (j * n + l) * N + i * n + k)
                end
            end
            constraint!(m, e, vars)
        end
    end

    return m
end

function sudoku(n, start, ::Val{:MOI})
    N = n^2
    m = Optimizer()
    MOI.add_variables(m, N^2)

    # Add domain to variables
    # TODO: make it accept starting solutions
    foreach(i -> MOI.add_constraint(m, VI(i), DiscreteSet(1:N)), 1:N^2)

    # Add constraints: line, columns; blocks
    foreach(i -> MOI.add_constraint(m, VOV(map(VI, (i * N + 1):((i + 1) * N))),
        MOIAllDifferent(N)), 0:(N - 1))
    foreach(i -> MOI.add_constraint(m, VOV(map(VI, [j * N + i for j in 0:(N - 1)])),
            MOIAllDifferent(N)), 1:N)

    for i in 0:(n - 1)
        for j in 0:(n - 1)
            vars = Vector{Int}()
            for k in 1:n
                for l in 0:(n - 1)
                    push!(vars, (j * n + l) * N + i * n + k)
                end
            end
            MOI.add_constraint(m, VOV(map(VI, vars)), MOIAllDifferent(N))
        end
    end

    return m
end

function sudoku(n, start, ::Val{:JuMP})
    N = n^2
    m = JuMP.Model(CBLS.Optimizer)

    @variable(m, 1 ≤ X[1:N, 1:N] ≤ N, Int)
    if !isnothing(start)
        for i in 1:N, j in 1:N
            v_ij = start[i,j]
            if 1 ≤ v_ij ≤ N
                @constraint(m, X[i,j] == v_ij)
            end
        end
    end

    for i in 1:N
        @constraint(m, X[i,:] in AllDifferent()) # rows
        @constraint(m, X[:,i] in AllDifferent()) # columns
    end
    for i in 0:(n-1), j in 0:(n-1)
        @constraint(m, vec(X[(i*n+1):(n*(i+1)), (j*n+1):(n*(j+1))]) in AllDifferent()) # blocks
    end
    return m, X
end

"""
    sudoku(n; start= Dictionary{Int, Int}(), modeler = :JuMP)
Create a model for the sudoku problem of domain `1:n²` with optional starting values. The `modeler` argument accepts :raw, :MOI, and :JuMP (default), which refer respectively to the solver internal model, the MathOptInterface model, and the JuMP model.

```julia
# Construct a JuMP model `m` and its associated matrix `grid` for sudoku 9×9
m, grid = sudoku(3)

# Same with a starting instance
instance = [
    9  3  0  0  0  0  0  4  0
    0  0  0  0  4  2  0  9  0
    8  0  0  1  9  6  7  0  0
    0  0  0  4  7  0  0  0  0
    0  2  0  0  0  0  0  6  0
    0  0  0  0  2  3  0  0  0
    0  0  8  5  3  1  0  0  2
    0  9  0  2  8  0  0  0  0
    0  7  0  0  0  0  0  5  3
]
m, grid = sudoku(3, start = instance)

# Run the solver
optimize!(m)

# Retrieve and display the values
solution = value.(grid)
display(solution, Val(:sudoku))
```
"""
sudoku(n; start=nothing, modeler = :JuMP) = sudoku(n, start, Val(modeler))

@doc raw"""
```julia
mutable struct SudokuInstance{T <: Integer} <: AbstractMatrix{T}
```
A `struct` for SudokuInstances, which is a subtype of `AbstractMatrix`.
---
```julia
SudokuInstance(A::AbstractMatrix{T})
SudokuInstance(::Type{T}, n::Int) # fill in blank sudoku of type T
SudokuInstance(n::Int) # fill in blank sudoku of type Int
SudokuInstance(::Type{T}) # fill in "standard" 9×9 sudoku of type T
SudokuInstance() # fill in "standard" 9×9 sudoku of type Int
SudokuInstance(n::Int, P::Pair{Tuple{Int, Int}, T}...) where {T <: Integer} # construct a sudoku given pairs of coordinates and values
SudokuInstance(P::Pair{Tuple{Int, Int}, T}...) # again, default to 9×9 sudoku, constructing given pairs
```
Constructor functions for the `SudokuInstance` `struct`.
"""
mutable struct SudokuInstance{T <: Integer} <: AbstractMatrix{T}
    A::AbstractMatrix{T} where {T <: Integer}

    function SudokuInstance(A::AbstractMatrix{T}) where {T <: Integer}
        size(A, 1) == size(A, 2) || throw(error("Sodokus must be square; received matrix of size $(size(A, 1))×$(size(A, 2))."))
        isequal(sqrt(size(A, 1)), isqrt(size(A, 1))) || throw(error("SudokuInstances must be able to split into equal boxes (e.g., a 9×9 SudokuInstance has three 3×3 squares).  Size given is $(size(A, 1))×$(size(A, 2))."))
        new{T}(A)
    end
    # # fill in blank sudoku if needed
    # SudokuInstance(::Type{T}, n::Int) where {T <: Integer} = new{T}(fill(zero(T), n, n))
    # SudokuInstance(n::Int) = new{Int}(SudokuInstance(Int, n))
    # # Use "standard" 9×9 if no size provided
    # SudokuInstance(::Type{T}) where {T <: Integer} = new{T}(SudokuInstance(T, 9))
    # SudokuInstance() = new{Int}(SudokuInstance(9))
    # # Construct a sudoku given coordinates and values
    # function SudokuInstance(n::Int, P::Pair{Tuple{Int,Int},T}...) where {T <: Integer}
    #     A = zeros(T, n, n)
    #     for (i, v) in P
    #         A[i...] = v
    #     end
    #     new{T}(A)
    # end
    # # again, default to 9×9
    # SudokuInstance(P::Pair{Tuple{Int,Int},T}...) where {T <: Integer} = new{T}(SudokuInstance(9, P...))
end

"""
    SudokuInstance(X::Dictionary)

Construct a `SudokuInstance` with the values `X` of a solver as input.
"""
function SudokuInstance(X::Dictionary)
    n = isqrt(length(X))
    A = zeros(Int, n, n)
    for (k,v) in enumerate(Base.Iterators.partition(X, n))
        A[k,:] = v
    end
    return SudokuInstance(A)
end

SudokuInstance(X::Matrix{Float64}) = SudokuInstance(map(Int, X))
# # abstract array interface for SudokuInstance struct
"""
    Base.size(S::SudokuInstance)

Extends `Base.size` for `SudokuInstance`.
"""
Base.size(S::SudokuInstance) = size(S.A)

"""
    Base.getindex(S::SudokuInstance, i::Int)
    Base.getindex(S::SudokuInstance, I::Vararg{Int,N}) where {N}

Extends `Base.getindex` for `SudokuInstance`.
"""
# Base.getindex(S::SudokuInstance, i::Int) = getindex(S.A, i)
Base.getindex(S::SudokuInstance, I::Vararg{Int,N}) where {N} = getindex(S.A, I...)

"""
    Base.setindex!(S::SudokuInstance, v, i::Int)
    Base.setindex!(S::SudokuInstance, v, I::Vararg{Int,N})

Extends `Base.setindex!` for `SudokuInstance`.
"""
# Base.setindex!(S::SudokuInstance, v, i) = setindex!(S.A, v, i)
# Base.setindex!(S::SudokuInstance, v, I::Vararg) = setindex!(S.A, v, I...)

const _rules = Dict(
    :up_right_corner => '┐',
    :up_left_corner => '┌',
    :bottom_left_corner => '└',
    :bottom_right_corner => '┘',
    :up_intersection => '┬',
    :left_intersection => '├',
    :right_intersection => '┤',
    :middle_intersection => '┼',
    :bottom_intersection => '┴',
    :column => '│',
    :row => '─',
    :blank => '⋅',  # this is the character used for 0s in a SudokuInstance puzzle
)

"""
    _format_val(a)

Format an integer `a` into a string for SudokuInstance.
"""
_format_val(a) = iszero(a) ? _rules[:blank] : string(a)

"""
    _format_line_segment(r, col_pos, M)

Format line segment of a sudoku grid.
"""
function _format_line_segment(r, col_pos, M)
    sep_length = length(r)

    line = string()
    for k in axes(r, 1)
        n_spaces = 1
        Δ = maximum((ndigits(i) for i in M[:, (col_pos * sep_length) + k])) - ndigits(r[k])
        if Δ ≥ 0
            n_spaces = Δ + 1
        end
        line *= repeat(' ', n_spaces) * _format_val(r[k])
    end

    return line * ' ' * _rules[:column]
end

"""
    _format_line(r, M)

Format line of a sudoku grid.
"""
function _format_line(r, M)
    sep_length = isqrt(length(r))

    line = _rules[:column]
    for i in 1:sep_length
        abs_sep_pos = sep_length * i
        line *= _format_line_segment(r[(abs_sep_pos - sep_length + 1):abs_sep_pos], i - 1, M)
    end

    return line
end

"""
    _get_sep_line(s, pos_row, M)

Return a line separator.
"""
function _get_sep_line(s, pos_row, M)
    sep_length = isqrt(s)

    # deal with left-most edges
    sep_line = string()
    if pos_row == 1
        sep_line *= _rules[:up_left_corner]
    elseif mod(pos_row, sep_length) == 0
        if pos_row == s
            sep_line *= _rules[:bottom_left_corner]
        else
            sep_line *= _rules[:left_intersection]
        end
    end

    # rest of row seperator; TODO: make less convoluted.  ATM it works, but I think it can be simplified.
    for pos_col in 1:s
        sep_line *= repeat(_rules[:row], maximum((ndigits(i) for i in M[:, pos_col])) + 1)
        if mod(pos_col, sep_length) == 0
            sep_line *= _rules[:row]
            if pos_col == s
                if pos_row == 1
                    sep_line *= _rules[:up_right_corner]
                elseif pos_row == s
                    sep_line *= _rules[:bottom_right_corner]
                else
                    sep_line *= _rules[:right_intersection]
                end
            elseif pos_row == 1
                sep_line *= _rules[:up_intersection]
            elseif pos_row == s
                sep_line *= _rules[:bottom_intersection]
            else
                sep_line *= _rules[:middle_intersection]
            end
        end
    end

    return sep_line
end

@doc raw"""
```julia
display(io::IO, S::SudokuInstance)
display(S::SudokuInstance) # default to stdout
```
Displays an ``n\times n`` SudokuInstance.
"""
function Base.display(io, S::SudokuInstance)
    sep_length = isqrt(size(S, 1))
    max_n_digits = maximum((ndigits(i) for i in S))

    println(io, _get_sep_line(size(S, 1), 1, S))
    for (i, r) in enumerate(eachrow(S))
        println(io, _format_line(r, S))
        if iszero(mod(i, sep_length))
            println(io, _get_sep_line(size(S, 1), i, S))
        end
    end

    return nothing
end

"""
    Base.display(S::SudokuInstance)

Extends `Base.display` to `SudokuInstance`.
"""
Base.display(S::SudokuInstance) = display(stdout, S)

"""
    Base.display(X::Dictionary)

Extends `Base.display` to a sudoku configuration.
"""
Base.display(X::Dictionary) = display(SudokuInstance(X))

"""
    Base.display(X, Val(:sudoku))

Extends `Base.display` to a sudoku configuration.
"""
Base.display(X, ::Val{:sudoku}) = display(SudokuInstance(X))


text:
using SeisRequests
using Test

@testset "All tests" begin
    include("util.jl")
    include("fdsnws.jl")
    include("irisws.jl")
    include("high_level.jl")
    include("get_stations.jl")
    include("get_events.jl")
    include("get_data.jl")
end


text:
function newBodyFromInertia(nameBody::String, mesh_inertia_info::MeshInertiaInfo)
    com = mesh_inertia_info.com
    m = mesh_inertia_info.mass
    I3 = mesh_inertia_info.tensor_I
    skew_com = Spatial.vector_to_skew_symmetric(com)
    term = SMatrix{3,3,Float64}(I3 + m * skew_com * transpose(skew_com))
    return RigidBody(SpatialInertia(CartesianFrame3D(nameBody), term, com * m, m))
end

function outputJointTransform_ParentChild(body_parent::RigidBody, body_child::RigidBody, evaluated_joint_type_in,
        dh::basic_dh{Float64}=one(basic_dh{Float64}) )

    rot, trans = dh_R_t(dh)
    rot = RotMatrix{3,Float64}(rot)
    j_parent_child = Joint(body_parent.name * "_" * body_child.name, evaluated_joint_type_in)
    x_parent_child = Transform3D(frame_before(j_parent_child), default_frame(body_parent), rot, trans)
    return j_parent_child, x_parent_child
end

function makeInertiaInfo(e_mesh::eMesh{Tri,Nothing}, i_prop::InertiaProperties{Tri})
    return makeInertiaTensor(getTriQuadRule(3), e_mesh.point, e_mesh.tri, i_prop)
end

function makeInertiaInfo(e_mesh::eMesh{Nothing,Tet}, i_prop::InertiaProperties{Tet})
    return makeInertiaTensor(getTetQuadRule(4), e_mesh.point, e_mesh.tet, i_prop)
end

function makeInertiaTensor(quad_rule::TriTetQuadRule{N_ζ,NQ}, point::Vector{SVector{3,Float64}},
        vec_vol_ind::Vector{SVector{N_ζ,Int64}}, i_prop::InertiaProperties) where {N_ζ, NQ}

    d = i_prop.d
    rho = i_prop.rho

    com, mesh_vol = centroidVolumeCombo(point, vec_vol_ind, d)
    tensor_I = zeros(3, 3)
    eye3 = SMatrix{3,3}(1.0I)
    for ind_k = vec_vol_ind
        points_simplex_k = point[ind_k]
        A = asMat(points_simplex_k)
        v_tet = equiv_volume(points_simplex_k, d)
        for k_quad_point = 1:NQ
            ζ_quad = quad_rule.zeta[k_quad_point]
            r = A * ζ_quad - com
            raw_tensor = eye3 * dot(r, r) - (r * transpose(r))
            the_mass = rho * quad_rule.w[k_quad_point] * v_tet
            tensor_I += the_mass * raw_tensor
        end
    end
    return MeshInertiaInfo(tensor_I, com, mesh_vol * rho, mesh_vol)
end

function centroidVolumeCombo(point::Vector{SVector{3,Float64}}, vec_vol_ind::Vector{SVector{N,Int64}},
        d::Union{Nothing,Float64}) where {N}

    v_cum, c_cum = 0.0, 0.0
    for ind_k = vec_vol_ind
        points_simplex_k = point[ind_k]
        v = equiv_volume(points_simplex_k, d)
        v_cum += v
        c_cum += v * centroid(points_simplex_k)
    end
    return c_cum / v_cum, v_cum
end

equiv_volume(v::SVector{4,SVector{3,Float64}}, ::Nothing) = volume(v)
equiv_volume(v::SVector{3,SVector{3,Float64}}, h::Float64) = area(v) * h