text:
# Importing Gurobi in this way avoids a warning with Requires
using .Gurobi: Gurobi

function run_scenario_gurobi(; solver_kwargs::Union{Dict,Nothing}=nothing, kwargs...)
    solver_kwargs = something(solver_kwargs, Dict("Method" => 2, "Crossover" => 0))
    try
        global env = Gurobi.Env()
    catch e
        println("Error encountered starting Gurobi.")
        throw(e)
    end
    try
        global m = run_scenario(;
            optimizer_factory=env, solver_kwargs=solver_kwargs, kwargs...
        )
        Gurobi.finalize(JuMP.backend(m))
    finally
        Gurobi.finalize(env)
        println("Connection to Gurobi closed successfully!")
    end
    # Return `nothing` to prevent `m` from the `try` block from being returned
    return nothing
end

function new_model(env::Gurobi.Env)
    return JuMP.direct_model(Gurobi.Optimizer(env))
end

ext: jl


text:
# FIXME: upstream https://github.com/JuliaGraphics/ColorVectorSpace.jl/issues/75
@inline _nan(::Type{HSV{Float16}}) = HSV{Float16}(NaN16,NaN16,NaN16)
@inline _nan(::Type{HSV{Float32}}) = HSV{Float32}(NaN32,NaN32,NaN32)
@inline _nan(::Type{HSV{Float64}}) = HSV{Float64}(NaN,NaN,NaN)
@inline _nan(::Type{T}) where {T} = nan(T)

# The default values used by extrapolation for off-domain points
const FillType = Union{Number,Colorant,Flat,Periodic,Reflect}
const FloatLike{T<:AbstractFloat} = Union{T,AbstractGray{T}}
const FloatColorant{T<:AbstractFloat} = Colorant{T}
@inline _default_fill(::Type{T}) where {T<:FloatLike} = convert(T, NaN)
@inline _default_fill(::Type{T}) where {T<:FloatColorant} = _nan(T)
@inline _default_fill(::Type{T}) where {T} = zero(T)

@inline _make_compatible(T, fill) = fill
@inline _make_compatible(::Type{T}, fill::Number) where {T} = T(fill)

box_extrapolation(etp::AbstractExtrapolation) = etp

function box_extrapolation(itp::AbstractInterpolation{T}, fill::FillType = _default_fill(T)) where T
    etp = extrapolate(itp, _make_compatible(T, fill))
    box_extrapolation(etp)
end

function box_extrapolation(parent::AbstractArray{T,N}, degree::D = Linear(), args...) where {T,N,D<:Union{Linear,Constant}}
    axs = axes(parent)
    itp = Interpolations.BSplineInterpolation{T,N,typeof(parent),BSpline{D},typeof(axs)}(parent, axs, BSpline(degree))
    box_extrapolation(itp, args...)
end

function box_extrapolation(parent::AbstractArray{T,N}, degree::Interpolations.Degree, args...) where {T,N}
    itp = interpolate(parent, BSpline(degree))
    box_extrapolation(itp, args...)
end

function box_extrapolation(parent::AbstractArray, fill::FillType)
    box_extrapolation(parent, Linear(), fill)
end

function box_extrapolation(itp::AbstractInterpolation, degree::Union{Linear,Constant}, args...)
    throw(ArgumentError("Boxing an interpolation in another interpolation is discouraged. Did you specify the parameter \"$degree\" on purpose?"))
end

function box_extrapolation(itp::AbstractInterpolation, degree::Interpolations.Degree, args...)
    throw(ArgumentError("Boxing an interpolation in another interpolation is discouraged. Did you specify the parameter \"$degree\" on purpose?"))
end

function box_extrapolation(itp::AbstractExtrapolation, fill::FillType)
    throw(ArgumentError("Boxing an extrapolation in another extrapolation is discouraged. Did you specify the parameter \"$fill\" on purpose?"))
end

# This is type-piracy, but necessary if we want Interpolations to be
# independent of OffsetArrays.
function AxisAlgorithms.A_ldiv_B_md!(dest::OffsetArray, F, src::OffsetArray, dim::Integer, b::AbstractVector)
    indsdim = axes(parent(src), dim)
    indsF = axes(F)[2]
    if indsF == indsdim
        AxisAlgorithms.A_ldiv_B_md!(parent(dest), F, parent(src), dim, b)
        return dest
    end
    throw(DimensionMismatch("indices $(axes(parent(src))) do not match $(axes(F))"))
end


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using LinearAlgebra: cross, norm, dot
using StaticArrays

struct CylindricalCoordinateSystem{T}
   ez::SVector{3, T}
   ex::SVector{3, T}
   ey::SVector{3, T}
   o::SVector{3, T}
end

struct CylindricalCoordinates{T}
   cosθ::T
   sinθ::T
   r::T
   z::T
end

function CylindricalCoordinateSystem(R::AbstractVector{T},
                                     o=zero(SVector{3,T})) where {T}
   @assert length(R) == 3
   # first coordinate
   ez = R / norm(R)
   # second coordinate
   if abs(ez[1] - 1) > 0.33
      ex = @SVector T[1, 0, 0]
   else
      ex = @SVector T[0, 1, 0]
   end
   ex = ex - dot(ez, ex) * ez
   ex /= norm(ex)
   # third coordinate
   ey = cross(ex, ez)
   return CylindricalCoordinateSystem(ez, ex, ey, o)
end

(C::CylindricalCoordinateSystem)(R) = cylindrical(C, R)

function cylindrical(C::CylindricalCoordinateSystem{T},
                     R::AbstractVector{T}) where {T}
   @assert length(R) == 3
   Rc = R - C.o
   x = dot(Rc, C.ex)
   y = dot(Rc, C.ey)
   z = dot(Rc, C.ez)
   r = sqrt(x^2 + y^2)
   cosθ = x / r
   sinθ = y / r
   return CylindricalCoordinates(cosθ, sinθ, r, z)
end

cartesian(C::CylindricalCoordinateSystem{T},
          c::CylindricalCoordinates{T}) where {T} =
      C.o + (c.r * c.cosθ * C.ex + c.r * c.sinθ * C.ey + c.z * C.ez)


text:
using StatisticalNonlocality, Distributions, Random, LinearAlgebra
import StatisticalNonlocality: ou_transition_matrix
import StatisticalNonlocality: uniform_phase
import Distributions: Uniform

# Set random seed for reproducibility
Random.seed!(10001)

# number_of_states
number_of_states = 15
# simulate a continuous time markov process
T = ou_transition_matrix(number_of_states - 1)
# T = uniform_phase(4)

γ = 0.01 # γ is basically the dt of the system let's say
eT = exp(γ * T) # get the transition probabilities
# each column gives the transition probability
# column i, means, given that I am in state i, each row j gives the probability to transition to state j

ceT = cumsum(eT, dims=1)

# Define the jump map
function next_state(current_state_index::Int, cT)
    vcT = view(cT, :, current_state_index)
    u = rand(Uniform(0, 1))
    # choose a random uniform variable and decide next state
    # depending on where one lies on the line with respect to the probability
    # of being found between given probabilities
    for i in eachindex(vcT)
        if u < vcT[i]
            return i
        end
    end

    return i
end

M = 1000000
markov_chain = zeros(Int, M)
markov_chain[1] = 4
for i = 2:M
    markov_chain[i] = next_state(markov_chain[i-1], ceT)
end

# construct transition probability from markov chain
constructed_transition_probability = zeros(number_of_states, number_of_states)
for i in 1:M-1
    local current_state = markov_chain[i]
    local next_state = markov_chain[i+1]
    constructed_transition_probability[next_state, current_state] += 1
end

normalization = sum(constructed_transition_probability, dims=1)
constructed_transition_probability = constructed_transition_probability ./ normalization

#=
using GLMakie
fig2 = scatter(markov_chain)
fig = hist(markov_chain)
=#

# dynamic container for holding times
holding_times = [[] for n in 1:number_of_states]
# dynamic container for states
state_list = []

# initialize
push!(holding_times[markov_chain[1]], γ)
push!(state_list, markov_chain[1])

# loop through steps
for i in 2:M
    current_state = markov_chain[i]
    previous_state = markov_chain[i-1]
    if current_state == previous_state
        holding_times[current_state][end] += γ
    else
        push!(state_list, current_state)
        push!(holding_times[current_state], γ)
    end
end

# construct transition matrix from state list
constructed_T = zeros(number_of_states, number_of_states)
normalization = sum(constructed_T, dims=1)
normalized_T = constructed_T ./ normalization

# first count how often state=i transitions to state=j
# the current state corresponds to a column of the matrix
# the next state corresponds to the row of the matrix
number_of_unique_states = length(state_list)
for i in 1:number_of_unique_states-1
    local current_state = state_list[i]
    local next_state = state_list[i+1]
    constructed_T[next_state, current_state] += 1
end

# now we need to normalize
normalization = sum(constructed_T, dims=1)
normalized_T = constructed_T ./ normalization
# now we account for holding times
holding_scale = 1 ./ mean.(holding_times)
for i in 1:number_of_states
    normalized_T[i, i] = -1.0
    normalized_T[:, i] *= holding_scale[i]
end
empirical_T = normalized_T

empirical_transition_probability = exp(empirical_T * γ)
empirical_T_2 = log(constructed_transition_probability) / γ
exact_transition_probability = exp(T * γ)
empirical_error = norm(empirical_transition_probability - exact_transition_probability) / norm(exact_transition_probability)
empirical_error_2 = norm(constructed_transition_probability - exact_transition_probability) / norm(exact_transition_probability)
empirical_error_T = norm(empirical_T - T) / norm(T)
empirical_error_2_T = norm(empirical_T_2 - T) / norm(T)
empirical_holding_scale_error = [abs(T[i, i] + holding_scale[i]) / abs(T[i, i]) for i in 1:number_of_states]
empirical_holding_scale_error_2 = [abs(T[i, i] - empirical_T_2[i, i]) / abs(T[i, i]) for i in 1:number_of_states]

println("For the default way of constructing the transition probability")
println("The error in constructed the matrix is ", empirical_error)
println("The relative error in holding scale is ", empirical_holding_scale_error)
println("note, cannot estimate holding scale to within an error less than ", γ, "ish", " for the first case")
println("For the second way of constructing the transition probability")
println("The error in constructed the matrix is ", empirical_error_2)
println("The relative error in holding scale is ", empirical_holding_scale_error_2)
println("This is independent of the choice of timestep, but may yield wonky transition matrices after taking the log")

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
using LargeCats
const LC = LargeCats
using Test

@testset "internals" begin
    arr1 = randn(1, 2, 3)
    arr2 = randn(1, 2, 4)

    dims = (3, 2)
    @test (1, 4, 7) === @inferred(LC.largecat_outsize([arr1, arr2], dims = (3, 2)))
    @test Float64 === @inferred(LC.largecat_eltype([arr1, arr2]))

    @test (2,) === @inferred(LC.largecat_outsize([[1, 2]], dims = 1))
    @test Int === @inferred(LC.largecat_eltype([[1, 2], [3, 4]]))
    @test Int === @inferred(LC.largecat_eltype(([1], [2])))
    @test Float64 === @inferred(LC.largecat_eltype(([1], [2.2])))

end

struct MyType
    value
end
@testset "exotic eltype" begin
    xs = [
        [MyType(1), MyType("a")], [MyType(:a)]
    ]
    @inferred largecat(xs, dims=1)
    @test largecat(xs, dims=1) == MyType[MyType(1), MyType("a"), MyType(:a)]
end

@testset "against cat" begin
    arr1 = randn(1, 2, 3)
    arr2 = randn(1, 2, 4)
    arrs = [arr1, arr2]
    dims = (2, 3)

    arr = @inferred largecat(arrs, dims = dims)
    @test eltype(arr) == Float64
    @test largecat(arrs, dims = dims) == cat(arrs..., dims = dims)

    function rand_setup()
        T = rand([Int8, Int, Float32, Float64])
        ndims = rand(1:5)
        dim = rand(1:ndims)
        base_shape = rand(1:10, ndims)
        nshapes = rand(1:50)
        shapes = NTuple{ndims, Int}[]
        for _ in 1:nshapes
            shape = copy(base_shape)
            shape[dim] = rand(1:10)
            push!(shapes, Tuple(shape))
        end
        return (T=T, shapes=shapes, dims=dim)
    end

    for _ in 1:100
        setup = rand_setup()
        arrs = [rand(setup.T, shape) for shape in setup.shapes]
        @test @inferred(largecat(arrs, dims=setup.dims)) == cat(arrs..., dims=setup.dims)
    end
end

if VERSION >= v"1.5"
    @testset "largecat perf " begin
        function printrow(n, rel_time, tot_time, mem, outsize)
            println(
                "|",
                rpad(string(n), 10),
                "|",
                rpad(string(rel_time), 15),
                "|",
                rpad(string(tot_time), 15),
                "|",
                rpad(string(mem), 15),
                "|",
                rpad(string(outsize), 15),
                "|",
            )
        end
        println("Performance counters")
        printrow(
            "n arrays",
            "time/length/ns",
            "total time/s",
            "allocs/bytes",
            "outsize/bytes",
        )

        for n in Int[1, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6]
            arrs = map(1:Int(n)) do i
                rand(Float32, rand(1:3), 1)
            end
            largecat(arrs[1:1], dims = (1,)) # warmup
            stats = @timed largecat(arrs, dims = (1,))
            outsize = Base.summarysize(stats.value)
            rel_time = round(1e9 * stats.time / length(stats.value), digits = 3)
            printrow(n, rel_time, stats.time, stats.bytes, outsize)
            # @info "$(stats.time)s $(stats.bytes) bytes"
            @test outsize < stats.bytes < outsize + 1000
        end

        println("Special situations")
        printrow("setup", "time/length/ns", "total time/s", "allocs/bytes", "outsize/bytes")
        for setup in [
            (
                T = Float32,
                shapes = [
                    (128, 128, 1, 1),
                    (128, 128, 2, 1),
                    (128, 128, 3, 1),
                    (128, 128, 4, 1),
                ],
                dims = 3,
            ),
            (
                T = Float32,
                shapes = [(128, 128, 128, 1, 1), (128, 128, 128, 1, 1)],
                dims = 4,
            ),
        ]
            arrs = map(shape -> randn(setup.T, shape), setup.shapes)
            out = largecat(arrs, dims = setup.dims)
            stats = @timed largecat(arrs, dims = setup.dims)
            outsize = Base.summarysize(stats.value)
            rel_time = round(1e9 * stats.time / length(stats.value), digits = 3)
            printrow(
                "$(length(arrs))x$(ndims(out))d",
                rel_time,
                stats.time,
                stats.bytes,
                outsize,
            )
            @test outsize < stats.bytes < outsize + 1000
        end
    end
end


text:
using Weber
using Base.Test

rng() = MersenneTwister(1983)

@testset "Oddball Design" begin
  pattern = oddball_paradigm(identity,20,150,rng=rng())
  @test sum(pattern) == 20
  @test length(pattern) == 170
  @test !any(pattern[2:end] .& (pattern[2:end] .== pattern[1:(end-1)]))
  @test any(pattern[3:end] .& (pattern[3:end] .== pattern[1:(end-2)]))
  @test !any(pattern[1:20])
  @test_throws AssertionError oddball_paradigm(identity,61,80)
  pattern2 = oddball_paradigm(identity,60,80,rng=rng())
  @test sum(pattern2) == 60
  @test length(pattern2) == 140
  pattern3 = oddball_paradigm(identity,60,140,oddball_spacing=2,rng=rng())
  @test sum(pattern3) == 60
  @test length(pattern3) == 200
  @test !any(pattern3[2:end] .& (pattern3[2:end] .== pattern3[1:(end-1)]))
  @test !any(pattern3[3:end] .& (pattern3[3:end] .== pattern3[1:(end-2)]))
  @test any(pattern3[4:end] .& (pattern3[4:end] .== pattern3[1:(end-3)]))
end


text:
# Required packages
using DataFrames, StatsBase, Printf, LinearAlgebra

# Reading in source files
include("dat.jl");
include("branchAndBound.jl");
include("utilities.jl");
include("multiComponents.jl");

# Pitprops example
k=5; # desired sparsity
prob = pitprops; #problem selection
obj, xVal, timetoBound, timetoConverge, timeOut,
	explored, toPrint = branchAndBound(prob, k) ;

print(@sprintf("objective value: %.4e\n", obj));
print(@sprintf("time to best solution: %.4e\n", timetoBound));
print(@sprintf("time to converge: %.4e\n", timetoConverge));
print(@sprintf("nodes explored: %d\n", explored));
println("solution vector:");
print(xVal);
println("\n");

# A synthetic example
dataset = randn(20,100);
Sigma = dataset'*dataset/20;
syntheticprob = problem(dataset, Sigma);

k=5; # desired sparsity
prob = syntheticprob; #problem selection
obj, xVal, timetoBound, timetoConverge, timeOut,
	explored, toPrint = branchAndBound(prob, k, outputFlag=0) ;

print(@sprintf("objective value: %.4e\n", obj));
print(@sprintf("time to best solution: %.4e\n", timetoBound));
print(@sprintf("time to converge: %.4e\n", timetoConverge));
print(@sprintf("nodes explored: %d\n", explored));
println("solution vector:");
print(xVal);
println("\n");

# Multiple Components
k=5; # desired sparsity
prob = pitprops; #problem selection
ncomp = 4;
all_x, timeOut = multiOptimalSPCA(prob, k, ncomp);

total_variance = adjVarExplained(prob.Sigma, all_x)

println("adjusted variance explained:");
print(total_variance);
print(@sprintf("total time: %.4e\n", timeOut));
println("solution vectors:");
print(all_x);

max_stars_repo_name: ryancorywright/Optimal-SPCA

max_issues_repo_name: ryancorywright/Optimal-SPCA

max_forks_repo_name: ryancorywright/Optimal-SPCA


text:
#*************************************************************************#
#
# Linear Type
#
#*************************************************************************#
"""
    Linear(ninputs, noutputs)

Initialize `ParametricModel` representing linear operator `z = A*x`.
Where `A` is a dense matrix with dimension `(noutputs, ninputs) `
cointaining the reshaped parameter vector `Θ`.
"""
struct Linear <: ParametricModel
    ninputs::Int
    noutputs::Int
    nparams::Int
end

# Constructor
function Linear(ninputs::Int, noutputs::Int)
    nparams = ninputs*noutputs
    Linear(ninputs, noutputs, nparams)
end

#*************************************************************************#
#
# Generate Buffers
#
#*************************************************************************#
function dΘ_buffer(lin::Linear)
    if lin.noutputs == 1
        return Matrix{Float64}(lin.noutputs, lin.nparams)
    else
        return hcat((speye(lin.noutputs) for i=1:lin.ninputs)...)
    end
end

#*************************************************************************#
#
# Compute Output and its Derivatives
#
#*************************************************************************#
function _evaluate!(lin::Linear,
                    x::AbstractVector{<:Real},
                    Θ::AbstractVector{<:Real},
                    z::AbstractVector{Float64},
                    dx::AbstractMatrix{Float64},
                    dΘ::AbstractMatrix{Float64})
    # Get Matrix
    matrix = reshape(Θ, lin.noutputs, lin.ninputs)

    # Compute z
    z .= matrix*x

    # Compute dz
    if !isempty(dx)
        dx .= matrix
    end

    # Compute dΘ
    if !isempty(dΘ)
        for i = 1:lin.ninputs
            for j = 1:lin.noutputs
                dΘ[j, (i-1)*lin.noutputs+j] = x[i]
            end
        end
    end
end

#*************************************************************************#
#
# Initializer
#
#*************************************************************************#
"""
    ScaledGaussianInitializer()

The `ScaledGaussianInitializer` will initialize the parameter vector
with of random values picked from normal distribution with standard
deviation of `1/sqrt(ninputs)`.

Valid only for Linear static model.
"""
immutable ScaledGaussianInitializer <: Initializer
end

function initial_guess(lin::Linear, init::ScaledGaussianInitializer)
    d = Distributions.Normal(0, 1/sqrt(lin.ninputs))

    return rand(d, lin.nparams)
end

function initial_guess(lin::Linear, init::DefaultInitializer)
    initial_guess(lin, ScaledGaussianInitializer())
end


text:
using JuliaSteps
using Test

@testset "JuliaSteps.jl" begin
    # Write your tests here.
end

max_stars_repo_name: hanjinliu/JuliaSteps

max_issues_repo_name: hanjinliu/JuliaSteps

max_forks_repo_name: hanjinliu/JuliaSteps


text:
module ZVSimulator

using Distributions
using MCMC

export
  ZVSample,
  ZVMean,
  sample_serialmc,
  sample_rand,
  zvmean,
  sample_serialmc_zvmean,
  sample_rand_zvmean,
  psim_vrf,
  psim_serialmc_vrf,
  psim_rand_vrf,
  mean_of_zvmean,
  var_of_zvmean,
  vrf,
  collect_diagnostic,
  writedlm_vrf,
  writedlm_reformat

include("zvsample.jl")
include("samplers.jl")
include("psim.jl")
include("stats.jl")
include("io.jl")

end


text:
# # Geometry

# The `WaveProp.Geometry` module handles everything related to the geometrical
# descriptions of curves/surfaces. Its main interface is exposed through
# `AbstractElement`, which represents an `M`-dimensional geometical object
# embedded in ℝᴺ. In abstract terms, an *element* τ is a subset of ℝᴺ described
# by
# * a reference element ``\tau \subset \mathbb{R}^N``
# * a *push-forward map* χ : τ̂ → τ ⊂ ℝᴺ.

# Concrete implementations of `AbstractElement`s are responsible for
# implementing a few crucial methods used to derive almost all functionality of
# an `el::AbstractElement` used throughout this package. In the [next section](#interface) we
# describe this interface.

# ## Interface

# The interface is composed of the following methods:
# - `domain(el)` : return the reference element `R::AbstractReferenceElement`
# - `el(u)` : evaluate the parametrization of the element at (parametric).
#   coordinate `u ∈ domain(el)`.
# - `jacobian(u)` : evaluate the jacobian of the parametrization at (parametric)
#   coordinate `u ∈ domain(el)`.

# ## Lagrangian elements

# An concrete implementation of `AbstractElement` is the parametric type
# `LagrangeElement{D,Np,N,T}`. The type parameters represent:
# - `D` : the domain of the element. Thus, `domain(el)` return an instance of
#   the singleton type `R`.
# - `Np` : the number of points/DOF which describe the element. For example, for
#   a flat triangle, `Np=3`, while for a flat quadrilateral element `Np=4`.
# - `N` : the dimension of the ambient/embedding space.
# - `T` : the data type used in the presentation of `el`. Defaults to `Float64`.

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
module MultDep

using Hecke
import Base.*

function multiplicative_group_mod_units_fac_elem(A::Vector{nf_elem}; use_max_ord::Bool = false)
  k = parent(A[1])
  @assert all(i->parent(i) === k, A)
  if use_max_ord
    zk = maximal_order(k)
    cp = coprime_base(A, zk)
  else
    cp = coprime_base(A)
  end
  sort!(cp, lt = (a,b) -> norm(a) > norm(b))
  M = sparse_matrix(FlintZZ)
  for a = A
    T = Tuple{Int, fmpz}[]
    for i = 1:length(cp)
      p = cp[i]
      v = valuation(a, p)
      if v != 0
        push!(T, (i, fmpz(v)))
      end
#      isone(I) && break
    end
    push!(M, sparse_row(FlintZZ, T))
  end
  h, t = Hecke.hnf_kannan_bachem(M, Val{true}, truncate = true)
  return h, t, cp
end

function units(h::SMat, t, b::Vector{nf_elem})
  u = FacElem{nf_elem, AnticNumberField}[]
  for i=nrows(h)+1:length(b)
    k = [fmpz(0) for i=b]
    k[i] = 1
    for i in length(t):-1:1
      Hecke.apply_right!(k, t[i])
    end
    push!(u, FacElem(b, k))
  end
  return u
end

function unit_group_mod_torsion_fac_elem(O::NfAbsOrd, u::Vector{FacElem{nf_elem, AnticNumberField}})
  U = Hecke._unit_group_init(O)
  s = signature(O)
  r = s[1] + s[2] - 1
  for y = u
    is_tors, p = is_torsion_unit(y, false, U.tors_prec)
    U.tors_prec = max(p, U.tors_prec)
    if is_tors
      continue
    end
    Hecke._add_unit(U, y)
    if length(U.units) >= r
      break
    end
  end
  if length(U.units) < r
    # maybe use pAdic stuff here...
    error("not complete yet")
  end

  U.full_rank = true

  U.units = Hecke.reduce(U.units, U.tors_prec)

  for y = u
    is_tors, p = is_torsion_unit(y, false, U.tors_prec)
    U.tors_prec = max(p, U.tors_prec)
    if is_tors
      continue
    end
    x = Hecke.reduce_mod_units([y], U)[1]
    is_tors, p = is_torsion_unit(x, false, U.tors_prec)
    U.tors_prec = max(p, U.tors_prec)
    if is_tors
      continue
    end
    Hecke._add_dependent_unit(U, x)
  end
  return U
end

function *(O1::NfAbsOrd, O2::NfAbsOrd)
  k = nf(O1)
  @assert k === nf(O2)
  b1 = basis(O1, k)
  n = length(b1)
  b2 = basis(O2, k)
  p = [x*y for (x,y) in Base.Iterators.ProductIterator((b1, b2))]
  d = reduce(lcm, [denominator(x) for x = p])
  M = zero_matrix(FlintZZ, n*n, n)
  z = fmpz()
  for i = 1:n*n
    a = p[i]*d
    Hecke.elem_to_mat_row!(M, i, z, a)
  end
  h = hnf(M)
  b = nf_elem[]
  for i=1:n
    push!(b, Hecke.elem_from_mat_row(k, h, i, d))
  end
  return Hecke.Order(k, b)
end

mutable struct GeIdeal
  a::NfAbsOrdIdl
  function GeIdeal(a::NfAbsOrdIdl)
    o =order(a)
    if o.is_maximal == 1
      return new(a)
    end
    #keep track of known maximality and use this to speed things up
    o = ring_of_multipliers(a)
    if o === order(a)
      return new(a)
    else
      return new(extend(a, o))
    end
  end
  function GeIdeal(a::nf_elem)
    o = Hecke.any_order(parent(a))
    d = denominator(a, o)
    return new(o(a*d)*o), new(o(d)*o)
  end
end

import Hecke.gcd, Hecke.isone, Hecke.*, Hecke.gcd_into!, Hecke.copy, Hecke.divexact,
       Hecke.is_unit, Hecke.coprime_base, Hecke.valuation

function make_compatible!(a::GeIdeal, b::GeIdeal)
  o1 = order(a.a)
  o2 = order(b.a)
  if o1 === o2
    return
  end
  o = o1*o2
  a.a = extend(a.a, o)
  b.a = extend(b.a, o)
end

#is_unit, divexact, gcd, isone, *, gcd_into!, copy
isone(a::GeIdeal) = isone(a.a)
is_unit(a::GeIdeal) = isone(a)
copy(a::GeIdeal) = GeIdeal(a.a)

function gcd(a::GeIdeal, b::GeIdeal)
  make_compatible!(a, b)
  return GeIdeal(a.a + b.a)
end

function gcd_into!(a::GeIdeal, b::GeIdeal, c::GeIdeal)
  make_compatible!(b, c)
  a.a = b.a + c.a
  return a
end

function *(a::GeIdeal, b::GeIdeal)
  make_compatible!(a, b)
  return GeIdeal(a.a * b.a)
end

function divexact(a::GeIdeal, b::GeIdeal)
  make_compatible!(a, b)
  return GeIdeal(divexact(a.a, b.a))
end

Hecke.norm(a::GeIdeal) = norm(a.a)

function coprime_base(A::Vector{nf_elem})
  c = Vector{GeIdeal}()
  for a = A
    n,d = GeIdeal(a)
    isone(n) || push!(c, n)
    isone(d) || push!(c, d)
  end
  return coprime_base(c)
end

function coprime_base(A::Vector{nf_elem}, O::NfAbsOrd)
  c = Vector{NfAbsOrdIdl{AnticNumberField, nf_elem}}()
  for a = A
    n,d = integral_split(a*O)
    isone(n) || push!(c, n)
    isone(d) || push!(c, d)
  end
  return coprime_base(c)
end

function valuation(a::nf_elem, p::GeIdeal)
  return valuation(a, p.a)
end

function mult_syzygies_units(A::Vector{FacElem{nf_elem, AnticNumberField}})
  p = next_prime(100)
  K = base_ring(parent(A[1]))
  m = maximum(degree, keys(factor(K.pol, GF(p)).fac))
  while m > 4
    p = next_prime(p)
    m = maximum(degree, keys(factor(K.pol, GF(p)).fac))
  end
         #experimentally, the runtime is dominated by log
  u = FacElem{nf_elem, AnticNumberField}[]
  prec = 10

  r1, r2 = signature(K)
  r = r1+r2 -1
  n = degree(K)
  C = qAdicConj(K, p)
  la = matrix(conjugates_log(A[1], C, prec, all = false, flat = true))'
  lu = zero_matrix(base_ring(la), 0, n)
  uu = []
  for a = A
    while true
      @vtime :qAdic 1 la = matrix(conjugates_log(a, C, prec, all = false, flat = true))'
      if iszero(la)
        @vtime :qAdic 1 @hassert :qAdic 1 verify_gamma([a], [fmpz(1)], fmpz(p)^prec)
        @vprint :qAdic 1 println("torsion found")
        break
      end
      lv = vcat(lu, la)
      #check_precision and change
      if false && any(x->precision(x) < prec, lv)
        println("loss of precision - not sure what to do")
        for i=1:rows(lv)
          for j = cols(lv) #seems to not do anything
            lv[i, j] = setprecision(lv[i, j], min_p)
            @assert precision(lv[i,j]) == min_p
          end
        end
      end
      @vtime :qAdic 1 k = Hecke.left_kernel_basis(lv)
      @assert length(k) < 2
      if length(k) == 0
        println("new ")
        push!(u, a)
        lu = vcat(lu, la)
        @assert length(u) <= r
      else # length == 1 extend the module
        s = fmpq[]
        for x = k[1]
          @vtime :qAdic 1 y = lift_reco(FlintQQ, x, reco = true)
          if y == nothing
            prec *= 2
            @vprint :qAdic 1  "increase prec to ", prec
            lu = matrix([conjugates_log(x, C, prec, all = false, flat = true) for x = u])'
            break
          end
          push!(s, y)
        end
        if length(s) < length(k[1])
          continue
        end
        d = reduce(lcm, map(denominator, s))
        gamma = fmpz[FlintZZ(x*d)::fmpz for x = s]
        @assert reduce(gcd, gamma) == 1 # should be a primitive relation
        @time if !verify_gamma(push!(copy(u), a), gamma, fmpz(p)^prec)
          prec *= 2
          @vprint :qAdic 1 "increase prec to ", prec
          lu = matrix([conjugates_log(x, C, prec, all = false, flat = true) for x = u])'
          continue
        end
        @assert length(gamma) == length(u)+1
        gamma = vcat(gamma[1:length(u)], [0 for i=length(u)+1:r+length(uu)], [gamma[end]])
        push!(uu, (a, gamma))
      end
      break
    end
  end
  #=
    let u_1, .., u_n be units and
       <u_i | i> has rank s and
        r_i in Z^n be such that
          prod u_i^r_i = 1  (OK, sum of the logs is zero)
          rank <r_i | i> = s as well
    so the r_i form a Q-basis for the relations.
    Essentially, the gamma of above are the r_i
    Let [H|0] = [r_i|i]*T be the hnf with trafo, so T in Gl(n, Z)
    Then
      <u_i|i> = <[u_i|i] T>
      [r_i|i] * [u_i|i]^t = 0 (by construction)
      [r_i|i] T inv(T) [u[i] | i] = 0
      [H | 0]   [v_i | i] = 0
      so, since H is triangular(!!) v_1, ... v_n-s = 0
      and <u_i |i> = <v_n-s+1, ..., v_n>

    for the case of n=s+1 this is mostly the "normal" construction.
    Note: as a side, the relations do not have to be primitive.
      If they are, (and n=s+1), then H = 1
  =#

  for i=1:length(uu)-1
    append!(uu[i][2], zeros(FlintZZ, length(uu[end][2])-length(uu[i][2])))
  end
  if length(uu) == 0
    U = matrix(FlintZZ, length(uu), length(uu[end][2]), reduce(vcat, [x[2] for x = uu]))
  else
    U = matrix(FlintZZ, length(uu), length(uu[end][2]), reduce(vcat, [x[2] for x = uu]))
  end
  _, U = hnf_with_transform(U')
  if false
    U = inv(U)
    V = sub(U, 1:rows(U), 1:cols(U)-length(u))
    U = sub(U, 1:rows(U), cols(U)-length(u)+1:cols(U))
    #U can be reduced modulo V...
    Z = zero_matrix(FlintZZ, cols(V), cols(U))
    I = identity_matrix(FlintZZ, cols(U)) * p^(2*prec)
    k = base_ring(A[1])
    A = [ Z V'; I U']
    l = lll(A)
    U = sub(l, cols(V)+1:rows(l), cols(U)+1:cols(l))
    U = lll(U)
  else
    U = lll(U')
  end
  return Hecke._transform(vcat(u, FacElem{nf_elem,AnticNumberField}[FacElem(k(1)) for i=length(u)+1:r], [x[1] for x = uu]), U')
end

function verify_gamma(a::Vector{FacElem{nf_elem, AnticNumberField}}, g::Vector{fmpz}, v::fmpz)
  #knowing that sum g[i] log(a[i]) == 0 mod v, prove that prod a[i]^g[i] is
  #torsion
  #= I claim N(1-a) > v^n for n the field degree:
   Let K be one of the p-adic fields involved, set b = a^g
   then log(K(b)) = 0 (v = p^l) by assumption
   so val(log(K(b))) >= l, but
   val(X) = 1/deg(K) val(norm(X)) for p-adics
   This is true for all completions involved, and sum degrees is n
 =#

  t = prod([a[i]^g[i] for i=1:length(a)])
  # t is either 1 or 1-t is large, norm(1-t) is div. by p^ln
  #in this case T2(1-t) is large, by the arguments above: T2>= (np^l)^2=:B
  # and, see the bottom, \|Log()\|_2^2 >= 1/4 arcosh((B-2)/2)^2
  B = ArbField(nbits(v)*2)(v)^2
  B = 1/2 *acosh((B-2)/2)^2
  p = Hecke.upper_bound(log(B)/log(parent(B)(2)), fmpz)
  @vprint :qAdic 1  "using", p, nbits(v)*2
  b = conjugates_arb_log(t, max(-Int(div(p, 2)), 2))
#  @show B , sum(x*x for x = b), is_torsion_unit(t)[1]
  @hassert :qAdic 1 (B > sum(x*x for x = b)) == is_torsion_unit(t)[1]
  return B > sum(x*x for x = b)
end

function lift_reco(::FlintRationalField, a::padic; reco::Bool = false)
  if reco
    u, v, N = getUnit(a)
    R = parent(a)
    fl, c, d = rational_reconstruction(u, prime(R, N-v))
    !fl && return nothing

    x = FlintQQ(c, d)
    if v < 0
      return x//prime(R, -v)
    else
      return x*prime(R, v)
    end
  else
    return lift(FlintQQ, a)
  end
end

Hecke.nrows(A::Matrix{T}) where {T} = size(A)[1]
Hecke.ncols(A::Matrix{T}) where {T} = size(A)[2]


text:
module CrystallineBravaisVendor

using ..Brillouin: BasisLike
using StaticArrays

# ---------------------------------------------------------------------------------------- #
# this is just a raw output of `bravaistype` from Crystalline, to avoid depending on all
# of Crystalline just for `bravaistype`; basically just hardcoded "memoization"; note that
# we do *not* normalize "oA" to "oC" (i.e., this is `bravaistype(_, _, normalize=false)`)

@inline function boundscheck_sgnum(sgnum::Integer, D::Integer)
    if D == 3
        if sgnum ∉ 1:230
            throw(DomainError(sgnum, "sgnum must be between 1 and 230 in dimension 3"))
        end
    elseif D == 2
        if sgnum ∉ 1:17
            throw(DomainError(sgnum, "sgnum must be between 1 and 17 in dimension 2"))
        end
    elseif D == 1
        if sgnum ∉ 1:2
            throw(DomainError(sgnum, "sgnum must be between 1 and 2 in dimension 1"))
        end
    else
        throw(DomainError(D, "dimensions must be 1, 2, or 3"))
    end
end

const BRAVAISTYPE_2D = (
    "mp", "mp", "op", "op", "oc", "op", "op", "op", "oc", "tp", "tp", "tp", "hp", "hp",
    "hp", "hp", "hp"
)
const BRAVAISTYPE_3D = (
    "aP", "aP", "mP", "mP", "mC", "mP", "mP", "mC", "mC", "mP", "mP", "mC", "mP", "mP",
    "mC", "oP", "oP", "oP", "oP", "oC", "oC", "oF", "oI", "oI", "oP", "oP", "oP", "oP",
    "oP", "oP", "oP", "oP", "oP", "oP", "oC", "oC", "oC", "oA", "oA", "oA", "oA", "oF",
    "oF", "oI", "oI", "oI", "oP", "oP", "oP", "oP", "oP", "oP", "oP", "oP", "oP", "oP",
    "oP", "oP", "oP", "oP", "oP", "oP", "oC", "oC", "oC", "oC", "oC", "oC", "oF", "oF",
    "oI", "oI", "oI", "oI", "tP", "tP", "tP", "tP", "tI", "tI", "tP", "tI", "tP", "tP",
    "tP", "tP", "tI", "tI", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tI", "tI",
    "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tI", "tI", "tI", "tI", "tP", "tP",
    "tP", "tP", "tP", "tP", "tP", "tP", "tI", "tI", "tI", "tI", "tP", "tP", "tP", "tP",
    "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tP", "tI", "tI",
    "tI", "tI", "hP", "hP", "hP", "hR", "hP", "hR", "hP", "hP", "hP", "hP", "hP", "hP",
    "hR", "hP", "hP", "hP", "hP", "hR", "hR", "hP", "hP", "hP", "hP", "hR", "hR", "hP",
    "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP",
    "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "hP", "cP", "cF",
    "cI", "cP", "cI", "cP", "cP", "cF", "cF", "cI", "cP", "cI", "cP", "cP", "cF", "cF",
    "cI", "cP", "cP", "cI", "cP", "cF", "cI", "cP", "cF", "cI", "cP", "cP", "cP", "cP",
    "cF", "cF", "cF", "cF", "cI", "cI"
)
function bravaistype(sgnum::Integer, Dᵛ::Val{D}=Val(3)) where D
    @boundscheck boundscheck_sgnum(sgnum, D)
    if D == 3
        return BRAVAISTYPE_3D[sgnum]
    elseif D == 2
        return BRAVAISTYPE_2D[sgnum]
    else
        return "lp"
    end
end
bravaistype(sgnum::Integer, D::Integer) = bravaistype(sgnum, Val(D))

# ---------------------------------------------------------------------------------------- #
# vendor the tools from Crystalline needed to convert back and forth between
# conventional/reciprocal lattice bases

# ... relies on `bravaistype` *not* normalizing "oA" to "oC"
centering(sgnum::Integer, Dᵛ::Val{D}) where D = last(bravaistype(sgnum, Dᵛ))

const PRIMITIVE_BASIS_MATRICES = (
    # 1D
    Base.ImmutableDict('p'=>SMatrix{1,1,Float64}(1)),                # primitive
    # 2D
    Base.ImmutableDict('p'=>SMatrix{2,2,Float64}([1 0; 0 1]),        # primitive/simple
                  'c'=>SMatrix{2,2,Float64}([1 1; -1 1]./2)),   # centered
    # 3D
    Base.ImmutableDict(
        'P'=>SMatrix{3,3,Float64}([1 0 0; 0 1 0; 0 0 1]),       # primitive/simple
        'F'=>SMatrix{3,3,Float64}([0 1 1; 1 0 1; 1 1 0]./2),    # face-centered
        'I'=>SMatrix{3,3,Float64}([-1 1 1; 1 -1 1; 1 1 -1]./2), # body-centered
        'R'=>SMatrix{3,3,Float64}([2 -1 -1; 1 1 -2; 1 1 1]./3), # rhombohedrally-centered
        'A'=>SMatrix{3,3,Float64}([2 0 0; 0 1 -1; 0 1 1]./2),   # base-centered (along x)
        'C'=>SMatrix{3,3,Float64}([1 1 0; -1 1 0; 0 0 2]./2))   # base-centered (along z)
    )
@inline function primitivebasismatrix(cntr::Char, ::Val{D}=Val(3)) where D
    D ∉ 1:3 && throw(DomainError(D, "dimension must be 1, 2, or 3"))
    return PRIMITIVE_BASIS_MATRICES[D][cntr]
end

function reciprocalbasis(Rs::BasisLike{D}) where D
    Rm = hcat(Rs...)
    Gm = 2π.*inv(transpose(Rm))
    return SVector{D}(ntuple(i->SVector{D,Float64}(Gm[:,i]), Val(D)))
end

function transform_reciprocal(Gs::BasisLike{D}, P::AbstractMatrix{<:Real}) where D
    Gm′ = hcat(Gs...)/P'
    return SVector{D}(ntuple(i->SVector{D}(Gm′[:,i]), Val(D)))
end

for (f, op_P) in ((:primitivize_reciprocal, :identity), (:conventionalize_reciprocal, :inv))
    @eval function $f(basis::BasisLike{D}, sgnum::Integer) where D
        cntr = centering(sgnum, Val(D))
        return $f(basis, cntr)
    end
    @eval function $f(Gs::BasisLike{D}, cntr::Char) where D
        if cntr == 'P' || cntr == 'p' # the conventional and primitive bases coincide
            return Gs
        else
            P = primitivebasismatrix(cntr, Val(D))
            return transform_reciprocal(Gs, $op_P(P))
        end
    end
end

end # module


text:
module GenerativeAD

using DrWatson
using StatsBase
using Flux
using Statistics
using ImageTransformations
using MLDataPattern

include("datasets/Datasets.jl")
using .Datasets: load_data

include("evaluation/Evaluation.jl")

include("experiments.jl")
include("bayesian_opt.jl")
include("models/Models.jl")

end #module


text:
using LogicCircuits
using Pkg.Artifacts

#####################
# Circuit loaders
#####################

zoo_lc(name, num_classes) =
    load_logistic_circuit(zoo_lc_file(name), num_classes)

zoo_clt_file(name) =
    artifact"circuit_model_zoo" * "/Circuit-Model-Zoo-0.1.2/clts/$name"

zoo_clt(name) =
    parse_clt(zoo_clt_file(name))

zoo_psdd_file(name) =
    artifact"circuit_model_zoo" * "/Circuit-Model-Zoo-0.1.2/psdds/$name"

zoo_psdd(name) =
    load_prob_circuit(zoo_psdd_file(name))


text:
function set_env!(vis, model::AbstractGameModel, i::Int, x::SVx) where {SVx}
    x = standardize(model, x[model.pz[i]])
    set_env!(vis, x)
    return nothing
end

function set_env!(vis, x::VehicleState)
    settransform!(vis["env"], Translation(-x.x, -x.y, 0.))
    return nothing
end

function set_camera!(vis)
    settransform!(vis["/Cameras/default"], Translation(-0.70, 0.0, 0.04))
end

function set_camera_birdseye!(vis; height::T=7.0) where {T}
    t = Translation(0.0, 0.0, height)
    r = compose(
        LinearMap(AngleAxis(pi/2, 0, 0, 1)),
        LinearMap(AngleAxis(1.0*pi, 0, 1, 0)))
    tr = compose(t,r)
    settransform!(vis["/Cameras/default"], tr)
    setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 5.0)
    return nothing
end


text:
using FilePaths
using Test

@testset "FilePaths" begin

include("test_uri.jl")
include("compat.jl")

end


text:
count=0;

function c(x)
    return sin(5000*x)*x;
end

function AdaptiveSimpson(f,a,b; abstol=1.0e-6)
    nf=3;
    ff=map(f,[a,(a+b)/2,b]);
    nf=3;
    I1=(b-a)*dot([1,4,1],ff)/6;
    function adaptrec(f,a,b,ff,I1,tol,nf)
        h=(b-a)/2;
        fm=map(f,[a+h/2,b-h/2]);
        nf=nf+2;

        fR=[ff[2],fm[2],ff[3]];
        fL=[ff[1],fm[1],ff[2]];
        IL=h*dot([1,4,1],fL)/6;
        IR=h*dot([1,4,1],fR)/6;
        I2=IL+IR;
        I=I2+(I2-I1)/15;

        if(abs(I-I2)>tol)
            IL,nf = adaptrec(f,a,a+h,fL,IL,tol/2,nf);
            IR,nf = adaptrec(f,b-h,b,fR,IR,tol/2,nf);
            I=IL+IR;
        end
        return I,nf;
    end
    return adaptrec(f,a,b,ff,I1,abstol,nf);
end

println(AdaptiveSimpson(c,-1,1));
println(count);

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
export FurutaPolicy, action

using ReinforcementLearning
using LinearAlgebra

struct FurutaPolicy{S} <: AbstractPolicy
    action_space::S
end

#Classical controller for the Furuta pendulum.

function FurutaPolicy(;s=nothing)
    FurutaPolicy(s)
end

function (p::FurutaPolicy)(env::AbstractEnv)
    x = state(env)

    sθ, cθ, θ̇ , ϕ, ϕ̇  = x
    θ = angle_norm(env.env.state[1])

    if abs(θ) < 0.5 && abs(θ̇ ) < 5
        L = 1/env.env.params.max_torque .* [-0.8139, -0.0541, -0.0226, -0.0258]

        u = -dot(L, [θ, θ̇ , ϕ, ϕ̇ ])

        if abs(ϕ̇ ) > 0.1
            uf = 0.0076*sign(ϕ̇ )/env.env.params.max_torque
        else
            uf = 0.008*sign(u)/env.env.params.max_torque
        end
        u += uf
    elseif abs(angle_norm(θ-π)) < 0.1 && abs(θ̇ ) < 1
        u = 1.0
    else
        En = (cθ - 1 + θ̇ ^2/(2*15.4520^2))
        u = θ̇ *cθ*En/(10 + abs(θ̇ *cθ*En))
    end

    u = clamp(u, -1.0, 1.0)
    return u
end

function action(p::FurutaPolicy, s; training = false)
    sθ, cθ, θ̇ , ϕ, ϕ̇  = s
    θ = angle_norm(atan(sθ, cθ))
    # θ, θ̇ , ϕ, ϕ̇  = s

    if abs(θ) < 0.5 && abs(θ̇ ) < 5
        L = 1/0.04 .* [-0.8139, -0.0541, -0.0226, -0.0258]

        u = -dot(L, [θ, θ̇ , ϕ, ϕ̇ ])

        if abs(ϕ̇ ) > 0.1
            uf = 0.0076*sign(ϕ̇ )/0.04
        else
            uf = 0.008*sign(u)/0.04
        end
        u += uf
    elseif abs(angle_norm(θ-π)) < 0.1 && abs(θ̇ ) < 1
        u = 1.0
    else
        En = (cθ - 1 + θ̇ ^2/(2*15.4520^2))
        u = θ̇ *cθ*En/(10 + abs(θ̇ *cθ*En))
    end

    u = clamp(u, -1.0, 1.0)
    return u
end

function (p::FurutaPolicy)(s)

    sθ, cθ, θ̇ , ϕ, ϕ̇  = s
    θ = atan(sθ, cθ)

    if abs(θ) < 0.5 && abs(θ̇ ) < 5
        L = 1/p.action_space.right .* [-0.8139, -0.0541, -0.0226, -0.0258]

        u = -dot(L, [θ, θ̇ , ϕ, ϕ̇ ]) + 0.0076*sign(ϕ̇ )/p.action_space.right

        if abs(ϕ̇ ) > 0.1
            uf = 0.0076*sign(ϕ̇ )/p.action_space.right
        else
            uf = 0.008*sign(u)/p.action_space.right
        end
        u += uf
    elseif abs(angle_norm(θ-π)) < 0.1 && abs(θ̇ ) < 1
        u = 1.0
    else
        En = (cθ - 1 + θ̇ ^2/(2*15.4520^2))
        u = θ̇ *cθ*En/(10 + abs(θ̇ *cθ*En))
    end

    u = clamp(u, -1.0, 1.0)
    return u
end

RLBase.update!(p::FurutaPolicy, x) = nothing
angle_norm(x) = Base.mod((x + Base.π), (2 * Base.π)) - Base.π

max_stars_repo_name: polhager/FurutaRL

max_issues_repo_name: polhager/FurutaRL

max_forks_repo_name: polhager/FurutaRL


text:
export PositiveSemidefinite

struct PositiveSemidefinite{T <: Real} <: ConvexSet
    r::T
end
PositiveSemidefinite() = PositiveSemidefinite(0.0)
```
This function projects the symmetric matrix Y onto the closest positive definite matrix.
```

function project!(s::PositiveSemidefinite{<:Real}, v::Matrix{T}, y::Matrix{T}) where {T <: Real}
   (D,V) = LinearAlgebra.eigen(y)
   X = zeros(size(y))
   for j = 1:size(y,1)
      if D[j] >= 0.0
         X = X .+ V[:,j] * D[j] * V[:,j]'
      end
   end
   v = X
end

function project(s::PositiveSemidefinite{<:Real}, y::Matrix{<:Real})
    project!(s, similar(y), y)
end


text:
__precompile__()
module kRPC
include("protobuf/krpc.jl")
using ProtoBuf
using Sockets
using Nullables
using Base64

mutable struct kRPCConnection
    conn::TCPSocket
    stream_conn::TCPSocket
    identifier::Array{UInt8, 1}
    str_listener::Task
    running::Vector{Bool}
end

Base.show(io::IO, z::kRPCConnection) = print(io, "kRPC connection " * base64encode(z.identifier))

struct kStream
    id::UInt64
end

mutable struct kStreamInfo
    listener::Function
    error::Nullable{String}
    typ::Type
    value::Any
end

mutable struct kStreamValue
    value::Any
end

struct kRPCException <: Exception
    message :: AbstractString
end

struct DelayedCall
    invocation::krpc.schema.ProcedureCall
    rtype::Type
end

struct DelayCalls
    calls::Array{DelayedCall, 1}
    DelayCalls()=new(DelayedCall[])
end

struct kEvent end
struct kPC
    pcall::krpc.schema.ProcedureCall
    rty::Union{Type,Nothing}
end

include("types.jl")
include("raw.jl")
include("codegen.jl")
include("streams.jl")
include("connection.jl")

export kRPCConnection, kEvent, kPC, kStream, DelayCalls, kStreamValue
export kRPCConnect, close
export add_stream, add_multiple_streams, clear_streams

macro checkbuilt()
    if isfile(normpath(joinpath(@__DIR__, "../deps/gen.jl")))
        return :(include(normpath(joinpath(@__DIR__, "../deps/gen.jl"))))
    elseif !@isdefined kRPCBuild
        error("kRPC not properly built; run Pkg.build(\"kRPC\") and follow prompts.")
    end
end

end


text:
module Fredholm
    using LinearAlgebra: I, ⋅, norm, Diagonal, mul!

    using ForwardDiff
    using NonNegLeastSquares
    abstract type AutoRegMethod end
    abstract type AbstractRegularization{T} end

    include("regularizations.jl")
    include("autoreg.jl")
    include("solve.jl")

    export invert, XuPei,LCurve,Tikhonov,SecondDerivative
end


text:
using RobotDynamics
using LinearAlgebra
using StaticArrays
using Test
using BenchmarkTools
using Rotations

function LinearQuad(model::Quadrotor, x0, u0, dt::Real)
    F = RobotDynamics.DynamicsJacobian(model)
    z = KnotPoint(x0,u0,dt)
    discrete_jacobian!(RK3, F, model, z)
    G = @MMatrix zeros(n,n-1)
    RobotDynamics.state_diff_jacobian!(G, model, z)
    A = RobotDynamics.get_static_A(F)
    B = RobotDynamics.get_static_B(F)
    A = G'A*G
    B = G'B
    LinearModel(SMatrix(A), SMatrix(B), dt=dt)
end

# Create model
quad = Quadrotor()
n = state_dim(quad)
x0,u0 = zeros(quad)
dt = 0.01
t = 0.0

model = LinearQuad(quad, x0, u0, dt)
A = model.A[1]
B = model.B[1]
x,u = rand(model)
z = KnotPoint(x,u,dt,t)

# Test basic properties
@test RobotDynamics.is_timevarying(model) == false
@test RobotDynamics.is_affine(model) == false
@test size(B) == (state_dim(model), control_dim(model))
@test size(A) == (state_dim(model), state_dim(model))

# Continuous dynamics shouldn't be implemented
@test_throws AssertionError dynamics(model, x, u)
@test_throws AssertionError dynamics(model, x, u, 0.)
@test_throws AssertionError dynamics(model, z)
F = RobotDynamics.DynamicsJacobian(model)
@test_throws MethodError jacobian!(F, model, x, u)
@test_throws MethodError jacobian!(F, model, x, u, 0)
@test_throws AssertionError jacobian!(F, model, z)

# Test discrete dynamics
x2 = discrete_dynamics(PassThrough, model, x, u, t, dt)
@test x2 ≈ A*x + B*u
@test @ballocated(discrete_dynamics(PassThrough, $model, $x, $u, $t, $dt),
    samples=2, evals=2) == 0

discrete_jacobian!(PassThrough, F, model, z)
@test F ≈ [A B]
@test @ballocated(discrete_jacobian!(PassThrough, $F, $model, $z),
    samples=2, evals=2) == 0

@test_throws AssertionError discrete_dynamics(RK3, model, x, u, t, dt)
@test_throws AssertionError discrete_jacobian!(RK3, F, model, z)


text:
using PathWeightSampling
using Test
using DiffEqJump
using StaticArrays
using Catalyst
using ModelingToolkit

system = PathWeightSampling.cooperative_chemotaxis_system(aggregator=Direct(), mmax=9)

joint = PathWeightSampling.reaction_network(system)

@test ModelingToolkit.tosymbol.(Catalyst.species(joint)) ==
      ModelingToolkit.tosymbol.(vcat((PathWeightSampling.independent_species(x) for x in [system.sn, system.rn, system.xn])...))

seed = 1234
conf = PathWeightSampling.generate_configuration(system; seed=seed)

@test conf.s_traj.t[end] == system.dtimes[end]
@test conf.x_traj.t[end] == system.dtimes[end]

full_conf = PathWeightSampling.generate_full_configuration(system; seed=seed)

@test conf.s_traj == full_conf.s_traj
@test conf.x_traj == full_conf.x_traj

cens = PathWeightSampling.ConditionalEnsemble(system)
full_conf2 = PathWeightSampling.sample(conf, cens)
@test full_conf.s_traj == full_conf2.s_traj
@test full_conf.r_traj != full_conf2.r_traj
@test full_conf.x_traj == full_conf2.x_traj


text:
using MultivariateOrthogonalPolynomials, FastTransforms, BlockBandedMatrices, Plots, LinearAlgebra, StaticArrays
pyplot()

####
# Solving
#
#   (Δ + k^2*I) * u = f
#
# in a unit disk |𝐱| ≤ 1,
# using orthogonal polynomials (in x and y) in a disk
# with the weight (1-|𝐱|^2) = (1-x^2-y^2)
####

Z = Zernike(1)
W = Weighted(Z) # w*Z
xy = axes(Z, 1);
x, y = first.(xy), last.(xy);
Δ = Z \ (Laplacian(xy) * W)
S = Z \ W # identity

k = 20
L = Δ + k^2 * S # discretisation of Helmholtz
f = @.(cos(x * exp(y)))

u = W * (L \ (Z \ f))
contourf(u)

# One can also fix the discretisation size

N = 20
Zₙ = Z[:,Block.(1:N)]
Wₙ = W[:,Block.(1:N)]
Lₙ = L[Block.(1:N),Block.(1:N)]

u = Wₙ * (Lₙ \ (Zₙ \ f))
contourf(u)

# We can also do eigenvalues of the Laplacian
N = 20
Δₙ = Δ[Block.(1:N),Block.(1:N)]
Sₙ = S[Block.(1:N),Block.(1:N)]

λ,Q = eigen(Symmetric(Matrix(Δₙ)), Symmetric(Matrix(Sₙ)))

contourf(Wₙ * Q[:,end-10])


text:
@naiveVariationalRule(:node_type     => NAutoregressiveX,
                      :outbound_type => Message{GaussianMeanPrecision},
                      :inbound_types => (Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARXOutNPPPPP)

@naiveVariationalRule(:node_type     => NAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARXIn1PNPPPP)

@naiveVariationalRule(:node_type     => NAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARXIn2PPNPPP)

@naiveVariationalRule(:node_type     => NAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution),
                      :name          => VariationalNARXIn3PPPNPP)

@naiveVariationalRule(:node_type     => NAutoregressiveX,
                      :outbound_type => Message{GaussianWeightedMeanPrecision},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution),
                      :name          => VariationalNARXIn4PPPPNP)

@naiveVariationalRule(:node_type     => NAutoregressiveX,
                      :outbound_type => Message{Gamma},
                      :inbound_types => (ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing),
                      :name          => VariationalNARXIn5PPPPPN)

max_stars_repo_name: wmkouw/NARX

max_issues_repo_name: wmkouw/NARX

max_forks_repo_name: wmkouw/NARX


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SkuInformation <: SwaggerModel
    name::Any # spec type: Union{ Nothing, SkuName } # spec name: name
    tier::Any # spec type: Union{ Nothing, Tier } # spec name: tier
    resourceType::Any # spec type: Union{ Nothing, String } # spec name: resourceType
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    locations::Any # spec type: Union{ Nothing, Vector{String} } # spec name: locations
    capabilities::Any # spec type: Union{ Nothing, Vector{SKUCapability} } # spec name: capabilities
    restrictions::Any # spec type: Union{ Nothing, Vector{Restriction} } # spec name: restrictions

    function SkuInformation(;name=nothing, tier=nothing, resourceType=nothing, kind=nothing, locations=nothing, capabilities=nothing, restrictions=nothing)
        o = new()
        validate_property(SkuInformation, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(SkuInformation, Symbol("tier"), tier)
        setfield!(o, Symbol("tier"), tier)
        validate_property(SkuInformation, Symbol("resourceType"), resourceType)
        setfield!(o, Symbol("resourceType"), resourceType)
        validate_property(SkuInformation, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(SkuInformation, Symbol("locations"), locations)
        setfield!(o, Symbol("locations"), locations)
        validate_property(SkuInformation, Symbol("capabilities"), capabilities)
        setfield!(o, Symbol("capabilities"), capabilities)
        validate_property(SkuInformation, Symbol("restrictions"), restrictions)
        setfield!(o, Symbol("restrictions"), restrictions)
        o
    end
end # type SkuInformation

const _property_map_SkuInformation = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("tier")=>Symbol("tier"), Symbol("resourceType")=>Symbol("resourceType"), Symbol("kind")=>Symbol("kind"), Symbol("locations")=>Symbol("locations"), Symbol("capabilities")=>Symbol("capabilities"), Symbol("restrictions")=>Symbol("restrictions"))
const _property_types_SkuInformation = Dict{Symbol,String}(Symbol("name")=>"SkuName", Symbol("tier")=>"Tier", Symbol("resourceType")=>"String", Symbol("kind")=>"String", Symbol("locations")=>"Vector{String}", Symbol("capabilities")=>"Vector{SKUCapability}", Symbol("restrictions")=>"Vector{Restriction}")
Base.propertynames(::Type{ SkuInformation }) = collect(keys(_property_map_SkuInformation))
Swagger.property_type(::Type{ SkuInformation }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SkuInformation[name]))}
Swagger.field_name(::Type{ SkuInformation }, property_name::Symbol) =  _property_map_SkuInformation[property_name]

const _allowed_SkuInformation_kind = ["Storage", "StorageV2", "BlobStorage", "FileStorage", "BlockBlobStorage"]

function check_required(o::SkuInformation)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    true
end

function validate_property(::Type{ SkuInformation }, name::Symbol, val)
    if name === Symbol("kind")
        Swagger.validate_param(name, "SkuInformation", :enum, val, _allowed_SkuInformation_kind)
    end
end


text:
function bt01()
#
# Brugnano and Trigiante
#
    p = poly([ones(6);-1;-1;im;im;im;-im;-im;-im;2])
    z = [1 6; -1  2; im  3; -im  3; 2  1];
    p, PolyZeros(z)
 end


text:
transD_GP.SkyTEM1DInversion.splitlineconvandlast(soundings, 50, 2; zstart=zstart,
    extendfrac=extendfrac, dz=dz, nlayers=nlayers, figsize=(15,7), prefix=zipsaveprefix)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA

max_forks_repo_name: GeoscienceAustralia/HiQGA


text:
"""
    Action_u1_dm <: Action

Struct containing self energy and vertex components for models with U(1) symmetric Dzyaloshinskii-Moriya interaction.
* `Σ :: Vector{Float64}` : negative imaginary part of the self energy
* `Γ :: Vector{Vertex}`  : Γxx, Γzz, ΓDM (i.e Γxy), Γdd, Γzd, Γdz component of the full vertex
"""
struct Action_u1_dm <: Action
    Σ :: Vector{Float64}
    Γ :: Vector{Vertex}
end

# generate action_u1_dm dummy
function get_action_u1_dm_empty(
    r :: Reduced_lattice,
    m :: Mesh,
    ) :: Action_u1_dm

    # init self energy
    Σ = zeros(Float64, length(m.σ))

    # init vertices
    Γ = Vertex[get_vertex_empty(r, m) for comp in 1 : 6]

    # build action
    a = Action_u1_dm(Σ, Γ)

    return a
end

# init action for symmetric u1 models
function init_action!(
    l :: Lattice,
    r :: Reduced_lattice,
    a :: Action_u1_dm
    ) :: Nothing

    # init bare action for Γxx, Γzz and ΓDM component
    ref_int = SVector{4, Int64}(0, 0, 0, 1)
    ref     = Site(ref_int, get_vec(ref_int, l.uc))

    for i in eachindex(r.sites)
        # get bond from lattice
        b = get_bond(ref, r.sites[i], l)

        # set Γxx bare according to spin exchange
        a.Γ[1].bare[i] = b.exchange[1, 1] / 4.0

        # set Γzz bare according to spin exchange
        a.Γ[2].bare[i] = b.exchange[3, 3] / 4.0

        # set ΓDM bare according to spin exchange
        a.Γ[3].bare[i] = b.exchange[1, 2] / 4.0
    end

    return nothing
end

# set repulsion for u1-dm symmetry
function set_repulsion!(
    A :: Float64,
    a :: Action_u1_dm
    ) :: Nothing

    # init bare action on-site with level repulsion for Γxx
    a.Γ[1].bare[1] = A / 4.0

    # init bare action on-site with level repulsion for Γzz
    a.Γ[2].bare[1] = A / 4.0

    return nothing
end

# helper function to disentangle flags during interpolation for u1 symmetric dm models
function apply_flags_u1_dm(
    b    :: Buffer,
    comp :: Int64
    )    :: Tuple{Float64, Int64}

    # clarify sign of contribution
    sgn = 1.0

    # ξ(μ) * ξ(ν) = -1 for Γzd & Γdz
    if comp in (5, 6) && b.sgn_μν
        sgn *= -1.0
    end

    # -ξ(μ) = -1 for Γdd & Γdz
    if comp in (4, 6) && b.sgn_μ
        sgn *= -1.0
    end

    # -ξ(ν) = -1 for Γdd & Γzd
    if comp in (4, 5) && b.sgn_ν
        sgn *= -1.0
    end

    # clarify which component to interpolate
    if b.exchange_flag
        # ΓDM -> -ΓDM for spin exchange (since ΓDM = Γxy = -Γyx)
        if comp == 3
            sgn *= -1.0
        # Γzd -> Γdz for spin exchange
        elseif comp == 5
            comp = 6
        # Γdz -> Γzd for spin exchange
        elseif comp == 6
            comp = 5
        end
    end

    return sgn, comp
end

# get all interpolated vertex components for u1 symmetric dm models
function get_Γ(
    site :: Int64,
    bs   :: Buffer,
    bt   :: Buffer,
    bu   :: Buffer,
    r    :: Reduced_lattice,
    a    :: Action_u1_dm
    ;
    ch_s :: Bool = true,
    ch_t :: Bool = true,
    ch_u :: Bool = true
    )    :: NTuple{6, Float64}

    Γxx = get_Γ_comp(1, site, bs, bt, bu, r, a, apply_flags_u1_dm, ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)
    Γzz = get_Γ_comp(2, site, bs, bt, bu, r, a, apply_flags_u1_dm, ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)
    ΓDM = get_Γ_comp(3, site, bs, bt, bu, r, a, apply_flags_u1_dm, ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)
    Γdd = get_Γ_comp(4, site, bs, bt, bu, r, a, apply_flags_u1_dm, ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)
    Γzd = get_Γ_comp(5, site, bs, bt, bu, r, a, apply_flags_u1_dm, ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)
    Γdz = get_Γ_comp(6, site, bs, bt, bu, r, a, apply_flags_u1_dm, ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)

    return Γxx, Γzz, ΓDM, Γdd, Γzd, Γdz
end

# get all interpolated vertex components for u1 symmetric dm models on all lattice sites
function get_Γ_avx!(
    r     :: Reduced_lattice,
    bs    :: Buffer,
    bt    :: Buffer,
    bu    :: Buffer,
    a     :: Action_u1_dm,
    temp  :: Array{Float64, 3},
    index :: Int64
    ;
    ch_s  :: Bool = true,
    ch_t  :: Bool = true,
    ch_u  :: Bool = true
    )     :: Nothing

    for comp in 1 : 6
        get_Γ_comp_avx!(comp, r, bs, bt, bu, a, apply_flags_u1_dm, view(temp, :, comp, index), ch_s = ch_s, ch_t = ch_t, ch_u = ch_u)
    end

    return nothing
end

# symmetrize full loop contribution and central part
function symmetrize!(
    r :: Reduced_lattice,
    a :: Action_u1_dm
    ) :: Nothing

    # get dimensions
    num_sites = size(a.Γ[1].ch_s.q2_1, 1)
    num_Ω     = size(a.Γ[1].ch_s.q2_1, 2)
    num_ν     = size(a.Γ[1].ch_s.q2_1, 3)

    # computation for q3
    for v in 1 : num_ν
        for vp in v + 1 : num_ν
            @turbo for w in 1 : num_Ω
                for i in 1 : num_sites
                    # get upper triangular matrix for (v, v') plane for s channel
                    a.Γ[1].ch_s.q3[i, w, v, vp] =  a.Γ[1].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[2].ch_s.q3[i, w, v, vp] =  a.Γ[2].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[3].ch_s.q3[i, w, v, vp] = -a.Γ[3].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[4].ch_s.q3[i, w, v, vp] =  a.Γ[4].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[5].ch_s.q3[i, w, v, vp] = -a.Γ[6].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[6].ch_s.q3[i, w, v, vp] = -a.Γ[5].ch_s.q3[r.exchange[i], w, vp, v]

                    # get upper triangular matrix for (v, v') plane for t channel
                    a.Γ[1].ch_t.q3[i, w, v, vp] =  a.Γ[1].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[2].ch_t.q3[i, w, v, vp] =  a.Γ[2].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[3].ch_t.q3[i, w, v, vp] = -a.Γ[3].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[4].ch_t.q3[i, w, v, vp] =  a.Γ[4].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[5].ch_t.q3[i, w, v, vp] = -a.Γ[6].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[6].ch_t.q3[i, w, v, vp] = -a.Γ[5].ch_t.q3[r.exchange[i], w, vp, v]

                    # get upper triangular matrix for (v, v') plane for u channel
                    a.Γ[1].ch_u.q3[i, w, v, vp] =  a.Γ[1].ch_u.q3[i, w, vp, v]
                    a.Γ[2].ch_u.q3[i, w, v, vp] =  a.Γ[2].ch_u.q3[i, w, vp, v]
                    a.Γ[3].ch_u.q3[i, w, v, vp] =  a.Γ[3].ch_u.q3[i, w, vp, v]
                    a.Γ[4].ch_u.q3[i, w, v, vp] =  a.Γ[4].ch_u.q3[i, w, vp, v]
                    a.Γ[5].ch_u.q3[i, w, v, vp] = -a.Γ[5].ch_u.q3[i, w, vp, v]
                    a.Γ[6].ch_u.q3[i, w, v, vp] = -a.Γ[6].ch_u.q3[i, w, vp, v]
                end
            end
        end
    end

    return nothing
end

# symmetrized addition for left part (right part symmetric to left part)
function symmetrize_add_to!(
    r   :: Reduced_lattice,
    a_l :: Action_u1_dm,
    a   :: Action_u1_dm
    )   :: Nothing

    # get dimensions
    num_sites = size(a_l.Γ[1].ch_s.q2_1, 1)
    num_Ω     = size(a_l.Γ[1].ch_s.q2_1, 2)
    num_ν     = size(a_l.Γ[1].ch_s.q2_1, 3)

    # computation for q1
    @turbo for w in 1 : num_Ω
        for i in 1 : num_sites
            # add q1 to s channel (right part from v <-> v' exchange)
            a.Γ[1].ch_s.q1[i, w] += a_l.Γ[1].ch_s.q1[i, w] + a_l.Γ[1].ch_s.q1[r.exchange[i], w]
            a.Γ[2].ch_s.q1[i, w] += a_l.Γ[2].ch_s.q1[i, w] + a_l.Γ[2].ch_s.q1[r.exchange[i], w]
            a.Γ[3].ch_s.q1[i, w] += a_l.Γ[3].ch_s.q1[i, w] - a_l.Γ[3].ch_s.q1[r.exchange[i], w]
            a.Γ[4].ch_s.q1[i, w] += a_l.Γ[4].ch_s.q1[i, w] + a_l.Γ[4].ch_s.q1[r.exchange[i], w]
            a.Γ[5].ch_s.q1[i, w] += a_l.Γ[5].ch_s.q1[i, w] - a_l.Γ[6].ch_s.q1[r.exchange[i], w]
            a.Γ[6].ch_s.q1[i, w] += a_l.Γ[6].ch_s.q1[i, w] - a_l.Γ[5].ch_s.q1[r.exchange[i], w]

            # add q1 to t channel (right part from v <-> v' exchange)
            a.Γ[1].ch_t.q1[i, w] += a_l.Γ[1].ch_t.q1[i, w] + a_l.Γ[1].ch_t.q1[r.exchange[i], w]
            a.Γ[2].ch_t.q1[i, w] += a_l.Γ[2].ch_t.q1[i, w] + a_l.Γ[2].ch_t.q1[r.exchange[i], w]
            a.Γ[3].ch_t.q1[i, w] += a_l.Γ[3].ch_t.q1[i, w] - a_l.Γ[3].ch_t.q1[r.exchange[i], w]
            a.Γ[4].ch_t.q1[i, w] += a_l.Γ[4].ch_t.q1[i, w] + a_l.Γ[4].ch_t.q1[r.exchange[i], w]
            a.Γ[5].ch_t.q1[i, w] += a_l.Γ[5].ch_t.q1[i, w] - a_l.Γ[6].ch_t.q1[r.exchange[i], w]
            a.Γ[6].ch_t.q1[i, w] += a_l.Γ[6].ch_t.q1[i, w] - a_l.Γ[5].ch_t.q1[r.exchange[i], w]

            # add q1 to u channel (right part from v <-> v' exchange)
            a.Γ[1].ch_u.q1[i, w] += a_l.Γ[1].ch_u.q1[i, w] + a_l.Γ[1].ch_u.q1[i, w]
            a.Γ[2].ch_u.q1[i, w] += a_l.Γ[2].ch_u.q1[i, w] + a_l.Γ[2].ch_u.q1[i, w]
            a.Γ[3].ch_u.q1[i, w] += a_l.Γ[3].ch_u.q1[i, w] + a_l.Γ[3].ch_u.q1[i, w]
            a.Γ[4].ch_u.q1[i, w] += a_l.Γ[4].ch_u.q1[i, w] + a_l.Γ[4].ch_u.q1[i, w]
            a.Γ[5].ch_u.q1[i, w] += a_l.Γ[5].ch_u.q1[i, w] - a_l.Γ[5].ch_u.q1[i, w]
            a.Γ[6].ch_u.q1[i, w] += a_l.Γ[6].ch_u.q1[i, w] - a_l.Γ[6].ch_u.q1[i, w]
        end
    end

    # computation for q2_1 and q2_2
    @turbo for v in 1 : num_ν
        for w in 1 : num_Ω
            for i in 1 : num_sites
                # add q2_1 and q2_2 to s channel (right part from v <-> v' exchange)
                a.Γ[1].ch_s.q2_1[i, w, v] += a_l.Γ[1].ch_s.q2_1[i, w, v] + a_l.Γ[1].ch_s.q2_2[r.exchange[i], w, v]
                a.Γ[2].ch_s.q2_1[i, w, v] += a_l.Γ[2].ch_s.q2_1[i, w, v] + a_l.Γ[2].ch_s.q2_2[r.exchange[i], w, v]
                a.Γ[3].ch_s.q2_1[i, w, v] += a_l.Γ[3].ch_s.q2_1[i, w, v] - a_l.Γ[3].ch_s.q2_2[r.exchange[i], w, v]
                a.Γ[4].ch_s.q2_1[i, w, v] += a_l.Γ[4].ch_s.q2_1[i, w, v] + a_l.Γ[4].ch_s.q2_2[r.exchange[i], w, v]
                a.Γ[5].ch_s.q2_1[i, w, v] += a_l.Γ[5].ch_s.q2_1[i, w, v] - a_l.Γ[6].ch_s.q2_2[r.exchange[i], w, v]
                a.Γ[6].ch_s.q2_1[i, w, v] += a_l.Γ[6].ch_s.q2_1[i, w, v] - a_l.Γ[5].ch_s.q2_2[r.exchange[i], w, v]
                a.Γ[1].ch_s.q2_2[i, w, v] += a_l.Γ[1].ch_s.q2_2[i, w, v] + a_l.Γ[1].ch_s.q2_1[r.exchange[i], w, v]
                a.Γ[2].ch_s.q2_2[i, w, v] += a_l.Γ[2].ch_s.q2_2[i, w, v] + a_l.Γ[2].ch_s.q2_1[r.exchange[i], w, v]
                a.Γ[3].ch_s.q2_2[i, w, v] += a_l.Γ[3].ch_s.q2_2[i, w, v] - a_l.Γ[3].ch_s.q2_1[r.exchange[i], w, v]
                a.Γ[4].ch_s.q2_2[i, w, v] += a_l.Γ[4].ch_s.q2_2[i, w, v] + a_l.Γ[4].ch_s.q2_1[r.exchange[i], w, v]
                a.Γ[5].ch_s.q2_2[i, w, v] += a_l.Γ[5].ch_s.q2_2[i, w, v] - a_l.Γ[6].ch_s.q2_1[r.exchange[i], w, v]
                a.Γ[6].ch_s.q2_2[i, w, v] += a_l.Γ[6].ch_s.q2_2[i, w, v] - a_l.Γ[5].ch_s.q2_1[r.exchange[i], w, v]

                # add q2_1 and q2_2 to t channel (right part from v <-> v' exchange)
                a.Γ[1].ch_t.q2_1[i, w, v] += a_l.Γ[1].ch_t.q2_1[i, w, v] + a_l.Γ[1].ch_t.q2_2[r.exchange[i], w, v]
                a.Γ[2].ch_t.q2_1[i, w, v] += a_l.Γ[2].ch_t.q2_1[i, w, v] + a_l.Γ[2].ch_t.q2_2[r.exchange[i], w, v]
                a.Γ[3].ch_t.q2_1[i, w, v] += a_l.Γ[3].ch_t.q2_1[i, w, v] - a_l.Γ[3].ch_t.q2_2[r.exchange[i], w, v]
                a.Γ[4].ch_t.q2_1[i, w, v] += a_l.Γ[4].ch_t.q2_1[i, w, v] + a_l.Γ[4].ch_t.q2_2[r.exchange[i], w, v]
                a.Γ[5].ch_t.q2_1[i, w, v] += a_l.Γ[5].ch_t.q2_1[i, w, v] - a_l.Γ[6].ch_t.q2_2[r.exchange[i], w, v]
                a.Γ[6].ch_t.q2_1[i, w, v] += a_l.Γ[6].ch_t.q2_1[i, w, v] - a_l.Γ[5].ch_t.q2_2[r.exchange[i], w, v]
                a.Γ[1].ch_t.q2_2[i, w, v] += a_l.Γ[1].ch_t.q2_2[i, w, v] + a_l.Γ[1].ch_t.q2_1[r.exchange[i], w, v]
                a.Γ[2].ch_t.q2_2[i, w, v] += a_l.Γ[2].ch_t.q2_2[i, w, v] + a_l.Γ[2].ch_t.q2_1[r.exchange[i], w, v]
                a.Γ[3].ch_t.q2_2[i, w, v] += a_l.Γ[3].ch_t.q2_2[i, w, v] - a_l.Γ[3].ch_t.q2_1[r.exchange[i], w, v]
                a.Γ[4].ch_t.q2_2[i, w, v] += a_l.Γ[4].ch_t.q2_2[i, w, v] + a_l.Γ[4].ch_t.q2_1[r.exchange[i], w, v]
                a.Γ[5].ch_t.q2_2[i, w, v] += a_l.Γ[5].ch_t.q2_2[i, w, v] - a_l.Γ[6].ch_t.q2_1[r.exchange[i], w, v]
                a.Γ[6].ch_t.q2_2[i, w, v] += a_l.Γ[6].ch_t.q2_2[i, w, v] - a_l.Γ[5].ch_t.q2_1[r.exchange[i], w, v]

                # add q2_1 and q2_2 to u channel (right part from v <-> v' exchange)
                a.Γ[1].ch_u.q2_1[i, w, v] += a_l.Γ[1].ch_u.q2_1[i, w, v] + a_l.Γ[1].ch_u.q2_2[i, w, v]
                a.Γ[2].ch_u.q2_1[i, w, v] += a_l.Γ[2].ch_u.q2_1[i, w, v] + a_l.Γ[2].ch_u.q2_2[i, w, v]
                a.Γ[3].ch_u.q2_1[i, w, v] += a_l.Γ[3].ch_u.q2_1[i, w, v] + a_l.Γ[3].ch_u.q2_2[i, w, v]
                a.Γ[4].ch_u.q2_1[i, w, v] += a_l.Γ[4].ch_u.q2_1[i, w, v] + a_l.Γ[4].ch_u.q2_2[i, w, v]
                a.Γ[5].ch_u.q2_1[i, w, v] += a_l.Γ[5].ch_u.q2_1[i, w, v] - a_l.Γ[5].ch_u.q2_2[i, w, v]
                a.Γ[6].ch_u.q2_1[i, w, v] += a_l.Γ[6].ch_u.q2_1[i, w, v] - a_l.Γ[6].ch_u.q2_2[i, w, v]
                a.Γ[1].ch_u.q2_2[i, w, v] += a_l.Γ[1].ch_u.q2_2[i, w, v] + a_l.Γ[1].ch_u.q2_1[i, w, v]
                a.Γ[2].ch_u.q2_2[i, w, v] += a_l.Γ[2].ch_u.q2_2[i, w, v] + a_l.Γ[2].ch_u.q2_1[i, w, v]
                a.Γ[3].ch_u.q2_2[i, w, v] += a_l.Γ[3].ch_u.q2_2[i, w, v] + a_l.Γ[3].ch_u.q2_1[i, w, v]
                a.Γ[4].ch_u.q2_2[i, w, v] += a_l.Γ[4].ch_u.q2_2[i, w, v] + a_l.Γ[4].ch_u.q2_1[i, w, v]
                a.Γ[5].ch_u.q2_2[i, w, v] += a_l.Γ[5].ch_u.q2_2[i, w, v] - a_l.Γ[5].ch_u.q2_1[i, w, v]
                a.Γ[6].ch_u.q2_2[i, w, v] += a_l.Γ[6].ch_u.q2_2[i, w, v] - a_l.Γ[6].ch_u.q2_1[i, w, v]
            end
        end
    end

    # computation for q3
    @turbo for vp in 1 : num_ν
        for v in 1 : num_ν
            for w in 1 : num_Ω
                for i in 1 : num_sites
                    # add q3 to s channel (right part from v <-> v' exchange)
                    a.Γ[1].ch_s.q3[i, w, v, vp] += a_l.Γ[1].ch_s.q3[i, w, v, vp] + a_l.Γ[1].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[2].ch_s.q3[i, w, v, vp] += a_l.Γ[2].ch_s.q3[i, w, v, vp] + a_l.Γ[2].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[3].ch_s.q3[i, w, v, vp] += a_l.Γ[3].ch_s.q3[i, w, v, vp] - a_l.Γ[3].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[4].ch_s.q3[i, w, v, vp] += a_l.Γ[4].ch_s.q3[i, w, v, vp] + a_l.Γ[4].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[5].ch_s.q3[i, w, v, vp] += a_l.Γ[5].ch_s.q3[i, w, v, vp] - a_l.Γ[6].ch_s.q3[r.exchange[i], w, vp, v]
                    a.Γ[6].ch_s.q3[i, w, v, vp] += a_l.Γ[6].ch_s.q3[i, w, v, vp] - a_l.Γ[5].ch_s.q3[r.exchange[i], w, vp, v]

                    # add q3 to t channel (right part from v <-> v' exchange)
                    a.Γ[1].ch_t.q3[i, w, v, vp] += a_l.Γ[1].ch_t.q3[i, w, v, vp] + a_l.Γ[1].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[2].ch_t.q3[i, w, v, vp] += a_l.Γ[2].ch_t.q3[i, w, v, vp] + a_l.Γ[2].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[3].ch_t.q3[i, w, v, vp] += a_l.Γ[3].ch_t.q3[i, w, v, vp] - a_l.Γ[3].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[4].ch_t.q3[i, w, v, vp] += a_l.Γ[4].ch_t.q3[i, w, v, vp] + a_l.Γ[4].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[5].ch_t.q3[i, w, v, vp] += a_l.Γ[5].ch_t.q3[i, w, v, vp] - a_l.Γ[6].ch_t.q3[r.exchange[i], w, vp, v]
                    a.Γ[6].ch_t.q3[i, w, v, vp] += a_l.Γ[6].ch_t.q3[i, w, v, vp] - a_l.Γ[5].ch_t.q3[r.exchange[i], w, vp, v]

                    # add q3 to u channel (right part from v <-> v' exchange)
                    a.Γ[1].ch_u.q3[i, w, v, vp] += a_l.Γ[1].ch_u.q3[i, w, v, vp] + a_l.Γ[1].ch_u.q3[i, w, vp, v]
                    a.Γ[2].ch_u.q3[i, w, v, vp] += a_l.Γ[2].ch_u.q3[i, w, v, vp] + a_l.Γ[2].ch_u.q3[i, w, vp, v]
                    a.Γ[3].ch_u.q3[i, w, v, vp] += a_l.Γ[3].ch_u.q3[i, w, v, vp] + a_l.Γ[3].ch_u.q3[i, w, vp, v]
                    a.Γ[4].ch_u.q3[i, w, v, vp] += a_l.Γ[4].ch_u.q3[i, w, v, vp] + a_l.Γ[4].ch_u.q3[i, w, vp, v]
                    a.Γ[5].ch_u.q3[i, w, v, vp] += a_l.Γ[5].ch_u.q3[i, w, v, vp] - a_l.Γ[5].ch_u.q3[i, w, vp, v]
                    a.Γ[6].ch_u.q3[i, w, v, vp] += a_l.Γ[6].ch_u.q3[i, w, v, vp] - a_l.Γ[6].ch_u.q3[i, w, vp, v]
                end
            end
        end
    end

    return nothing
end


text:
module Log4jl

export trace, debug, info, warn, error, fatal,
       # @trace, @debug, @info, @warn, @error, @fatal,
       @logger, @Log4jl

import Base: append!, serialize, show, in, delete!, string,
             trace, info, warn, error, filter

# To be set parameters
global LOG4JL_LINE_SEPARATOR
global LOG4JL_DEFAULT_STATUS_LEVEL
global LOG4JL_INTERNAL_STATUS_LEVEL
global LOG4JL_LOG_EVENT
global LOG4JL_CONTEXT_SELECTOR
global LOGGER

const LOG4JL_CONFIG_DEFAULT_PREFIX = "log4jl"
const LOG4JL_CONFIG_EXTS = Dict{Symbol, Vector{AbstractString}}()
const LOG4JL_CONFIG_TYPES = Dict{Symbol, DataType}()

# Imports
include("utils.jl")
include("levels.jl")
include("types.jl")
include("messages.jl")
include("event.jl")
include("filter.jl")
include("layouts.jl")
include("appenders.jl")
include("logconfig.jl")
include("config.jl")
include("logger.jl")
include("context.jl")
include("selector.jl")

# Constants
const LOG4JL_DEFAULT_MESSAGE = Messages.ParameterizedMessage
const ROOT_LOGGER_NAME = ""

function makemethods(fn::Symbol, lvl::Level.EventLevel)
    eval(Log4jl, quote
        $fn(l::AbstractLogger, marker::MARKER, msg, params...) = log(l, $lvl, marker, msg, params...)
        $fn(l::AbstractLogger, marker::Symbol, msg, params...) = log(l, $lvl, MARKER(marker), msg, params...)
        $fn(l::AbstractLogger, msg, params...)                 = log(l, $lvl, msg, params...)
    end)
end

# Logger methods
for (fn,lvl) in ((:trace, Level.TRACE),
                 (:debug, Level.DEBUG),
                 (:info,  Level.INFO),
                 (:warn,  Level.WARN),
                 (:error, Level.ERROR),
                 (:fatal, Level.FATAL))
    makemethods(fn, lvl)
end

"`logger` macro parameters parser"
function parseargs(params, cmname, defmsg=LOG4JL_DEFAULT_MESSAGE)
    config_eval = Nullable{Expr}()
    name = cmname
    msgfactory = :($defmsg)
    config_file = Nullable{AbstractString}()

    length(params) > 0 && for p in params
        trace(LOGGER, "Logger parameter: $p, $(typeof(p))")
        if isa(p, AbstractString)
            name = p
        elseif isa(p, Expr) && p.head == :block
            config_eval = Nullable(p)
        elseif isa(p, Expr) && (p.head == :(=) || p.head ==:kw)
            if p.args[1] == :URI
                config_file = NAME(p.args[2])
            elseif p.args[1] == :MSG
                msgfactory = p.args[2]
            end
        end
    end

    return name, msgfactory, config_file, config_eval
end

"""Create logger

Creates logger instance. It accepts following parameters:

1. `name`: a logger name as string from the configuration
2. `MSG=<message_type>`: a message type used in a configured logger
3. `URI=<config_location>`: a configuration location
4. `begin <configuration> end`: a configuration program (must return `Configuration` object)
"""
macro logger(params...)
    # get current module and its locations
    cm = current_module()
    cmname = string(cm)
    cmdir = moduledir(cm)

    # Parse arguments
    logname, msgfactory, config_file, config_eval = parseargs(params, cmname)

    # find context using a module name
    ctx = getcontext(cmname, cfgloc = config_file, cfgexp = config_eval)

    # return logger
    quote
        logger($ctx, $logname, $msgfactory)
    end
end

macro rootlogger(params...)
    :(@logger $ROOT_LOGGER_NAME $(params...))
end

""" This macro is used to define logger macros for the specified logger: @info, @error, etc.

Macro take a logger assignment expression and creates a variety of macros for the created logger.
```julia
    @Log4jl LOG = @logger # Create logger `LOG` and corresponding macros
    @info "AAA"           # All calls will be made with `LOG` logger
```
"""
macro Log4jl(expr)
    lvar = expr.args[1]
    rval = expr.args[2]
    mod = current_module()
    for fn in [:trace, :debug, :info, :warn, :error, :fatal]
        eval(mod, :( macro $fn(msg...)
            Expr(:call, esc($fn), $lvar, msg...)
        end))
    end
    :(const $(esc(lvar)) = $rval)
end

macro register(lvar)
    mod = current_module()
    for fn in [:trace, :debug, :info, :warn, :error, :fatal]
        eval(mod, :( macro $fn(msg...)
            Expr(:call, esc($fn), $lvar, msg...)
        end))
    end
    :()
end

"Return logging context by name"
function getcontext(ctxname::AbstractString;
                    cfgloc::Nullable{AbstractString}=Nullable{AbstractString}(),
                    cfgexp::Nullable{Expr}=Nullable{Expr}())

    debug(LOGGER, """Configuration file is $(get(cfgloc, "not provided"))""")

    # get context (and set configuration location if available)
    ctx = context(LOG4JL_CONTEXT_SELECTOR, ctxname, get(cfgloc, ""))

    # start context if necessary
    if ctx.state == LifeCycle.INITIALIZED
        # configuration defined as expression
        if !isnull(cfgexp)
            debug(LOGGER, "Evaluating configuration block")
            cfg = evalConfiguration(cfgexp)
            start(ctx, cfg)
        elseif !isnull(cfgloc)
            cfg = getconfig(get(cfgloc), ctxname)
            start(ctx, cfg)
        else
            start(ctx)
        end
    end

    return ctx
end

"Shutdown logging system by stopping all logger contexts"
function shutdown()
    for (ctxname, ctx) in Log4jl.contexts(Log4jl.LOG4JL_CONTEXT_SELECTOR)
        stop(ctx)
    end
end

function __init__()
    # Default line separator
    global const LOG4JL_LINE_SEPARATOR = "LOG4JL_LINE_SEPARATOR" in keys(ENV) ?
                                         convert(Vector{UInt8}, ENV["LOG4JL_LINE_SEPARATOR"]) :
                                         @windows? [0x0d, 0x0a] : [0x0a]
    eval(Layouts, parse("import ..Log4jl.LOG4JL_LINE_SEPARATOR"))

    # Default logger status level
    global const LOG4JL_DEFAULT_STATUS_LEVEL = "LOG4JL_DEFAULT_STATUS_LEVEL" in keys(ENV) ?
                                               get(getlevel(ENV["LOG4JL_DEFAULT_STATUS_LEVEL"]), Level.ERROR) :
                                               Level.ERROR

    # Default internal logger status level
    global const LOG4JL_INTERNAL_STATUS_LEVEL = "LOG4JL_INTERNAL_STATUS_LEVEL" in keys(ENV) ?
                                               get(getlevel(ENV["LOG4JL_INTERNAL_STATUS_LEVEL"]), Level.WARN) :
                                               Level.WARN

    # Default logger event generator
    global const LOG4JL_LOG_EVENT = "LOG4JL_LOG_EVENT" in keys(ENV) ?
                                    eval(parse(ENV["LOG4JL_LOG_EVENT"])) :
                                    Log4jlEvent

    # Default logger context selector
    global const LOG4JL_CONTEXT_SELECTOR = "LOG4JL_CONTEXT_SELECTOR" in keys(ENV) ?
                                           eval(Expr(:call, parse(ENV["LOG4JL_CONTEXT_SELECTOR"]))) :
                                           ModuleContextSelector()

    # Internal status logger
    global const LOGGER = SimpleLogger("StatusLogger", LOG4JL_DEFAULT_MESSAGE, LOG4JL_INTERNAL_STATUS_LEVEL)
    LOGGER.showname = false # do not show
end

end


text:
@testset "Checking converging data (Real)" begin
    # block of maximally correlated values:
    N_corr = 16

    # number of blocks
    N_blocks = 131_072 # 2^17

    uncorrelated = rand(N_blocks) # rand(rng, N_blocks)
    av = mean(uncorrelated)
    _var = var(uncorrelated)
    stderr = sqrt(_var/N_blocks)

    correlated = [x for x in uncorrelated for _ in 1:N_corr]
    @test BinningAnalysis.unbinned_tau(correlated) ≈ 7.5 rtol = 0.05

    @testset "FullBinner" begin
        B = FullBinner(correlated)

        # Strongly/Maximally correlated

        # Note:
        # With xs = [sample, sample] we have
        # 1/N sum(xs) = 1/N (sum(xs[1:N/2] + sum(xs[N/2+1:N])) = 2/N sum(xs[1:N/2])
        # thus both the mean and variance should not vary when repeating data
        # this doesn't seem to be exactly true for Statistics.var(vcat(sample, sample))
        # this isn't exactly true for BinningAnalysis due to truncation at odd binning levels

        # Expectations:
        # - mean matches
        # - var matches up to binsize = N_corr for power binsizes
        # - varN = var / N
        # - std_error = sqrt(var / N) grows with sqrt(N)
        # - autocorrelation time should be 0.5 * (binsize-1) (i.e. maximal)

        for binsize in 2 .^ (0:3)
            target = 0.5binsize - 0.5

            @test mean(B, binsize) ≈ av
            @test var(B, binsize) ≈ _var rtol=0.01
            @test varN(B, binsize) * N_blocks ≈ _var * binsize / N_corr rtol=0.01
            @test std_error(B, binsize) * sqrt(N_corr / binsize) ≈ stderr rtol=0.01
            @test tau(B, binsize) ≈ target rtol=0.01
        end

        # Not correlated

        # Note:
        # If we expand the variance of pair binned data
        # \frac{1}{N} \sum_{i=1}^N (\frac{x_{2i-1} + x_{2i}}{2} - \mu)^2
        # and pull out jalf the unbinned variance
        # \frac{1}{2} \frac{1}{2N} \sum_{i=1}^{2N} (x_i - \mu)^2
        # we are left with an additional term
        # \frac{1}{2} \frac{1}{N} \sum_{i=1}^N (x_{2i-1} x_{2i} - 2\mu (x_{2i-1} x_{2i}) + \mu^2)
        # if we assume no correlation we can simplify the first term as
        # \langle x_{2i-1} x_{2i} \rangle = \langle x_{2i-1} \rangle \langle x_{2i} \rangle = \mu^2
        # and the term drops, leaving us with
        # var(x_{pair binned}) ≈ \frac{1}{2} var(x)

        # At and after binning out correlations
        # - mean matches
        # - var ≈ unbinned var / 2
        # - varN = var / N with N = N_blocks * N_corr / binsize
        # - std_error = sqrt(var / N) should be constant as the change in N is
        #                             cancelled by the change in var
        # - tau should be roughly constant

        for binsize in 2 .^ (4:10)
            f = N_corr / binsize

            @test mean(B, binsize) ≈ av
            @test var(B, binsize) ≈ f * _var rtol=0.075
            @test varN(B, binsize) * N_blocks ≈ _var rtol=0.075
            @test std_error(B, binsize) ≈ stderr rtol=0.075
            @test tau(B, binsize) ≈ 0.5(N_corr - 1) rtol = 0.075
        end
    end

    @testset "LogBinner" begin
        B = LogBinner(correlated)

        # Strongly/Maximally correlated
        for lvl in 1:4
            binsize = 2^(lvl-1)
            target = 0.5 * (binsize - 1)

            @test mean(B, lvl) ≈ av
            @test var(B, lvl) ≈ _var rtol=0.01
            @test varN(B, lvl) * N_blocks ≈ _var * binsize / N_corr rtol=0.01
            @test std_error(B, lvl) * sqrt(N_corr / binsize) ≈ stderr rtol=0.01
            @test tau(B, lvl) ≈ target rtol=0.01
        end

        # Not correlated
        for lvl in 5:11
            binsize = 2.0 ^ (lvl-1)
            f = N_corr / binsize

            @test mean(B, lvl) ≈ av
            @test var(B, lvl) ≈ f * _var rtol=0.075
            @test varN(B, lvl) * N_blocks ≈ _var rtol=0.075
            @test std_error(B, lvl) ≈ stderr rtol=0.075
            @test tau(B, lvl) ≈ 0.5(N_corr - 1) rtol = 0.075
        end
    end

    @testset "ErrorPropagator" begin
        B = ErrorPropagator(correlated)

        # Strongly/Maximally correlated
        for lvl in 1:4
            binsize = 2^(lvl-1)
            target = 0.5 * (binsize - 1)

            @test mean(B, 1, lvl) ≈ av
            @test var(B, 1, lvl) ≈ _var rtol=0.01
            @test varN(B, 1, lvl) * N_blocks ≈ _var * binsize / N_corr rtol=0.01
            @test std_error(B, 1, lvl) * sqrt(N_corr / binsize) ≈ stderr rtol=0.01
            @test tau(B, 1, lvl) ≈ target rtol=0.01
        end

        # Not correlated
        for lvl in 5:11
            binsize = 2.0 ^ (lvl-1)
            f = N_corr / binsize

            @test mean(B, 1, lvl) ≈ av
            @test var(B, 1, lvl) ≈ f * _var rtol=0.075
            @test varN(B, 1, lvl) * N_blocks ≈ _var rtol=0.075
            @test std_error(B, 1, lvl) ≈ stderr rtol=0.075
            @test tau(B, 1, lvl) ≈ 0.5(N_corr - 1) rtol = 0.075
        end
    end
end

@testset "Ising Energies" begin
    # This data comes from an AFM Ising model on a triangular lattice, L = 8, T = 2.0
    Es = open("test.data", "r") do f
        Float64[read(f, Float64) for _ in 1:2^14]
    end

    av = mean(Es)
    τ = BinningAnalysis.unbinned_tau(Es)
    stderr = sqrt(var(Es) * (1 + 2τ) / length(Es))

    @testset "FullBinner" begin
        B = FullBinner(Es)
        N = BinningAnalysis._reliable_level(B)

        @test mean(B) ≈ av
        @test tau(B) ≈ τ rtol=0.5
        @test std_error(B) ≈ stderr rtol=0.3
    end

    @testset "LogBinner" begin
        B = LogBinner(Es)

        @test mean(B) ≈ av
        @test tau(B) ≈ τ rtol=0.5
        @test std_error(B) ≈ stderr rtol=0.3
    end

    @testset "ErrorPropagator" begin
        B = ErrorPropagator(Es)

        @test mean(B, 1) ≈ av
        @test tau(B, 1) ≈ τ rtol=0.5
        @test std_error(B, 1) ≈ stderr rtol=0.3
    end
end

max_stars_repo_name: crstnbr/BinningAnalysis

max_issues_repo_name: crstnbr/BinningAnalysis

max_forks_repo_name: crstnbr/BinningAnalysis


text:
function resolve_import(x, state::State)
    if typof(x) === Using || typof(x) === Import
        u = typof(x) === Using
        i = 2
        n = length(x)

        root = par = getsymbolserver(state.server)
        bindings = []
        while i <= n
            arg = x[i]
            if isidentifier(arg) || typof(arg) === CSTParser.MacroName
                if refof(arg) !== nothing
                    par = refof(arg)
                else
                    par = _get_field(par, arg, state)
                end
            elseif typof(arg) === PUNCTUATION && kindof(arg) == CSTParser.Tokens.COMMA
                # end of chain, make available
                if i > 2
                    _mark_import_arg(x[i - 1], par, state, u)
                end
                par = root
            elseif typof(arg) === OPERATOR && kindof(arg) == CSTParser.Tokens.COLON
                root = par
                if par !== nothing && i > 2 && isidentifier(x[i - 1]) && refof(x[i - 1]) === nothing
                    setref!(x[i - 1], par)
                end
            elseif typof(arg) === PUNCTUATION && kindof(arg) == CSTParser.Tokens.DOT
                # dot between identifiers
                if par !== nothing && i > 2 && isidentifier(x[i - 1]) && refof(x[i - 1]) === nothing
                    setref!(x[i - 1], par)
                end
            elseif typof(arg) === OPERATOR && kindof(arg) == CSTParser.Tokens.DOT
                # dot prexceding identifiser
                if par == root == getsymbolserver(state.server)
                    par = state.scope
                elseif par isa Scope && parentof(par) !== nothing
                    par = parentof(par)
                else
                    return
                end
            else
                return
            end
            if i == n
                _mark_import_arg(x[i], par, state, u)
            end
            i += 1
        end
    end
end

function _mark_import_arg(arg, par, state, u)
    if par !== nothing && (typof(arg) === IDENTIFIER || typof(arg) === MacroName)
        if par isa Binding # mark reference to binding
            push!(par.refs, arg)
        end
        if par isa SymbolServer.VarRef
            par = SymbolServer._lookup(par, getsymbolserver(state.server), true)
            !(par isa SymbolServer.SymStore) && return
        end
        if bindingof(arg) === nothing
            if !hasmeta(arg)
                arg.meta = Meta()
            end
            arg.meta.binding = Binding(arg, par, _typeof(par, state), [], nothing, nothing)
        end
        if u && par isa SymbolServer.ModuleStore
            if state.scope.modules isa Dict
                state.scope.modules[Symbol(valof(arg))] = par
            else
                state.scope.modules = Dict(Symbol(valof(arg)) => par)
            end
        elseif u && par isa Binding && par.val isa SymbolServer.ModuleStore
            if state.scope.modules isa Dict
                state.scope.modules[Symbol(valof(arg))] = par.val
            else
                state.scope.modules = Dict(Symbol(valof(arg)) => par.val)
            end
        elseif u && par isa Binding && par.val isa EXPR && (typof(par.val) === CSTParser.ModuleH || typof(par.val) === CSTParser.BareModule)
            if state.scope.modules isa Dict
                state.scope.modules[Symbol(valof(arg))] = scopeof(par.val)
            else
                state.scope.modules = Dict(Symbol(valof(arg)) => scopeof(par.val))
            end
        elseif u && par isa Binding && par.val isa Binding && par.val.val isa EXPR && (typof(par.val.val) === CSTParser.ModuleH || typof(par.val.val) === CSTParser.BareModule)
            if state.scope.modules isa Dict
                state.scope.modules[Symbol(valof(arg))] = scopeof(par.val.val)
            else
                state.scope.modules = Dict(Symbol(valof(arg)) => scopeof(par.val.val))
            end
        end
    end
end

function _get_field(par, arg, state)
    arg_str_rep = CSTParser.str_value(arg)
    if par isa SymbolServer.EnvStore && haskey(par, Symbol(arg_str_rep))
        return par[Symbol(arg_str_rep)]
    elseif par isa SymbolServer.ModuleStore # imported module
        if Symbol(arg_str_rep) === par.name.name
            return par
        elseif haskey(par, Symbol(arg_str_rep))
            par = par[Symbol(arg_str_rep)]
            if par isa SymbolServer.VarRef # reference to dependency
                return SymbolServer._lookup(par, getsymbolserver(state.server), true)
            end
            return par
        end
    elseif par isa Scope && scopehasbinding(par, arg_str_rep)
        return par.names[arg_str_rep]
    elseif par isa Binding
        if par.val isa Binding
            return _get_field(par.val, arg, state)
        elseif par.val isa EXPR && CSTParser.defines_module(par.val) && scopeof(par.val) isa Scope
            return _get_field(scopeof(par.val), arg, state)
        elseif par.val isa SymbolServer.ModuleStore
            return _get_field(par.val, arg, state)
        end
    end
    return
end


text:
module PkgSegments

import TOML
import UUIDs

"""
    PackageKey

A helper struct for capturing and parsing bacic information
about a package. Either using the naming convention of
PkgName or PkgName:UUID to uniquely identify the package within
the project
"""
struct PackageKey
    name::String
    uuid::Union{Nothing, UUIDs.UUID}

    function PackageKey(name::AbstractString)
        parts = split(name, ":")
        if length(parts) == 1
            PackageKey(parts[1], nothing)
        elseif length(parts) == 2
            PackageKey(parts[1], parts[2])
        else
            throw("Unable to construct PackageKey from $name")
        end
    end
    PackageKey(name::AbstractString, uuid::Nothing) = new(name, nothing)
    PackageKey(name::AbstractString, uuid::AbstractString) = new(name, UUIDs.UUID(uuid))
    PackageKey(name::AbstractString, uuid::UUIDs.UUID) = new(name, uuid)
    PackageKey(params::Pair{String, <:Any}) = PackageKey(params.first, params.second)
end

"""
    Base.:(==)(k1::PackageKey, k2::PackageKey)

Overrides equality of the PackageKey so that if the UUID is not defined
then the name is used solely for identifying equality.
"""
function Base.:(==)(k1::PackageKey, k2::PackageKey)
    if k1.uuid === nothing || k2.uuid === nothing
        k1.name == k2.name
    else
        k1.name == k2.name && k1.uuid == k2.uuid
    end
end

"""
    Base.hash(k::PackageKey, h::UInt64)

Override the hash inline with the equality function, so that we just
generate the hash from the name.
"""
Base.hash(k::PackageKey, h::UInt64) = hash(k.name, h)

"""
    generatesegment!(directory::AbstractString, deps::Set{PackageKey}; subdir::AbstractString="seg")

Generate the segmented Manifest.toml and Project.toml in given directory
based on the desired dependencies.
"""
function generatesegment!(directory::AbstractString, deps::Set{PackageKey}; subdir::AbstractString="seg")
    mkpath(joinpath(directory, subdir))
    open(joinpath(directory, subdir, "Project.toml"), "w") do io
        TOML.print(io, projectsegment(directory, deps); sorted=true)
    end

    open(joinpath(directory, subdir, "Manifest.toml"), "w") do io
        TOML.print(io, manifestsegment(directory, deps); sorted=true)
    end
end

"""
    genfromsegfile!(directory::AbstractString)

Uses a PkgSegments.toml file in the directory to specify the segments
to generate. Each section should specify the `deps` as an array and the
subdir to create the files as `subdir`
"""
function genfromsegfile!(directory::AbstractString)
    segments = TOML.parsefile(joinpath(directory, "PkgSegments.toml"))
    for (name, segment) in segments
        deps = Set([PackageKey(k) for k in segment["deps"]])
        generatesegment!(directory, deps; subdir=segment["subdir"])
    end
end

manifesttokey(data::Dict{String,<:Any}) = Dict([PackageKey(n,v["uuid"]) => v for (n,vs) in data for v in vs])
keytomanifest(data::Dict{PackageKey,<:Any}) = reduce(data; init=Dict{String,Any}()) do val, (key, d)
    entries = get!(val, key.name, [])
    push!(entries, d)
    val
end

manifestkey(entry::Pair, manifest) = PackageKey(entry)
manifestkey(entry, manifest) = PackageKey(entry, manifest[PackageKey(entry)]["uuid"])

function segmentdata(project::Dict{String,<:Any}, manifest::Dict{String,<:Any}, deps::Set{PackageKey})
    keymanifest = manifesttokey(manifest)
    newdeps = Set{PackageKey}()

    depqueue = Vector{PackageKey}()
    for dep in deps
        if dep ∈ keys(keymanifest)
            push!(newdeps, PackageKey(dep.name, keymanifest[dep]["uuid"]))
            push!(depqueue, dep)
        end
    end

    while length(depqueue) > 0
        key = pop!(depqueue)
        entry = keymanifest[key]

        entrydeps = [manifestkey(d, keymanifest) for d in get(entry, "deps", [])]

        append!(depqueue, entrydeps)
        newdeps = ∪(newdeps, entrydeps)
    end

    for key in keys(keymanifest)
        if key ∉ newdeps
            delete!(keymanifest, key)
        end
    end

    segmentmanifest = keytomanifest(keymanifest)
    # Now let's get the segment project. We only maintain a specific set of fields
    segmentproject = Dict{String,Any}()
    for field in ["name", "uuid", "deps"]
        segmentproject[field] = project[field]
    end
    currentdeps = Set([PackageKey(key) for key in segmentproject["deps"]])
    for dep in newdeps
        if dep ∉ currentdeps
            segmentproject["deps"][dep.name] = string(dep.uuid)
        end
    end

    for entry in segmentproject["deps"]
        key = PackageKey(entry)
        if key ∉ newdeps
            delete!(segmentproject["deps"], entry[1])
        end
    end

    (segmentmanifest, segmentproject)
end

function projectsegment(directory::AbstractString, deps::Set{PackageKey})
    alldeps = copy(deps)
    # Make sure we include julia, which is pretty much a guarantee.
    push!(alldeps, PackageKey("julia"))
    data = TOML.parsefile(joinpath(directory, "Project.toml"))
    # Now we work through the dependencies and associated compats.
    for entry in data["deps"]
        key = PackageKey(entry)
        if key ∉ deps
            delete!(data["deps"], entry[1])
            if entry[1] in keys(data["compat"])
                delete!(data["compat"], entry[1])
            end
        end
    end

    data
end

function manifestsegment(directory::AbstractString, deps::Set{PackageKey})
    alldeps = copy(deps)
    data = TOML.parsefile(joinpath(directory, "Manifest.toml"))

    depqueue = collect(alldeps)
    while length(depqueue) > 0
        key = pop!(depqueue)
        entries = data[key.name]
        if key.uuid !== nothing
            entries = [e for e in entries if UUIDs.UUID(e["uuid"]) == key.uuid]
        end
        length(entries) != 1 && error("Unable to identify unique package for $key")
        entry = entries[1]
        deps = [PackageKey(d) for d in get(entry, "deps", [])]
        append!(depqueue, deps)
        for dep in deps
            push!(alldeps, dep)
        end
    end

    for (name, info) in data
        subentries = [s for s in info if PackageKey(name, s["uuid"]) ∈ alldeps]
        if length(subentries) > 0
            data[name] = subentries
        else
            delete!(data, name)
        end
    end

    data
end

end # module


text:
@testset "exports" begin
    # Ensure that every exported name is actually defined
    for name in names(PressureFieldContact)
        @test isdefined(PressureFieldContact, name)
        !isdefined(PressureFieldContact, name) && println("MISSING: ", name)
    end
end


text:
using Test
using ChemistryFeaturization.Utils.GraphBuilding
using Xtals

@testset "GraphBuilding" begin
    path1 = abspath(@__DIR__, "..", "test_data", "strucs", "mp-195.cif")
    adj, els = build_graph(path1; use_voronoi = true)
    wm_true = [0.0 1.0 1.0 1.0; 1.0 0.0 1.0 1.0; 1.0 1.0 0.0 1.0; 1.0 1.0 1.0 0.0]
    els_true = ["Ho", "Pt", "Pt", "Pt"]

    @test adj == wm_true
    @test els == els_true

    adj, els = build_graph(path1; use_voronoi = false)
    @test adj == wm_true
    @test els == els_true

    # test that we get the same results building from a Crsytal object
    adjc, elsc = build_graph(Crystal(path1))
    @test adjc == wm_true
    @test elsc == els_true
end


text:
# ---
# title: 214. Shortest Palindrome
# id: problem214
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String
# link: <https://leetcode.com/problems/shortest-palindrome/description/>
# hidden: true
# ---
#
# Given a string _**s**_ , you can convert it to a palindrome by adding
# characters in front of it. Find and return the shortest palindrome you can
# find by performing this transformation.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aacecaaa"
#     Output: "aaacecaaa"
#
#
# **Example 2:**
#
#
#
#     Input: s = "abcd"
#     Output: "dcbabcd"
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 5 * 104`
#   * `s` consists of lowercase English letters only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=
Template Method allows an algorithm to be built from different parts
where subclasses can provide implementations without changing the
structure of the algorithm.
=#

# Julia's first class functions can be used for that purpose.
module MLPipeline

using GLM, DataFrames, LinearAlgebra

"""
Run machine learning pipeline
"""
function run(data::DataFrame, response::Symbol, predictors::Vector{Symbol})
    train, test = split_data(data, 0.7)
    model = fit(train, response, predictors)
    validate(test, model, response)
end

"""
Split data randomly and return a tuple of data frames for training and testing
respectively.  Specify `keep` as a number between 0 and 1.0 for the percentage
of records kept for testing.
"""
split_data(df::DataFrame, keep::Float64) = let b = rand(size(df,1)) .< keep
    (df[b, :], df[.!b, :])
end

"""
Return a model that can be used for predicting the `response` variable
from the `predictors` variables.
"""
function fit(df::DataFrame, response::Symbol, predictors::Vector{Symbol})
    formula = Term(response) ~ +(Term.(predictors)...)
    return lm(formula, df)
end

"""
Predict response variable using the model and validate accuracy of model.
"""
function validate(df, model, response)
    yhat = predict(model, df)
    y = df[:, response]
    return (result = [y yhat], rmse = rmse(yhat .- y))
end

"Root mean square error."
rmse(ys) = norm(ys) / sqrt(length(ys))

"""
Run machine learning pipeline (template version w/ keyword args)
"""
function run2(data::DataFrame, response::Symbol, predictors::Vector{Symbol};
            fit = fit, split_data = split_data, validate = validate)
    train, test = split_data(data, 0.7)
    model = fit(train, response, predictors)
    validate(test, model, response)
end

end #module

using RDatasets, GLM
using .MLPipeline

# Using original template
boston = dataset("MASS", "Boston");
result, rmse = MLPipeline.run(boston, :MedV, [:Rm, :Tax, :Crim]);
println(rmse)

# Custom fitting function
function fit_glm(df::DataFrame, response::Symbol, predictors::Vector{Symbol})
    formula = Term(response) ~ +(Term.(predictors)...)
    return glm(formula, df, Normal(), IdentityLink())
end

result, rmse = MLPipeline.run2(
    boston, :MedV, [:Rm, :Tax, :Crim],
    fit = fit_glm);

println(rmse)