text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module TerminalMenus

terminal = nothing  # The user terminal

import REPL

function __init__()
    global terminal
    terminal = REPL.Terminals.TTYTerminal(get(ENV, "TERM", Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)
end

include("util.jl")
include("config.jl")

include("AbstractMenu.jl")
include("RadioMenu.jl")
include("MultiSelectMenu.jl")

export
    RadioMenu,
    MultiSelectMenu,
    request

# TODO: remove in Julia 2.0
# While not exported, AbstractMenu documented these as an extension interface
@deprecate printMenu printmenu
function writeLine(buf::IOBuffer, m::AbstractMenu, idx::Int, cursor::Bool)
    Base.depwarn("`writeLine` is deprecated, use `writeline` instead.", :writeLine)
    error("unimplemented")
end

end # module

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using WNNMDenoise
using MAT, FileIO
using ImageQualityIndexes

@assert VERSION >= v"1.7.0-beta4"

using LinearAlgebra
using MKL
@info BLAS.get_config() # MKL is faster on svd

for σₙ = [20, 20, 40, 60, 80, 100]
    # we run 20 twice to get rid of the compilation time
    img = Float32.(load("house.png")) * 255;
    noisy_img = matopen(joinpath(@__DIR__, "data", "house_AWGN_$(σₙ).mat")) do io
        read(io, "N_Img")
    end

    start_psnr = assess_psnr(noisy_img, img, 255)

    out = copy(noisy_img)
    f = WNNM(σₙ)
    @time f(out, noisy_img, clean_img=img);

    end_psnr = assess_psnr(out, img, 255)
end

ext: jl


text:
#include("../src/ColonyRuntime.jl")
#include("../src/Crypto.jl")

#import .ColonyRuntime
#import .Crypto

using Pkg
using Test

Pkg.activate("..")

using Test
using ColonyRuntime
using ColonyRuntime.Crypto

Base.exit_on_sigint(false)

global server_prvkey = "09545df1812e252a2a853cca29d7eace4a3fe2baad334e3b7141a98d43c31e7b"

server = ColonyRuntime.ColoniesServer("localhost", 8080)

function test_addcolony()
  colonyid = Crypto.id(Crypto.prvkey())
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  added_colony == colony
end

function test_addruntime()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid, "amd", 1, 1024, "nvidia", 1, ColonyRuntime.PENDING)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)

  added_runtime == runtime
end

function test_approveruntime()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid, "amd", 1, 1024, "nvidia", 1, ColonyRuntime.PENDING)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)
  true # an exception will be raised if it does not work
end

function test_submitprocess()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid, "amd", 1, 1024, "nvidia", 1, ColonyRuntime.PENDING)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)

  conditions = ColonyRuntime.Conditions(colonyid, [], "test_runtime_type", 1, 1024, 1)
  env = Dict()
  env["args"] = "test_args"
  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  added_process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)

  # we don't bother testing all process attrbiutes, at least the process was assing a 64 chars process id
  length(added_process.processid) == 64
end

function test_getprocess()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid, "amd", 1, 1024, "nvidia", 1, ColonyRuntime.PENDING)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)

  conditions = ColonyRuntime.Conditions(colonyid, [], "test_runtime_type", 1, 1024, 1)
  env = Dict()
  env["args"] = "test_args"
  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  added_process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)

  process_from_server = ColonyRuntime.getprocess(server, added_process.processid, runtime_prvkey)
  process_from_server.processid == added_process.processid
end

function test_assignprocess()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid, "amd", 1, 1024, "nvidia", 1, ColonyRuntime.PENDING)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)

  conditions = ColonyRuntime.Conditions(colonyid, [], "test_runtime_type", 1, 1024, 1)
  env = Dict()
  env["args"] = "test_args"
  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  added_process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)

  assigned_process = ColonyRuntime.assignprocess(server, colonyid, runtime_prvkey)
  assigned_process.processid == added_process.processid
end

function test_addattribute()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)

  conditions = ColonyRuntime.Conditions(colonyid, [], "test_runtime_type", 1, 1024, 1)
  env = Dict()
  env["args"] = "test_args"
  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  added_process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)

  assigned_process = ColonyRuntime.assignprocess(server, colonyid, runtime_prvkey)
  assigned_process.processid == added_process.processid

  attribute = ColonyRuntime.Attribute(assigned_process.processid, "test_result", "test_result_value")
  added_attribute = ColonyRuntime.addattribute(server, attribute, runtime_prvkey)

  added_attribute.key == "test_result"
end

function test_closeprocess_successful()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)

  conditions = ColonyRuntime.Conditions(colonyid, [], "test_runtime_type", 1, 1024, 1)
  env = Dict()
  env["args"] = "test_args"
  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  added_process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)

  assigned_process = ColonyRuntime.assignprocess(server, colonyid, runtime_prvkey)
  assigned_process.processid == added_process.processid

  attribute = ColonyRuntime.Attribute(assigned_process.processid, "test_result", "test_result_value")
  added_attribute = ColonyRuntime.addattribute(server, attribute, runtime_prvkey)

  ColonyRuntime.closeprocess(server, assigned_process.processid, true, runtime_prvkey)

  process_from_server = ColonyRuntime.getprocess(server, assigned_process.processid, runtime_prvkey)
  process_from_server.state == ColonyRuntime.SUCCESS
end

function test_closeprocess_failed()
  colony_prvkey = Crypto.prvkey()
  colonyid = Crypto.id(colony_prvkey)
  colony = ColonyRuntime.Colony(colonyid, "test_name2")
  added_colony = ColonyRuntime.addcolony(server, colony, server_prvkey)

  runtime_prvkey = Crypto.prvkey()
  runtime = ColonyRuntime.Runtime(Crypto.id(runtime_prvkey), "test_runtime_type", "test_runtime_name", colonyid)
  added_runtime = ColonyRuntime.addruntime(server, runtime, colony_prvkey)
  ColonyRuntime.approveruntime(server, added_runtime.runtimeid, colony_prvkey)

  conditions = ColonyRuntime.Conditions(colonyid, [], "test_runtime_type", 1, 1024, 1)
  env = Dict()
  env["args"] = "test_args"
  processpec = ColonyRuntime.ProcessSpec(-1, -1, conditions, env)
  added_process = ColonyRuntime.submitprocess(server, processpec, runtime_prvkey)

  assigned_process = ColonyRuntime.assignprocess(server, colonyid, runtime_prvkey)
  assigned_process.processid == added_process.processid

  attribute = ColonyRuntime.Attribute(assigned_process.processid, "test_result", "test_result_value")
  added_attribute = ColonyRuntime.addattribute(server, attribute, runtime_prvkey)

  ColonyRuntime.closeprocess(server, assigned_process.processid, false, runtime_prvkey)

  closed_process = ColonyRuntime.getprocess(server, assigned_process.processid, runtime_prvkey)
  closed_process.state == ColonyRuntime.FAILED
end

@testset begin
try
  @test test_addcolony()
  @test test_addruntime()
  @test test_approveruntime()
  @test test_submitprocess()
  @test test_getprocess()
  @test test_assignprocess()
  @test test_addattribute()
  @test test_closeprocess_successful()
  @test test_closeprocess_failed()
catch err
  typeof(err) == InterruptException && rethrow(err)
  print(err)
end

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const _DataType = (;[
    Symbol("DT_INVALID") => Int32(0),
    Symbol("DT_FLOAT") => Int32(1),
    Symbol("DT_DOUBLE") => Int32(2),
    Symbol("DT_INT32") => Int32(3),
    Symbol("DT_UINT8") => Int32(4),
    Symbol("DT_INT16") => Int32(5),
    Symbol("DT_INT8") => Int32(6),
    Symbol("DT_STRING") => Int32(7),
    Symbol("DT_COMPLEX64") => Int32(8),
    Symbol("DT_INT64") => Int32(9),
    Symbol("DT_BOOL") => Int32(10),
    Symbol("DT_QINT8") => Int32(11),
    Symbol("DT_QUINT8") => Int32(12),
    Symbol("DT_QINT32") => Int32(13),
    Symbol("DT_BFLOAT16") => Int32(14),
    Symbol("DT_QINT16") => Int32(15),
    Symbol("DT_QUINT16") => Int32(16),
    Symbol("DT_UINT16") => Int32(17),
    Symbol("DT_COMPLEX128") => Int32(18),
    Symbol("DT_HALF") => Int32(19),
    Symbol("DT_RESOURCE") => Int32(20),
    Symbol("DT_VARIANT") => Int32(21),
    Symbol("DT_UINT32") => Int32(22),
    Symbol("DT_UINT64") => Int32(23),
    Symbol("DT_FLOAT_REF") => Int32(101),
    Symbol("DT_DOUBLE_REF") => Int32(102),
    Symbol("DT_INT32_REF") => Int32(103),
    Symbol("DT_UINT8_REF") => Int32(104),
    Symbol("DT_INT16_REF") => Int32(105),
    Symbol("DT_INT8_REF") => Int32(106),
    Symbol("DT_STRING_REF") => Int32(107),
    Symbol("DT_COMPLEX64_REF") => Int32(108),
    Symbol("DT_INT64_REF") => Int32(109),
    Symbol("DT_BOOL_REF") => Int32(110),
    Symbol("DT_QINT8_REF") => Int32(111),
    Symbol("DT_QUINT8_REF") => Int32(112),
    Symbol("DT_QINT32_REF") => Int32(113),
    Symbol("DT_BFLOAT16_REF") => Int32(114),
    Symbol("DT_QINT16_REF") => Int32(115),
    Symbol("DT_QUINT16_REF") => Int32(116),
    Symbol("DT_UINT16_REF") => Int32(117),
    Symbol("DT_COMPLEX128_REF") => Int32(118),
    Symbol("DT_HALF_REF") => Int32(119),
    Symbol("DT_RESOURCE_REF") => Int32(120),
    Symbol("DT_VARIANT_REF") => Int32(121),
    Symbol("DT_UINT32_REF") => Int32(122),
    Symbol("DT_UINT64_REF") => Int32(123),
]...)

const SpecializedType = (;[
    Symbol("ST_INVALID") => Int32(0),
    Symbol("ST_TENSOR_LIST") => Int32(1),
]...)

export _DataType, SpecializedType


text:
# Comparison between the Short Pulse Aproximation and Pulse Gradient Spin Echo (with Laplacian Eigen Functions)

using MRIsim, Plots
using MRIsim: γ, Planes, SignalE, DIF_base, plot_grads

#Constants
𝒊 = 1im; D = 2e-9 # m2/s
#DIF Sequence
T = 80e-3; δ = 10e-3;
DIF = DIF_base(1,T-δ,δ); #plot_grads(DIF)
#Restricted planes with Short Pulse Approximation (SPA) - from my thesis
Erest(L,q) = sum([(n==0 ? 1 : 2)*exp(-D*(T-δ)*(n*π/L)^2)*sinc.(L*q .+ n/2).^2 ./(1 .- n./(2*L*q)).^2 for n = 0:100])
#Free diffusion - SPA
Efree(q) = exp.(-4π^2*D*q^2)

## Comparison
gr()
plot()
g = (0:1000)*1e-3
q = γ*g*δ #q-values
color = [:blue, :red, :green]
for (i, L) = enumerate([2, 7, 90]*1e-6)
    Ll = round(Int,L*1e6)
    μ = Planes(L,D)
    #Laplacian Eigen Functions for different gradient strengths
    E = [SignalE(μ, gi*DIF) for gi ∈ g]
    plot!(q*1e-6, abs.(E), linecolor=color[i], label="L = "*string(Ll)*" μm", yaxis=:log, legend=:outertopright)
    plot!(q*1e-6, abs.(Erest.(L,q)), linecolor=color[i], label="L = "*string(Ll)*" μm", linestyle=:dash,yaxis=:log, legend=:outertopright)
end
plot!(q*1e-6, abs.(Efree.(q)),linecolor=:purple,label="L = ∞ μm",yaxis=:log, legend=:outertopright)
ylims!(1e-6,1)
ylabel!("E(q)")
xlabel!("q [μm⁻¹]")
title!("Laplacian Eigen Functions vs Short Pulse Approximation")
#savefig("./src/examples/Figures/PGSEvsSPA.pdf")


text:
mutable struct RingPolymerArray{T} <: AbstractArray{T,3}
    data::Array{T,3}
    quantum_atoms::Vector{Int}
    classical_atoms::Vector{Int}
    normal::Bool
    function RingPolymerArray{T}(data::Array{T,3}, quantum, classical, normal) where {T}
        constrain_classical_atoms!(data, classical)
        new(data, quantum, classical, normal)
    end
end

function RingPolymerArray(data::Array{T,3}; quantum=nothing, normal=false) where {T}
    quantum === nothing && (quantum = collect(axes(data,2)))
    classical = collect(axes(data,2))
    setdiff!(classical, quantum)
    RingPolymerArray{T}(data, quantum, classical, normal)
end

Base.size(A::RingPolymerArray) = size(A.data)
Base.getindex(A::RingPolymerArray, i::Int, j::Int, k::Int) = A.data[i,j,k]

function Base.setindex!(A::RingPolymerArray, v, i::Int, j::Int, k::Int)
    if j in A.classical_atoms
        if k == 1
            setindex!(A.data, v, i, j, k)
            set_classical!(A, i, j)
        end
    else
        setindex!(A.data, v, i, j, k)
    end
end

function Base.similar(A::RingPolymerArray, ::Type{S}, dims::Dims) where {S}
    if length(dims) == 3
        return RingPolymerArray{S}(Array{S}(undef, dims), A.quantum_atoms, A.classical_atoms, A.normal)
    else
        return similar(A.data, S, dims)
    end
end

struct RingPolymerStyle <: Broadcast.AbstractArrayStyle{3} end
Base.BroadcastStyle(::Type{<:RingPolymerArray}) = RingPolymerStyle()
RingPolymerStyle(::Val{3}) = RingPolymerStyle()
RingPolymerStyle(::Val{N}) where {N} = Broadcast.DefaultArrayStyle{N}()

function Base.similar(bc::Broadcast.Broadcasted{RingPolymerStyle}, ::Type{Eltype}) where {Eltype}
    A = find_rpa(bc)
    if A === nothing
        return similar(Array{Eltype}, axes(bc))
    else
        return similar(A)
    end
end

"`A = find_rpa(As)` return the first RingPolymerArray among the arguments."
find_rpa(bc::Base.Broadcast.Broadcasted) = find_rpa(bc.args)
find_rpa(args::Tuple) = find_rpa(find_rpa(args[1]), Base.tail(args))
find_rpa(x) = x
find_rpa(::Tuple{}) = nothing
find_rpa(a::RingPolymerArray, rest) = a
find_rpa(::Any, rest) = find_rpa(rest)

"""
    constrain_classical_atoms!(A::AbstractArray{T,3}, classical_atoms) where {T}

Make sure that the classical atoms all have the same values.
"""
function constrain_classical_atoms!(A::AbstractArray{T,3}, classical_atoms) where {T}
    for j in classical_atoms
        for i in axes(A,1)
            set_classical!(A, i, j)
        end
    end
end

"Copy data from first index into replicas. All classical atoms have repeated values."
set_classical!(A::AbstractArray{T,3}, i, j) where {T} = A[i,j,2:end] .= A[i,j,1]
set_classical!(A::RingPolymerArray, i, j) = set_classical!(A.data, i, j)

"""
Evaluate centroid of ring polymer.
"""
function get_centroid(A::AbstractArray{T,3}) where {T}
    out = zeros(T, size(A)[1:2])
    get_centroid!(out, A)
    return out
end

function get_centroid!(out::Matrix{T}, A::RingPolymerArray{T}) where {T}
    if A.normal
        copyto!(out, @view A[:,:,1])
        out[:,A.quantum_atoms] ./= sqrt(size(A,3))
    else
        get_centroid!(out, A.data)
    end
    return out
end

function get_centroid!(out::Matrix{T}, A::AbstractArray{T,3}) where {T}
    fill!(out, zero(eltype(A)))
    for i in axes(A,3)
        for j in axes(A,2)
            for k in axes(A,1)
                @inbounds out[k,j] += A[k,j,i]
            end
        end
    end
    for I in eachindex(out)
        @inbounds out[I] /= size(A,3)
    end
    return out
end


text:
#==========================================================================================#

# ONE SAMPLE: INTERFACE

function hausman_1s(obj₁::ParModel, obj₂::ParModel)
    return hausman_1s(obj₁, obj₂, intersect(coefnames(obj₁), coefnames(obj₂)))
end

function hausman_1s(obj₁::ParModel, obj₂::ParModel, name::String)

    return hausman_1s(obj₁, obj₂, [name])
end

function hausman_1s(obj₁::ParModel, obj₂::ParModel, names::Vector{String})

    i₁    = findall((in)(names), coefnames(obj₁))
    i₂    = findall((in)(names), coefnames(obj₂))
    w₁    = getweights(obj₁)
    w₂    = getweights(obj₂)
    corr₁ = getcorr(obj₁)
    corr₂ = getcorr(obj₂)

    (iszero(i₁) | iszero(i₂)) && throw("missing coefficients in at least one model")
    isequal(corr₁, corr₂) || throw("different correlation structures")
    isequal(getweights(obj₁), getweights(obj₂)) || throw("different weighting schemes")

    output       = ParEstimate()
    output.names = copy(names)

    _hausman_1s!(output, obj₁, i₁, obj₂, i₂, w₁, corr₁)

    return output
end

# ONE SAMPLE: ESTIMATION

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::UnitWeights,
        corr::Heteroscedastic
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::FrequencyWeights,
        corr::Heteroscedastic
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::AbstractWeights,
        corr::Heteroscedastic
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁) ; lmul!(Diagonal(w), ψ₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::UnitWeights,
        corr::Clustered
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = corr.mat * view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = corr.mat * view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::AbstractWeights,
        corr::Clustered
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = corr.mat * view(Ψ₁, :, i₁) ; lmul!(Diagonal(w), ψ₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = corr.mat * view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::UnitWeights,
        corr::CrossCorrelated
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * corr.mat * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_1s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        obj₂::ParModel,
        i₂::Vector{Int},
        w::AbstractWeights,
        corr::CrossCorrelated
    )

    Ψ₁  = influence(obj₁, w)
    ψ₁  = view(Ψ₁, :, i₁) ; lmul!(Diagonal(w), ψ₁)
    Ψ₂  = influence(obj₂, w)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w), ψ₂)
    V₁₂ = ψ₁' * corr.mat * ψ₂

    adjfactor!(V₁₂, obj₁, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

#==========================================================================================#

# TWO INDEPENDENT SAMPLES

function hausman_2s(obj₁::ParModel, obj₂::ParModel)
    return hausman_2s(obj₁, obj₂, intersect(coefnames(obj₁), coefnames(obj₂)))
end

function hausman_2s(obj₁::ParModel, obj₂::ParModel, name::String)
    return hausman_2s(obj₁, obj₂, [name])
end

function hausman_2s(obj₁::ParModel, obj₂::ParModel, names::Vector{String})

    i₁ = findall((in)(names), coefnames(obj₁))
    i₂ = findall((in)(names), coefnames(obj₂))

    (iszero(i₁) | iszero(i₂)) && throw("missing coefficients in at least one model")

    output       = ParEstimate()
    output.names = copy(names)
    output.β     = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V     = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂)

    return output
end

#==========================================================================================#

# TWO DEPENDENT SAMPLES: INTERFACE

function hausman_2s(obj₁::ParModel, obj₂::ParModel, corr::CorrStructure)
    return hausman_2s(obj₁, obj₂, corr, intersect(coefnames(obj₁), coefnames(obj₂)))
end

function hausman_2s(obj₁::ParModel, obj₂::ParModel, corr::CorrStructure, name::String)
    return hausman_2s(obj₁, obj₂, corr, [name])
end

function hausman_2s(obj₁::ParModel, obj₂::ParModel, corr::CorrStructure, names::Vector{String})

    i₁ = findall((in)(names), coefnames(obj₁))
    i₂ = findall((in)(names), coefnames(obj₂))
    w₁ = getweights(obj₁)
    w₂ = getweights(obj₂)

    (iszero(i₁) | iszero(i₂)) && throw("missing coefficients in at least one model")

    if !(typeof(corr) == typeof(getcorr(obj₁)) == typeof(getcorr(obj₂)))
        throw("incompatible correlation structures")
    end

    output       = ParEstimate()
    output.names = copy(names)

    _hausman_2s!(output, obj₁, i₁, w₁, obj₂, i₂, w₂, corr)

    return output
end

# TWO DEPENDENT SAMPLES: ESTIMATION

function _hausman_2s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::UnitWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::UnitWeights,
        corr::Heteroscedastic
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)
    touse       = nonmissing₁ .* nonmissing₂
    touse₁      = touse[nonmissing₁]
    touse₂      = touse[nonmissing₂]

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, touse₁, i₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, touse₂, i₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, obj₂, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_2s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::AbstractWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::AbstractWeights,
        corr::Heteroscedastic
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)
    touse       = nonmissing₁ .* nonmissing₂
    touse₁      = touse[nonmissing₁]
    touse₂      = touse[nonmissing₂]

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, touse₁, i₁) ; lmul!(Diagonal(w₁), ψ₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, touse₂, i₂) ; lmul!(Diagonal(w₂), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, obj₂, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_2s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::UnitWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::UnitWeights,
        corr::Clustered
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(corr.mat, :, nonmissing₁) * view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(corr.mat, :, nonmissing₂) * view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, obj₂, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_2s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::AbstractWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::AbstractWeights,
        corr::Clustered
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(corr.mat, :, nonmissing₁) * view(Ψ₁, :, i₁) ; lmul!(Diagonal(w₁), ψ₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(corr.mat, :, nonmissing₂) * view(Ψ₂, :, i₂) ; lmul!(Diagonal(w₂), ψ₂)
    V₁₂ = ψ₁' * ψ₂

    adjfactor!(V₁₂, obj₁, obj₂, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_2s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::UnitWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::UnitWeights,
        corr::CorrStructure
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, :, i₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, :, i₂)
    V₁₂ = ψ₁' * view(corr.mat, nonmissing₁, nonmissing₂) * ψ₂

    adjfactor!(V₁₂, obj₁, obj₂, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end

function _hausman_2s!(
        output::ParEstimate,
        obj₁::ParModel,
        i₁::Vector{Int},
        w₁::AbstractWeights,
        obj₂::ParModel,
        i₂::Vector{Int},
        w₂::AbstractWeights,
        corr::CorrStructure
    )

    nonmissing₁ = getnonmissing(obj₁)
    nonmissing₂ = getnonmissing(obj₂)

    Ψ₁  = influence(obj₁, w₁)
    ψ₁  = view(Ψ₁, :, i₁) ; lmul!(Diagonal(w₁), ψ₁)
    Ψ₂  = influence(obj₂, w₂)
    ψ₂  = view(Ψ₂, :, i₂) ; lmul!(Diagonal(w₂), ψ₂)
    V₁₂ = ψ₁' * view(corr.mat, nonmissing₁, nonmissing₂) * ψ₂

    adjfactor!(V₁₂, obj₁, obj₂, corr)

    output.β = view(coef(obj₁), i₁) - view(coef(obj₂), i₂)
    output.V = view(vcov(obj₁), i₁, i₁) + view(vcov(obj₂), i₂, i₂) - (V₁₂' + V₁₂)
end


text:
using CircularArrays

import Base.~

# First Half

const maxlbl₁ = 9
~(x::Int, y::Int) = mod1(x - y, maxlbl₁) # Circular subtraction

function simulate1!(cups, moves)

    function place(cups, pickup, destination)
        t = Vector(cups)
        for p in pickup
            pckpos = findfirst(x -> x == p, t)
            deleteat!(t, pckpos)
        end
        destpos = findfirst(x -> x == destination, t)
        CircularVector([t[1:destpos]; pickup; t[destpos + 1:end]])
    end

    current = cups[1]
    destination = 0
    for m = 1:moves
        currpos = findfirst(x -> x == current, cups)
        pickup = cups[currpos + 1:currpos + 3]
        # Step 1 - The crab picks up the three cups that are immediately
        # clockwise of the current cup. They are removed from the circle;

        destination = current ~ 1
        # Step 2 - The crab selects a destination cup: the cup with a label
        # equal to the current cup's label minus one.

        while destination ∈ pickup
            # If this would select one of the cups that was just picked up,
            # the crab will keep subtracting one until it finds a cup that
            # wasn't just picked up.

            destination = destination ~ 1
        end

        println("-- move $m --")
        println("cups: $cups")
        println("current: $current")
        println("pick up: $pickup")
        println("destination: $destination")
        println()

        cups = place(cups, pickup, destination)
        # Step 3 - The crab places the cups it just picked up so that they are
        # immediately clockwise of the destination cup. They keep the same
        # order as when they were picked up.

        currpos = findfirst(x -> x == current, cups)
        current = cups[currpos + 1]
        # Step 4 - The crab selects a new current cup: the cup which is
        # immediately clockwise of the current cup.
    end
    cups
end

puzzleinput = CircularVector([6, 8, 5, 9, 7, 4, 2, 1, 3])
puzzletestinput = CircularVector([3, 8, 9, 1, 2, 5, 4, 6, 7])
simulate1!(puzzleinput, 101)

# Second Half

const maxlbl₂ = 1000000
~(x::Int, y::Int) = mod1(x - y, maxlbl₂) # Circular subtraction

function simulate2!(labeling, moves)

    cups = collect([2:maxlbl2 + 1]...)
    labeling = [labeling; length(labeling) + 1]
    for i ∈ 1:9
        n = labeling[i]
        nx = labeling[i + 1]
        cups[n] = nx
    end
    cups[maxlbl2] = labeling[1]
    # The cups array store at position n the label of the cup immediately
    # clockwise of cup labeled n, that makes all so simple!

    cur = labeling[1]
    des = 0
    pck = [0, 0, 0]
    for m = 1:moves
        pck[1] = cups[cur]s
        pck[2] = cups[pck[1]]
        pck[3] = cups[pck[2]]
        # Step 1 - The crab picks up the three cups that are immediately
        # clockwise of the current cup. They are removed from the circle;

        des = cur ~ 1
        # Step 2 - The crab selects a destination cup: the cup with a label
        # equal to the current cup's label minus one.

        while des ∈ pck
            # If this would select one of the cups that was just picked up,
            # the crab will keep subtracting one until it finds a cup that
            # wasn't just picked up.

            des = des ~ 1
        end

        nexttodes = cups[des]
        cups[des] = pck[1]
        cups[cur] = cups[pck[3]]
        cups[pck[3]] = nexttodes
        # Step 3 - The crab places the cups it just picked up so that they are
        # immediately clockwise of the destination cup. They keep the same
        # order as when they were picked up.

        cur = cups[cur]
        # Step 4 - The crab selects a new current cup: the cup which is
        # immediately clockwise of the current cup.
    end

    cups[1] * cups[cups[1]]
end

puzzletestinput = [3, 8, 9, 1, 2, 5, 4, 6, 7]
puzzleinput = [6, 8, 5, 9, 7, 4, 2, 1, 3]

@time simulate2!(puzzleinput, 10000000)


text:
using NeuralVerification
using NeuralVerification: Layer, n_nodes

_vec(A::AbstractMatrix) = vec(A)
_vec(A::Number) = [A]
_vec(A::AbstractVector) = A

"""
   @modelpath(name)

Return the absolute path to file `name` relative to the executing script.

### Input

- `model_path` -- folder name (ignored by default)
- `name`       -- filename

### Output

A string.

### Notes

This macro is equivalent to `joinpath(@__DIR__, name)`.
The `@modelpath` macro is used in model scripts to load data files relative to the
location of the model, without having to change the directory of the Julia session.
For instance, suppose that the folder `/home/projects/models` contains the script
`my_model.jl`, and suppose that the data file `my_data.dat` located in the same
directory is required to be loaded by `my_model.jl`. Then,

```julia
# suppose the working directory is /home/julia/ and so we ran the script as
# julia -e "include("../projects/models/my_model.jl")"
# in the model file /home/projects/models/my_model.jl we write:
d = open(@modelpath("", "my_data.dat"))
# do stuff with d
```

In this example, the macro `@modelpath("", "my_data.dat")` evaluates to the string
`/home/projects/models/my_data.dat`. If the script `my_model.jl` only had
`d = open("my_data.dat")`, without `@modelpath`, this command would fail as julia
would have looked for `my_data.dat` in the *working* directory, resulting in an
error that the file `/home/julia/my_data.dat` is not found.
"""
macro modelpath(model_path::String, name::String)
    __source__.file === nothing && return nothing
    _dirname = dirname(String(__source__.file))
    dir = isempty(_dirname) ? pwd() : abspath(_dirname)
    return joinpath(dir, name)
end

# ================================================
# Reading a network in MAT format
# ================================================

"""
    read_nnet_mat(file; key=nothing, act_key="activation_fcns")

Read a neural network stored in a `.mat` file and return the corresponding network
in the format of `NeuralVerification.jl`.

### Input

- `file`    -- string indicating the location of the `.mat` file containing the neural network
- `key`     -- (optional, default: `nothing`) key used to search the dictionary containing the controller;
               by default we search the top-level dictionary; a typical value is `"controller"`
- `act_key` -- (optional, default: `"activation_fcns"`) key used to search the activation
               functions; typical values are `"activation_fcns"` or `"act_fcns"`

### Output

A `Network` struct.

### Notes

The following activation functions are supported:

- RELU: "relu" (`ReLU`)
- Identity: "linear" (`Id`)
"""
function read_nnet_mat(file::String; key=nothing, act_key="activation_fcns")
    isdefined(@__MODULE__, :MAT) || error("package 'MAT' is required")

    vars = matread(file)

    # some models store the controller under a specified key
    if !isnothing(key)
        !haskey(vars, key) && throw(ArgumentError("didn't find key $key, existing keys are $(keys(vars))"))
        dic = vars[key]
    else
        dic = vars
    end

    # get number of layers either from a dictionary entry or from the length of the weights array
    if haskey(dic, "number_of_layers")
        m = Int(dic["number_of_layers"])
    else
        m = length(dic["W"])
    end
    layers = Vector{Layer}(undef, m)
    aF = dic[act_key]

    for n = 1:m

        # weights matrix
        W = dic["W"][n]
        s = size(W)
        if length(s) == 4
            # models sometimes are stored as a multi-dim array
            # with two dimensions which are flat
            if s[3] == 1 && s[4] == 1
                W = reshape(W, s[1], s[2])
            else
                throw(ArgumentError("unexpected dimension of the weights matrix: $s"))
            end
        end

        # bias
        b = dic["b"][n]

        # activation function
        if aF[n] == "relu"
            act = ReLU()
        elseif aF[n] == "linear"
            act = Id()
        else
            error("error, aF = $(aF[n]), nLayer = $n")
        end
        layers[n] = Layer(W, _vec(b), act)
    end

    return Network(layers)
end

# ================================================
# Reading a network in YAML format
# (load the data with YAML.jl)
# ================================================

const ACT_YAML = Dict("Id"=>Id(),
                      "ReLU"=>ReLU(),
                      "Sigmoid"=>Sigmoid(),
                      "Tanh"=>Tanh())

"""
    read_nnet_yaml(data::Dict)

Read a neural network from a file in YAML format (see `YAML.jl`) and convert it

Read a neural network stored in a `.mat` file and return the corresponding network
in the format of `NeuralVerification.jl`.

### Input

- `file` -- string indicating the location of the `.mat` file containing the neural network
- `key`  -- (optional, default: `"controller"`) key used to search the dictionary

### Output

A `Network` struct.

### Notes

The following activation functions are supported: identity, relu, sigmoid and tanh;
see `NeuralNetworkAnalysis.ACT_YAML`.
"""
function read_nnet_yaml(data::Dict)
    NLAYERS = length(data["offsets"])
    layers = []
    for k in 1:NLAYERS
        weights = data["weights"][k]
        W = copy(reduce(hcat, weights)')
        b = data["offsets"][k]
        a = ACT_YAML[data["activations"][k]]
        L = Layer(W, b, a)
        push!(layers, L)
    end
    return Network(layers)
end

# ================================================
# Reading a network in Sherlock format
# ================================================

"""
    read_nnet_sherlock(file::String)

Read a neural network from a file in Sherlock format.

### Input

- `file` -- string indicating the location of the input file

### Output

A `Network`.

### Notes

The Sherlock format is the format used by the tool `Sherlock`.
The format is documented
[here](https://github.com/souradeep-111/sherlock/blob/981fbf8d3bc99f76e0135f2f518e97d2a318cb7c/sherlock-network-format.pdf).
All layers including the output layer implicitly use a ReLU activation function.
"""
function read_nnet_sherlock(file::String)
    layers = nothing
    open(file, "r") do io
        n_inputs = parse(Int, readline(io))  # number of neurons in input layer
        n_outputs = parse(Int, readline(io))  # number of neurons in output layer
        n_hlayers = parse(Int, readline(io))  # number of hidden layers
        layers = Vector{Layer}(undef, n_hlayers + 1)

        # one line for each number of neurons in the hidden layers
        n_neurons = Vector{Int}(undef, n_hlayers)
        @inbounds for i in 1:n_hlayers
            n_neurons[i] = parse(Int, readline(io))
        end

        # one line for each weight and bias in the following order:
        # - all incoming weights to neuron 1 in hidden layer 1
        # - bias term of neuron 1 in hidden layer 1
        # - all incoming weights to neuron 2 in hidden layer 1
        # - bias term of neuron 2 in hidden layer 1
        # - ...
        # - all incoming weights to the last neuron in hidden layer 1
        # - bias term of the last neuron in hidden layer 1
        # continue with hidden layer 2 until the output layer
        @inbounds for layer in 1:(n_hlayers + 1)
            m = layer > n_hlayers ? n_outputs : n_neurons[layer]
            n = layer == 1 ? n_inputs : n_neurons[layer - 1]
            W, b = _read_weights_biases_sherlock(io, m, n)
            # the Sherlock format implicitly uses ReLU activation functions
            layers[layer] = Layer(W, b, ReLU())
        end
    end

    return Network(layers)
end

function _read_weights_biases_sherlock(io, m, n)
    W = Matrix{Float32}(undef, m, n)
    b = Vector{Float32}(undef, m)
    for i in 1:m
        for j in 1:n
            W[i, j] = parse(Float32, readline(io))
        end
        b[i] = parse(Float32, readline(io))
    end
    return W, b
end

"""
    write_nnet_sherlock(nnet::Network, file::String)

Write a neural network to a file in Sherlock format.

### Input

- `nnet` -- neural network
- `file` -- string indicating the location of the output file

### Notes

See [`read_nnet_sherlock`](@ref) for information about the Sherlock format.
"""
function write_nnet_sherlock(nnet::Network, file::String)
    layers = nnet.layers
    n_inputs = size(layers[1].weights, 2)
    n_outputs = n_nodes(layers[end])
    n_hlayers = length(layers) - 1  # includes the output layer
    open(file, "w") do io
        println(io, string(n_inputs))  # number of neurons in input layer
        println(io, string(n_outputs))  # number of neurons in output layer
        println(io, string(n_hlayers))  # number of hidden layers

        # one line for each number of neurons in the hidden layers
        @inbounds for i in 1:n_hlayers
            println(io, string(n_nodes(layers[i])))
        end

        # one line for each weight and bias of the hidden and output layers
        @inbounds for layer in layers
            if !(layer.activation isa ReLU)
                @info "the Sherlock format requires ReLU activations, but " *
                    "received a $(typeof(layer.activation)) activation"
            end
            _write_weights_biases_sherlock(io, layer)
        end
    end
    nothing
end

function _write_weights_biases_sherlock(io, layer)
    W = layer.weights
    b = layer.bias
    m, n = size(W)
    for i in 1:m
        for j in 1:n
            println(io, W[i, j])
        end
        println(io, b[i])
    end
end

# ====================
# toy model generation
# ====================

# create a random affine IVP
function random_affine_ivp(n, m)
    A_rand = rand(n, n) .- 0.5  # random values around zero
    b_rand = rand(m) .- 0.5  # random values around zero
    function random_system!(dx, x, params, t; A=A_rand, b=b_rand)
        for i in 1:n
            dx[i] = zero(x[i])
            for j in 1:n
                dx[i] += A[i, j] * x[j]
            end
            for j in 1:m
                dx[i] += b[j] * x[n + j]
            end
        end
        for i in (n+1):(n+m)
            dx[i] = zero(x[i])
        end
    end

    X0 = rand(BallInf, dim=n)
    U0 = ZeroSet(m)  # needed but ignored

    system = BlackBoxContinuousSystem(random_system!, n + m)
    ivp = InitialValueProblem(system, X0 × U0)

    return ivp
end

# create a family of random control problems with the same IVP
function random_control_problems(n::Integer, m::Integer; ivp=nothing, period=0.1)
    if ivp == nothing
        ivp = random_affine_ivp(n, m)
    end

    vars_idx = Dict(:state_vars=>1:n, :control_vars=>(n+1):(n+m))
    problems = ControlledPlant[]

    # create random constant controller
    W = zeros(m, n)
    b = 0.1 .* (rand(m) .- 0.5)  # small random values around zero
    controller = Network([Layer(W, b, Id())])
    push!(problems, ControlledPlant(ivp, controller, vars_idx, period))

    # create random affine controller
    W = 0.1 .* (rand(m, n) .- 0.5)  # small random values around zero
    b = zeros(m)
    controller = Network([Layer(W, b, Id())])
    push!(problems, ControlledPlant(ivp, controller, vars_idx, period))

    # create random affine controller with larger matrix values
    W = 1.1 * (rand(m, n) .- 0.5)  # small random values around zero
    b = zeros(m)
    controller = Network([Layer(W, b, Id())])
    push!(problems, ControlledPlant(ivp, controller, vars_idx, period))

    # create random affine controller with large bias values
    W = 0.1 .* (rand(m, n) .- 0.5)  # small random values around zero
    b = 10 * ones(m)
    controller = Network([Layer(W, b, Id())])
    push!(problems, ControlledPlant(ivp, controller, vars_idx, period))

    # create random controller with a single layer and ReLU activation function
    W = 1.1 * (rand(m, n) .- 0.5)  # small random values around zero
    b = zeros(m)
    controller = Network([Layer(W, b, ReLU())])
    push!(problems, ControlledPlant(ivp, controller, vars_idx, period))

    # create random controller with five layers
    l = 5
    n_max_neurons = 10
    layers = Layer[]
    k_in = n
    for i in 1:l
        k_out = i == l ? m : rand(1:n_max_neurons)
        W = 1.1 * (rand(k_out, k_in) .- 0.5)  # small random values around zero
        b = 0.1 .* (rand(k_out) .- 0.5)  # small random values around zero
        af = i == l ? Id() : ReLU()  # ReLU activation except for the last layer
        push!(layers, Layer(W, b, af))
        k_in = k_out
    end
    controller = Network(layers)
    push!(problems, ControlledPlant(ivp, controller, vars_idx, period))

    return problems
end

output_dim(controller::Network) = size(controller.layers[end].weights, 1)

# relative size between the set-based output and the (CH of) sampled output
function relative_size(X0, nsamples, controller, solver=Ai2())
    @assert output_dim(controller) == 1 "the dimension of the output of the network needs to be 1, but is $output_dim(controller)"
    o = overapproximate(forward_network(solver, controller, X0), Interval)
    u = forward_network(SampledApprox(nsamples), controller, X0)
    return diameter(o)/diameter(u)
end

# print the result `@timed`
function print_timed(stats::NamedTuple)
    Base.time_print(stats.time * 1e9, stats.bytes, stats.gctime * 1e9,
                    Base.gc_alloc_count(stats.gcstats), 0, true)
end


text:
# BEGIN v0.2

# The implementation of NCC is actually ZNCC (#63)
@deprecate NCC ZNCC
@deprecate ncc zncc

# END v0.2


text:
using Printf
# using Plots
using Statistics
using LinearAlgebra
using BenchmarkTools

using PyPlot
pygui(true);

#machine parameters
p=8 #order of matrix power calculation
m=14 #discretization number
mx=0.01986
my=0.02008
kx=1.60312
ky=1.155697
sx=408866
sy=413445

#test change

#cutting parameters
kt=644*10^6
kr=0.368
d=0.008
ae=0.0004
z=1

#calculating cutting teeth
aepu=ae/d
aepd=2-ae/d
phi_inu=0; phi_outu=acos(1-2*aepu);
phi_ind=acos(1-2(aepd-1)); phi_outd=pi;

function phi(j,t,Omega)
    mod(t*Omega+j*2*pi/z,2*pi)
end

function gu(phi::Float64)
    gu_val::Int8=0;
    if phi < phi_inu
        gu_val=0
    elseif phi >= phi_inu && phi <= phi_outu
        gu_val=1
    elseif phi > phi_outu
        gu_val=0
    else
end
return(gu_val)
end

function gd(phi::Float64)
    gd_val::Int8=0;
    if phi < phi_ind
        gd_val=0
    elseif phi >= phi_ind && phi <= phi_outd
        gd_val=1
    elseif phi > phi_outd
        gd_val=0
    else
end
return (gd_val)
end

#delay magnitude
function tau(Omega)
    2*pi/(Omega*z)
end
#timestep
function deltat(Omega)
    tau(Omega)/m
end

#how to plot
x=[n for n in 0:0.2:2*pi]
y=[gu(n) for n in 0:0.2:2*pi]
println(y)
scatter(x,y)

#tool coefficient matrices
Mm=[mx 0; 0 my];
Km=[kx 0; 0 ky];
Sm=[sx 0; 0 sy];
println(Mm," ",Km ," ",Sm)

function W(t,Omega,w)
    W_f=zeros(Float64,2,2)
    for j = 0:(z-1)
        W_f=W_f+gu(phi(j,t,Omega))*[-sin(2*phi(j,t,Omega))-kr+kr*cos(2*phi(j,t,Omega)) -1-cos(2*phi(j,t,Omega))-kr*sin(2*phi(j,t,Omega)); 1-cos(2*phi(j,t,Omega))-kr*sin(2*phi(j,t,Omega)) sin(2*phi(j,t,Omega))-kr-kr*cos(2*phi(j,t,Omega))]
    end
    return ((w*kt/2)*W_f)
end

println(W(0.001,200.0,100.0))

function A(t,Omega,w)
    vcat(hcat(-inv(Mm)*Km,-inv(Mm)*(Sm-W(t,Omega,w))),hcat(Matrix{Float64}(I, 2, 2),zeros(2,2)))
end

function B(t,Omega,w)
    vcat(hcat(zeros(2,2),-inv(Mm)*W(t,Omega,w)),hcat(zeros(2,2),zeros(2,2)))
end

function expM(M)
n=size(M)[1]
Mval=zeros(n,n)
    for j = 0:p
        Mval=Mval+(M^j)/(factorial(j))
    end
return(Mval)
end

function P(i,Omega,w)
    expM(A(i*deltat(Omega),Omega,w)*deltat(Omega))
end

function R(i,Omega,w)
    0.5*((expM(A(i*deltat(Omega),Omega,w)*deltat(Omega))-Matrix{Float64}(I, 4, 4))*inv(A(i*deltat(Omega),Omega,w))*B(i*deltat(Omega),Omega,w))
end

function Zw(i,Omega,w)
    vcat(hcat(P(i,Omega,w),zeros(4,2*(m+2)-8),R(i,Omega,w)[:,3:4],R(i,Omega,w)[:,3:4]),hcat(zeros(2*(m+2)-4,2),Matrix{Float64}(I, 2*(m+2)-4, 2*(m+2)-4),zeros(2*(m+2)-4,2)))
end

function Zws(Omega,w)
    Z=Zw(0,Omega,w)
    for j = 1:(m-1)
        Z=Zw(j,Omega,w)*Z
    end
    return(Z)
end

eigv=eigvals(Zws(1700,0.001))

eigRe=[real(eigv[n]) for n in 1:length(eigv)]
eigIm=[imag(eigv[n]) for n in 1:length(eigv)]
scatter(eigRe,eigIm)
Re_circle=[cos(n) for n in 0:0.05:2*pi]
Im_circle=[sin(n) for n in 0:0.05:2*pi]
plot(ploteig,Re_circle,Im_circle)
display(eigv)
display(ploteig)


text:
module AtBackslash

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end AtBackslash

export @\

using Base.Meta: isexpr

function is_interpolation(ex)
    isexpr(ex, :$) || return false
    if length(ex.args) != 1 || isexpr(ex.args[1], :tuple)
        error("Only single-argument \$ is supported. Got: ", ex)
    end
    return true
end

"""
    @\\(NCT1, NTC2, ..., NTCn) :: NamedTuple
    @\\(; NCT1, NTC2, ..., NTCn) :: NamedTuple
    @\\(expr)
"""
macro \(args...)
    @gensym input

    replace_symbols(x) = x
    replace_symbols(x::QuoteNode) =
        if x.value isa Symbol
            :($input.$(x.value))
        else
            x
        end
    replace_symbols(x::Symbol) =
        if x === :_
            input
        else
            esc(x)
        end
    replace_symbols(ex::Expr) =
        if isexpr(ex, :...) && ex.args == [:_]
            :($input...)
        elseif isexpr(ex, :parameters)
            args = map(ex.args) do x
                if x isa Symbol
                    Expr(:(=), esc(x), replace_symbols(x))
                elseif x isa QuoteNode && x.value isa Symbol
                    Expr(:(=), esc(x.value), replace_symbols(x))
                else
                    replace_symbols(x)
                end
            end
            Expr(ex.head, args...)
        elseif is_interpolation(ex)
            # Should I still replace `_`?
            esc(ex.args[1])
        elseif isexpr(ex, :.)
            Expr(:., replace_symbols(ex.args[1]), ex.args[2:end]...)
        elseif isexpr(ex, :macrocall)
            esc(ex)
        elseif isexpr(ex, :quote)
            ex
        else
            Expr(ex.head, replace_symbols.(ex.args)...)
        end

    if args == (:_,)
        # Should I?
        error("`@\\_` not supported. Use `identity`.")
    end

    #! format: off
    if (
        length(args) == 1 &&
        ((args[1] isa QuoteNode && args[1].value isa Symbol) || args[1] isa Symbol)
    )
        ex, = args
        error(
            "Ambiguous expression: `@\\($ex)`\n",
            "Use `@\\(_.$ex)` if you mean `x -> x.$ex`.\n",
            "Use `@\\(; $ex)` if you mean `x -> ($ex = x.$ex,)`.",
        )
    end
    #! format: on

    # Handle: @\(_..., a=x+1, b) :: NamedTuple
    all(args) do x
        (x isa Symbol ||
         (x isa QuoteNode && x.value isa Symbol) ||
         (is_interpolation(x) && x.args[1] isa Symbol) ||
         (isexpr(x, :(=)) && x.args[1] isa Symbol) ||
         (isexpr(x, :...) && x.args == [:_]))
    end && begin
        return :($input -> $(replace_symbols(Expr(:tuple, Expr(:parameters, args...)))))
    end

    # Handle: @\(; KEY => VAL) :: NamedTuple
    if length(args) == 1 && isexpr(args[1], :parameters)
        return :($input -> $(replace_symbols(Expr(:tuple, args[1]))))
    end

    # Handle: @\[a, b] :: Array
    # Handle: @\(x > y) :: Bool
    if length(args) == 1
        return :($input -> $(replace_symbols(args[1])))
    end

    error("Unsupported expression:\n", "@\\", Expr(:tuple, args...))
end

end # module


text:
"""
    ∇²f!(output, f, γ)

Finite difference operator for a second derivative in two dimensions.

Computes the laplacian of an input `f` times a scalar `γ` and stores the result in `output`.

# Mathematics

The laplacian operator in two dimensions can be written as

`` \\nabla^2 f = \\frac{\\partial^2 f}{\\partial x^2} + \\frac{\\partial^2 f}{\\partial y}. ``

For the discretization of this operator we use a nine point stencil, such the neighbors as well as the diagonal elements.
The concrete derivation can be found in the references below, we just show the final result

``\\nabla^2 f = \\frac{1}{6}\\bigg[4(f_{i+1,j} + f_{i,j+1} + f_{i-1,j} + f_{i,j-1}) \\newline
                \\qquad\\qquad +(f_{i+1,j+1} + f_{i-1,j+1} + f_{i-1,j-1} + f_{i+1,j-1}) \\newline
                \\qquad\\qquad -20f_{i,j}\\bigg]  ,``

where we have used Julia conventions, downwards (left) is positive.
The whole expression can be multiplied with a scalar `γ` if needed.

# Examples

```jldoctest
julia> using Swalbe, Test

julia> arg = reshape(collect(1.0:25),5,5)
5×5 Matrix{Float64}:
 1.0   6.0  11.0  16.0  21.0
 2.0   7.0  12.0  17.0  22.0
 3.0   8.0  13.0  18.0  23.0
 4.0   9.0  14.0  19.0  24.0
 5.0  10.0  15.0  20.0  25.0

julia> res = zeros(5,5); Swalbe.∇²f!(res, arg, -1.0)

julia> analytics = [-30.0 -5.0 -5.0 -5.0 20;
                    -25.0 0.0 0.0 0.0 25.0;
                    -25.0 0.0 0.0 0.0 25.0;
                    -25.0 0.0 0.0 0.0 25.0;
                    -20.0 5.0 5.0 5.0 30.0];

julia> for i in eachindex(analytics)
           @test analytics[i] ≈ res[i] atol=1e-10
       end
```

# References

- [Junk & Klar](https://epubs.siam.org/doi/10.1137/S1064827599357188)
- [Succi et al.](https://doi.org/10.1016/j.jcp.2012.07.037)

See also: [`Swalbe.∇f!`](@ref)
"""
function ∇²f!(output, f, γ)
    # Straight elements j+1, i+1, i-1, j-1
    hip = circshift(f, (1,0))
    hjp = circshift(f, (0,1))
    him = circshift(f, (-1,0))
    hjm = circshift(f, (0,-1))
    # Diagonal elements
    hipjp = circshift(f, (1,1))
    himjp = circshift(f, (-1,1))
    himjm = circshift(f, (-1,-1))
    hipjm = circshift(f, (1,-1))
    # In the end it is just a weighted sum...
    output .= γ .* (2/3 .* (hjp .+ hip .+ him .+ hjm)
                 .+ 1/6 .* (hipjp .+ himjp .+ himjm .+ hipjm)
                 .- 10/3 .* f)
    return nothing
end

function ∇²f!(output, f::Vector, dgrad)
    hip, him = viewneighbors_1D(dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(hip, f, 1)
    circshift!(him, f, -1)
    # In the end it is just a weighted sum...
    output .= hip .- 2 .* f .+ him
    return nothing
end

"""
    ∇f!(outputx, outputy, f)

Gradient calculation with finite differences.

Computes both spatial first derivatives with a nine point stencil from an input `f` and writes the result to `outputx` and `outputy`.
Since broadcasting is simple on the **GPU** we make use of `circshift` for the neighbors.

# Mathematics

The gardient in two dimensions is given as

`` \\nabla f = \\big(\\frac{\\partial f}{\\partial x}, \\frac{\\partial f}{\\partial y}\\big)^T .``

Again with the nine point stencil this reduces to

`` \\frac{\\partial f}{\\partial x} = \\frac{1}{3} (f_{i+1,j} - f_{i-1,j}) + \\frac{1}{12}(f_{i+1,j+1} - f_{i-1,j+1} - f_{i-1,j-1} + f_{i+1,j-1}) ,``

and for the `y` component we get

`` \\frac{\\partial f}{\\partial y} = \\frac{1}{3} (f_{i,j+1} - f_{i,j-1}) + \\frac{1}{12}(f_{i+1,j+1} + f_{i-1,j+1} - f_{i-1,j-1} - f_{i+1,j-1}) .``

For the exact derivation feel free to read the reference by Junk and Klar.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> arg = reshape(collect(1.0:25),5,5)
5×5 Matrix{Float64}:
 1.0   6.0  11.0  16.0  21.0
 2.0   7.0  12.0  17.0  22.0
 3.0   8.0  13.0  18.0  23.0
 4.0   9.0  14.0  19.0  24.0
 5.0  10.0  15.0  20.0  25.0

julia> resx = zeros(5,5); resy = zeros(5,5); Swalbe.∇f!(resx, resy, arg)

julia> whatXshouldbe = [-1.5 -1.5 -1.5 -1.5 -1.5;
                         1.0 1.0 1.0 1.0 1.0;
                         1.0 1.0 1.0 1.0 1.0;
                         1.0 1.0 1.0 1.0 1.0;
                        -1.5 -1.5 -1.5 -1.5 -1.5];

julia> for i in eachindex(resx) # Test the x-component
           @test resx[i] ≈ whatXshouldbe[i] atol=1e-10
       end

julia> whatYshouldbe = [-7.5 5.0 5.0 5.0 -7.5;
                        -7.5 5.0 5.0 5.0 -7.5;
                        -7.5 5.0 5.0 5.0 -7.5;
                        -7.5 5.0 5.0 5.0 -7.5;
                        -7.5 5.0 5.0 5.0 -7.5];

julia> for i in eachindex(resy) # Test the y-component
           @test resy[i] ≈ whatYshouldbe[i] atol=1e-10
       end
```

# References

- [Junk & Klar](https://epubs.siam.org/doi/10.1137/S1064827599357188)
- [Succi et al.](https://doi.org/10.1016/j.jcp.2012.07.037)

See also: [`Swalbe.∇²f!`](@ref)
"""
function ∇f!(outputx, outputy, f)
    # Straight elements j+1, i+1, i-1, j-1
    fip = circshift(f, (1,0))
    fjp = circshift(f, (0,1))
    fim = circshift(f, (-1,0))
    fjm = circshift(f, (0,-1))
    # Diagonal elements
    fipjp = circshift(f, (1,1))
    fimjp = circshift(f, (-1,1))
    fimjm = circshift(f, (-1,-1))
    fipjm = circshift(f, (1,-1))
    # In the end it is just a weighted sum...
    outputx .= -1/3 .* (fip .- fim) .- 1/12 .* (fipjp .- fimjp .- fimjm .+ fipjm)
    outputy .= -1/3 .* (fjp .- fjm) .- 1/12 .* (fipjp .+ fimjp .- fimjm .- fipjm)

    return nothing
end

function ∇f!(outputx, outputy, f, a)
    # Straight elements j+1, i+1, i-1, j-1
    fip = circshift(f, (1,0))
    fjp = circshift(f, (0,1))
    fim = circshift(f, (-1,0))
    fjm = circshift(f, (0,-1))
    # Diagonal elements
    fipjp = circshift(f, (1,1))
    fimjp = circshift(f, (-1,1))
    fimjm = circshift(f, (-1,-1))
    fipjm = circshift(f, (1,-1))
    # In the end it is just a weighted sum...
    outputx .= a .* (-1/3 .* (fip .- fim) .- 1/12 .* (fipjp .- fimjp .- fimjm .+ fipjm))
    outputy .= a .* (-1/3 .* (fjp .- fjm) .- 1/12 .* (fipjp .+ fimjp .- fimjm .- fipjm))

    return nothing
end

function ∇f!(outputx, outputy, f, dgrad, a)
    fip, fjp, fim, fjm, fipjp, fimjp, fimjm, fipjm = viewneighbors(dgrad)
    # Straight elements j+1, i+1, i-1, j-1
    circshift!(fip, f, (1,0))
    circshift!(fjp, f, (0,1))
    circshift!(fim, f, (-1,0))
    circshift!(fjm, f, (0,-1))
    # Diagonal elements
    circshift!(fipjp, f, (1,1))
    circshift!(fimjp, f, (-1,1))
    circshift!(fimjm, f, (-1,-1))
    circshift!(fipjm, f, (1,-1))
    # In the end it is just a weighted sum...
    outputx .= a .* (-1/3 .* (fip .- fim) .- 1/12 .* (fipjp .- fimjp .- fimjm .+ fipjm))
    outputy .= a .* (-1/3 .* (fjp .- fjm) .- 1/12 .* (fipjp .+ fimjp .- fimjm .- fipjm))

    return nothing
end

function ∇f!(output::Vector, f, dgrad, a)
    fip, fim = viewneighbors_1D(dgrad)
    # One dim case, central differences
    circshift!(fip, f, 1)
    circshift!(fim, f, -1)

    # In the end it is just a weighted sum...
    output .= a .* -0.5 .* (fip .- fim)

    return nothing
end

function ∇f!(output, f::Vector, dgrad)
    fip, fim = viewneighbors_1D(dgrad)
    # One dim case, central differences
    circshift!(fip, f, 1)
    circshift!(fim, f, -1)

    # In the end it is just a weighted sum...
    output .= -0.5 .* (fip .- fim)

    return nothing
end

"""
    viewneighbors(f)

Generates a view for all nine populations of a **D2Q9** distribution function.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> ftest = reshape(collect(1.0:5*5*8),5,5,8);

julia> f1, f2, f3, f4, f5, f6, f7, f8 = Swalbe.viewneighbors(ftest);

julia> @test all(f3 .== ftest[:,:,3])
Test Passed

```

See also: [`Swalbe.∇f!`](@ref), [`Swalbe.filmpressure!`](@ref)
"""
function viewneighbors(f)
    f1 = view(f, :, :, 1)
    f2 = view(f, :, :, 2)
    f3 = view(f, :, :, 3)
    f4 = view(f, :, :, 4)
    f5 = view(f, :, :, 5)
    f6 = view(f, :, :, 6)
    f7 = view(f, :, :, 7)
    f8 = view(f, :, :, 8)

    return f1, f2, f3, f4, f5, f6, f7, f8
end

"""
    viewneighbors_1D(f)

Generates a view for the two neighbors of a **D1Q3** distribution function.

# Examples
```jldoctest
julia> using Swalbe, Test

julia> ftest = reshape(collect(1.0:30),10,3);

julia> f1, f2 = Swalbe.viewneighbors_1D(ftest);

julia> @test all(f2 .== ftest[:,2])
Test Passed

```

See also: [`Swalbe.∇f!`](@ref), [`Swalbe.filmpressure!`](@ref)
"""
function viewneighbors_1D(f)
    f1 = view(f, :, 1)
    f2 = view(f, :, 2)

    return f1, f2
end


text:
# Monte Carlo solution for arbitrary systems.

"""
Sampling parameters for a particular system.
"""
struct SamplingParameters{S,M,P}
    sys::DiagonalSystem{S,M}

    "Imaginary time step."
    tau::Float64

    "Surface weights."
    weights::Weights
    "Precision (inverse covariance) matrices (M, S)."
    precs::Matrix{Matrix{Float64}}
    "Multivariate normal distributions (M, S)."
    mvns::Matrix{MvNormal}

    "Precomputed cosh factors (M, S)."
    Cs::Matrix{Float64}
    "Precomputed sinh factors (M, S)."
    Ss::Matrix{Float64}
    "Precomputed products of sinh factors (S)."
    S_prods::Vector{Float64}
end

"""
    SamplingParameters(sys::DiagonalSystem, beta::Float64, P::Int)

Generate sampling parameters for `sys` at `beta` with `P` links.
"""
function SamplingParameters(sys::DiagonalSystem{S,M}, beta::Float64, P::Int) where {S,M}
    issimple(sys) || throw(DomainError(:sys, "System for sampling must be simple."))

    tau = beta / P

    Zas = ones(S)
    for s in 1:S
        for m in 1:M
            Zas[s] *= 1.0 / (2*sinh(0.5 * beta * sys.freq[m, s]))
        end
    end
    preweights = exp.(-beta * (diag(sys.energy) + sys.deltas))
    weights = Weights(preweights .* Zas)

    Cs = cosh.(sys.freq * tau)
    Ss = sinh.(sys.freq * tau)
    S_prods = Float64[prod(Ss[:, s]) for s in 1:S]

    precs = Matrix{Matrix{Float64}}(undef, M, S)
    mvns = Matrix{MvNormal}(undef, M, S)
    for s in 1:S
        for m in 1:M
            mean = sys.ds[m, s] * ones(P)

            # Precision matrix.
            prec = diagm(0 => 2 * Cs[m, s] * ones(P), -1 => -ones(P-1), 1 => -ones(P-1))
            prec[1, end] = -1.0
            prec[end, 1] = -1.0
            prec ./= Ss[m, s]

            # Convariance matrix.
            cov = inv(prec)
            maximum(abs.(cov' - cov)) < 1e-13 || @warn "Asymmetric cov: $(maximum(abs.(cov' - cov)))"
            # Force it to be symmetric.
            cov = (cov + cov') / 2

            precs[m, s] = prec
            mvns[m, s] = MvNormal(mean, cov)
        end
    end

    SamplingParameters{S,M,P}(sys, tau, weights, precs, mvns, Cs, Ss, S_prods)
end

function Base.show(io::IO, sp::SamplingParameters{S,M,P}) where {S,M,P}
    println(io, typeof(sp))
    println(io, "Sampling parameters with $(P) link$(P == 1 ? "" : "s") for a system with $(S) surface$(S == 1 ? "" : "s") and $(M) mode$(M == 1 ? "" : "s").")
    println(io, "imaginary time step (tau): $(sp.tau)")
    println(io, "sampling weights:")
    show_vector(io, weights(sp))
    nothing
end

"""
    weights(sp::SamplingParameters)

Compute the normalized surface weights for `sp`.
"""
weights(sp::SamplingParameters) = sp.weights/sum(sp.weights)

"""
    weights(sys::DiagonalSystem, beta::Float64)

Compute the normalized surface weights for `sys` at `beta`.
"""
weights(sys::DiagonalSystem, beta::Float64) = weights(SamplingParameters(sys, beta, 2))

"""
Monte Carlo solution for an arbitrary system.
"""
abstract type Sampling <: Solution end

function sampling_matrix_free_particle(sp::SamplingParameters{S,M,P}, qs1::Vector{Float64}, qs2::Vector{Float64}; scaling=nothing) where {S,M,P}
    FM = zeros(S, S)
    for s in 1:S
        x = -sp.tau * (sp.sys.energy[s, s] + sp.sys.deltas[s])
        for m in 1:M
            q_disp1 = qs1[m] - sp.sys.ds[m, s]
            q_disp2 = qs2[m] - sp.sys.ds[m, s]
            x += -0.5 * ((q_disp1^2 + q_disp2^2) * sp.Cs[m, s] - 2 * q_disp1 * q_disp2) / sp.Ss[m, s]
        end
        FM[s, s] += exp(x) / sqrt(sp.S_prods[s])
    end
    scaling === nothing && (scaling = maximum(FM))
    FM/scaling, scaling
end

function sampling_matrix_interaction(sys::System{S,M}, sp::SamplingParameters{S,M,P}, qs::Vector{Float64}) where {S,M,P}
    preIM = zeros(S, S)
    for s1 in 1:S
        # Only build the upper triangle.
        for s2 in 1:s1
            if s1 != s2
                preIM[s2, s1] += sys.energy[s2, s1]
                for m in 1:M
                    preIM[s2, s1] += sys.lin[m, s2, s1] * qs[m]
                end
            end
            for m1 in 1:M
                for m2 in 1:M
                    preIM[s2, s1] += 0.5 * sys.quad[m2, m1, s2, s1] * qs[m2] * qs[m1]
                end
            end
        end
    end
    Symmetric(preIM), exp(Symmetric(-sp.tau * preIM))
end

function sampling_matrix_energy(sp::SamplingParameters{S,M,P}, qs1::Vector{Float64}, qs2::Vector{Float64}) where {S,M,P}
    EM2 = zeros(S, S)
    for s in 1:S
        EM2[s, s] -= sp.sys.energy[s, s] + sp.sys.deltas[s]
        for m in 1:M
            q_disp1 = qs1[m] - sp.sys.ds[m, s]
            q_disp2 = qs2[m] - sp.sys.ds[m, s]
            EM2[s, s] -= sp.sys.freq[m, s] * 0.5 * sp.Cs[m, s] / sp.Ss[m, s]
            EM2[s, s] -= sp.sys.freq[m, s] * q_disp1 * q_disp2 * sp.Cs[m, s] / sp.Ss[m, s]^2
            EM2[s, s] += sp.sys.freq[m, s] * 0.5 * (q_disp1^2 + q_disp2^2) / sp.Ss[m, s]^2
        end
    end
    EM2
end

function sampling_matrix_heat_capacity(sp::SamplingParameters{S,M,P}, qs1::Vector{Float64}, qs2::Vector{Float64}) where {S,M,P}
    CM1 = zeros(S, S)
    for s in 1:S
        for m in 1:M
            q_disp1 = qs1[m] - sp.sys.ds[m, s]
            q_disp2 = qs2[m] - sp.sys.ds[m, s]
            CM1[s, s] += sp.sys.freq[m, s]^2 * 0.5 / sp.Ss[m, s]^2
            CM1[s, s] += sp.sys.freq[m, s]^2 * 0.5 * q_disp1 * q_disp2 * (sp.Ss[m, s]^2 + sp.Cs[m, s]^2 + 3) / sp.Ss[m, s]^3
            CM1[s, s] -= sp.sys.freq[m, s]^2 * (q_disp1^2 + q_disp2^2) * sp.Cs[m, s] / sp.Ss[m, s]^3
        end
    end
    CM1
end


text:
pts = invariantize([rand(3) for i = 1:15])
triang = triangulate(pts)

μ_exact = invariantmeasure(pts, TriangulationBinning(), ExactIntersection())
μ_approx = invariantmeasure(pts, TriangulationBinning(), ApproximateIntersection())

@test μ_exact isa TransferOperatorTriangulationExact
@test μ_approx isa TransferOperatorTriangulationApprox


text:
using AbstractAlgebra # implements arbitrary precision rationals

tanplus(x,y) = (x + y) / (1 - x * y)

function taneval(coef, frac)
    if coef == 0
        return 0
    elseif coef < 0
        return -taneval(-coef, frac)
    elseif isodd(coef)
        return tanplus(frac, taneval(coef - 1, frac))
    else
        x = taneval(div(coef, 2), frac)
        return tanplus(x, x)
    end
end

taneval(tup::Tuple) = taneval(tup[1], tup[2])

tans(v::Vector{Tuple{BigInt, Rational{BigInt}}}) = foldl(tanplus, map(taneval, v), init=0)

const testmats = Dict{Vector{Tuple{BigInt, Rational{BigInt}}}, Bool}([
    ([(1, 1//2), (1, 1//3)], true), ([(2, 1//3), (1, 1//7)], true),
    ([(12, 1//18), (8, 1//57), (-5, 1//239)], true),
    ([(88, 1//172), (51, 1//239), (32, 1//682), (44, 1//5357), (68, 1//12943)], true),
    ([(88, 1//172), (51, 1//239), (32, 1//682), (44, 1//5357), (68, 1//12944)], false)])

function runtestmats()
    println("Testing matrices:")
    for (k, m) in testmats
        ans = tans(k)
        println((ans == 1) == m ? "Verified as $m: " : "Not Verified as $m: ", "tan $k = $ans")
    end
end

runtestmats()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Cornerback

include("ReceiverLogic.jl")
include("CornerbackLogic.jl")

export fade, post

end


text:
module MCMC

using LinearAlgebra
using Distributions

export update_mean, update_var, metropolis, adaptive_metropolis

eye(n) = diagm(0 => ones(Float64, n))

## posterior interval
pint(theta::Array{Float64,2}) = [quantile(theta[:,j], [0.025, 0.975]) for j in 1:size(theta,2)]

update_mean(t::Int, Xt::Float64, m0::Float64) = (Xt + m0*(t-1)) / t
update_mean(t::Int, Xt::Array{Float64,1}, m0::Array{Float64,1}) = @. (Xt + m0*(t-1)) / t

function update_var(t::Int, Xt::Float64, mt::Float64, m0::Float64, V0::Float64)
    return V0*(t-2)/(t-1) + t*(mt-m0)^2
end

function update_var(t::Int, Xt::Array{Float64,1}, mt::Array{Float64,1},
                   m0::Array{Float64,1}, V0::Array{Float64,2}, eps::Matrix{Float64})
    return V0.*(t-2)/(t-1) + 1/(t-1).*((t-1).*(m0 * transpose(m0)) .- t.* (mt * transpose(mt)) .+ (Xt * transpose(Xt)) + eps)
end

## single Metropolis-Hastings step
function metropolis(target_log::Function, init, proposal::Function, proposal_lpdf::Function)
    prop = proposal(init)
    q1 = proposal_lpdf(prop, init)
    q2 = proposal_lpdf(init, prop)

    r = q1 + target_log(prop) - (q2 + target_log(init))

    if log(rand(1)[1]) < r
        return prop, 1
    else
        return init, 0
    end
end

function metropolis(target_log::Function, init::T, prop::T) where T <: Real
    r = target_log(prop) - target_log(init)

    if log(rand()) < r
        return prop, 1
    else
        return init, 0
    end
end

metropolis(target_log::Function, init, proposal::Function) = metropolis(target_log, init, proposal(init))

function metropolis(target_log::Function, init::T, proposal::Function, init_loglik::T) where T <: Real
    prop = proposal(init)
    tlp = target_log(prop)
    tli = target_log(init, init_loglik)

    if log(rand()) < (tlp[1] - tli[1])
        return prop, 1, tlp[2]
    else
        return init, 0, tlp[2]
    end
end

## single adaptive Metropolis step
function adaptive_metropolis(target_log::Function, init::Vector{Float64}, V::Matrix{Float64},
    t::Int, t0::Int=2*length(init), s::Float64=(2.38)^2/length(init), beta::Float64=0.05)

    d = length(init)
    proposal = zeros(Float64, d)

    if (t <= t0) | (rand() < beta)
        proposal = rand(MvNormal(init, (0.1)^2*eye(d)/d))
    else
        proposal = rand(MvNormal(init, s.*V))
    end

    r = target_log(proposal) - target_log(init)

    if log(rand()) < r
        return proposal, 1
    else
        return init, 0
    end
end

function adaptive_metropolis(target_log::Function, init::Float64, V::Float64, t::Int, t0::Int=2, s::Float64=2.38, beta::Float64=0.05)

    proposal::Float64 = 0

    if (t <= t0) | (rand() < beta)
        proposal = rand(Normal(init, 0.1))
    else
        if(s*V == 0)
            proposal = init
        else
            proposal = rand(Normal(init, s*sqrt(V)))
        end
    end

    r = target_log(proposal) - target_log(init)

    if log(rand()) < r
        return proposal, 1
    else
        return init, 0
    end
end

end


text:
function   js()
    println("Well done ma boi")
end

function last_update()
    println("Last update: March 03 2021")
end


text:
using Pickle
using BSON

@doc """
    collect_samples(path::String)

Collect samples stored in pickle format in given 'path'.
""" ->
function collect_samples(path::String)
    names::Array{String} = readdir(path)
    dataset::Array{Float64, 3} = zeros(4, 256, 1)
    for name::String in names
        d::Matrix{Float64} = hcat(Pickle.load(open(join([path, name])))...)
        dataset = cat(dataset, d[:, 1:256, :]; dims=3)
    end
    return dataset[:,:,2:end];
end

@doc """
    save_dataset(mal_path::String, clean_path::String, path::String)

Load malware ('mal_path') and cleanware ('clean_path') samples. Collect them and save
complete dataset in BSON format in given 'path'.
""" ->
function save_dataset(mal_path::String, clean_path::String, path::String)
    malwares = collect_samples(mal_path)
    cleanwares = collect_samples(clean_path)
    bson(path, m = malwares, c = cleanwares)
end

@doc """
    load_dataset(path::String; dataset_split::Float64=0.8)

Load BSON dataset from given 'path' and split it to train and validation parts
(default split coefficient 'dataset_split' is 0.8). Returns two tuples.
""" ->
function load_dataset(path::String; dataset_split::Float64=0.8)
    dataset = BSON.load(path)
    dataset[:m] = permutedims(dataset[:m], (2,1,3))
    dataset[:m] = dataset[:m][:,:,shuffle(1:end)]
    dataset[:c] = permutedims(dataset[:c], (2,1,3))
    dataset[:c] = dataset[:c][:,:,shuffle(1:end)]

    mal_cnt::Int64 = size(dataset[:m])[3]
    clean_cnt::Int64 = size(dataset[:c])[3]
    trn_mal = dataset[:m][:,:,1:Int64(floor(dataset_split*mal_cnt))]
    val_mal = dataset[:m][:,:,Int64(floor(dataset_split*mal_cnt)):end]
    trn_clean = dataset[:c][:,:,1:Int64(floor(dataset_split*clean_cnt))]
    val_clean = dataset[:c][:,:,Int64(floor(dataset_split*clean_cnt)):end]

    trn_data = Array{Float32, 3}(cat(trn_mal, trn_clean; dims=3))
    val_data = Array{Float32, 3}(cat(val_mal, val_clean; dims=3))
    trn_labels = Vector{Int64}(cat(ones(size(trn_mal)[3]),
                        zeros(size(trn_clean)[3]); dims=1))
    val_labels = Vector{Int64}(cat(ones(size(val_mal)[3]),
                        zeros(size(val_clean)[3]); dims=1))

    trn_dataset = (trn_data, trn_labels)
    val_dataset = (val_data, val_labels)

    return trn_dataset, val_dataset
end

max_stars_repo_name: kozvojtex/MMWavenet

max_issues_repo_name: kozvojtex/MMWavenet

max_forks_repo_name: kozvojtex/MMWavenet


text:
type BrentSolver <: Solver1D
  solverInfo::SolverInfo
end
BrentSolver(maxEvals::Int = 100, lowerBoundEnforced::Bool = false, upperBoundEnforced::Bool = false, lowerBound::Float64 = 0.0, upperBound::Float64 = 0.0) =
  BrentSolver(SolverInfo(maxEvals, lowerBoundEnforced, upperBoundEnforced, lowerBound, upperBound))

# brent solver function
function _solve(solver::BrentSolver, f::Function, accuracy::Float64, xMin::Float64, xMax::Float64, fxMin::Float64, fxMax::Float64, root::Float64, eval_num::Int)
  froot = f(root)
  eval_num += 1
  max_evals = solver.solverInfo.maxEvals
  if froot * fxMin < 0
    xMax = xMin
    fxMax = fxMin
  else
    xMin = xMax
    fxMin = fxMax
  end

  d = root - xMax
  e = d

  while eval_num < max_evals
    if (froot > 0.0 && fxMax > 0.0) || (froot < 0.0 && fxMax < 0.0)
      # rename xMin, root, xMax, and adjust the bounds
      xMax = xMin
      fxMax = fxMin
      e = d = root - xMin
    end

    if abs(fxMax) < abs(froot)
      xMin = root
      root = xMax
      xMax = xMin
      fxMin = froot
      froot = fxMax
      fxMax = fxMin
    end

    # convergence check
    xAcc1 = 2.0 * eps(Float64) * abs(root) + 0.5 * accuracy
    xMid = (xMax - root) / 2.0
    if abs(xMid) <= xAcc1 || is_close(froot, 0.0)
      _ = f(root)
      eval_num += 1
      return root
    end

    if abs(e) >= xAcc1 && abs(fxMin) > abs(froot)
      # attempting inverse quadratic interpolation
      s = froot / fxMin
      if is_close(xMin, xMax)
        p = 2.0 * xMid * s
        q = 1.0 - s
      else
        q = fxMin / fxMax
        r = froot / fxMax
        p = s * (2.0 * xMid * q * (q - r) - (root - xMin) * (r - 1.0))
        q = (q - 1.0) * (r - 1.0) * (s - 1.0)
      end

      if p > 0.0
        q = -q # check if in bounds
      end

      p = abs(p)
      min1 = 3.0 * xMid * q - abs(xAcc1 * q)
      min2 = abs(e * q)

      if (2.0 * p < (min1 < min2 ? min1 : min2))
        e = d
        d = p / q
      else
        # interpolation failed, used bisection
        d = xMid
        e = d
      end
    else
      # bounds are increasing too slowly
      d = xMid
      e = d
    end
    xMin = root
    fxMin = froot
    if abs(d) > xAcc1
      root += d
    else
      root += sign(xMid) * xAcc1
    end

    froot = f(root)
    eval_num += 1
  end

  error("Maximum number of function evals exceeded!")
end


text:
using StructC14N
using Test

template = (xrange=NTuple{2,Number},
            yrange=NTuple{2,Number},
            title="A string")

c = canonicalize(template,
              (xr=(1,2), tit="Foo"))

@test c.xrange == (1, 2)
@test ismissing(c.yrange)
@test c.title == "Foo"

c = canonicalize(template,
              ((1,2), (3.3, 4.4), "Foo"))

@test c.xrange == (1, 2)
@test c.yrange == (3.3, 4.4)
@test c.title == "Foo"

c = canonicalize(template,
              xr=(21,22), tit="Bar")
@test c.xrange == (21, 22)
@test ismissing(c.yrange)
@test c.title == "Bar"

mutable struct AStruct
    xrange::NTuple{2, Number}
    yrange::NTuple{2, Number}
    title::String
end

template = AStruct((0,0), (0., 0.), "A string")

c = canonicalize(template,
              (xr=(1,2), tit="Foo"))

@test c.xrange == (1, 2)
@test c.yrange == (0., 0.)
@test c.title == "Foo"

c = canonicalize(template,
              ((1,2), (3.3, 4.4), "Foo"))

@test c.xrange == (1, 2)
@test c.yrange == (3.3, 4.4)
@test c.title == "Foo"

c = canonicalize(template,
              ((11,12), (13.3, 14.4), "Foo"))
@test c.xrange == (11, 12)
@test c.yrange == (13.3, 14.4)
@test c.title == "Foo"

c = canonicalize(template,
              xr=(21,22), tit="Bar")
@test c.xrange == (21, 22)
@test c.yrange == (0., 0.)
@test c.title == "Bar"

function wrapper(template; kwargs...)
    return canonicalize(template; kwargs...)
end
c = wrapper(template; xr=(31,32), tit="BAZ")
@test c.xrange == (31, 32)
@test c.yrange == (0., 0.)
@test c.title == "BAZ"

configtemplate = (optStr=String,
                  optInt=Int,
                  optFloat=Float64)

configentry = "aa, 1, 2"
c = canonicalize(configtemplate, (split(configentry, ",")...,))
@test c.optStr == "aa"
@test c.optInt == 1
@test c.optFloat == 2.0

configentry = "optFloat=20, optStr=\"aaa\", optInt=10"
c = canonicalize(configtemplate, eval(Meta.parse("($configentry)")))
@test c.optStr == "aaa"
@test c.optInt == 10
@test c.optFloat == 20.0


text:
# Module containing code for Interplanetary Navigation and Guidance class
# Christopher Binz
# University of Maryland

module Interplanetary

using Sundials,
    Roots,
    LinearAlgebra

export muSun, muEarth, AU, RE

const muSun   = 1.3271244004193938e11 # km^3/s^2
const muEarth = 3.98600440e5          # km^3/s^2
const AU      = 1.495978707e8         # km
const RE      = 6378.137              # km

include("utilities.jl")
include("frames.jl")
include("dynamics.jl")
include("transfers.jl")

end # end module


text:
using Koma, JLD2, MAT
fantoma = matread("/home/ccp/Downloads/fantoma.mat")["fantoma"]

ρ = fantoma["PD"]
Δx = 2e-3
M, N = size(ρ)
FOVx = (M-1)*Δx #[m]
FOVy = (N-1)*Δx #[m]
xx = -FOVx/2:Δx:FOVx/2 #spin coordinates
yy = -FOVy/2:Δx:FOVy/2 #spin coordinates
x = xx   .+ yy'*0
y = xx*0 .+ yy' #grid points
z = y*0.
v = 1 # m/s
phantom = Phantom(name="brain",
                  x= y[ρ.!=0],
                  y=-x[ρ.!=0],
                  z= z[ρ.!=0],
                  ρ=fantoma["PD"][ρ.!=0],
                  T1=fantoma["T1"][ρ.!=0]*1e-3,
                  T2=fantoma["T2"][ρ.!=0]*1e-3,
                  Δw=2π*fantoma["df"][ρ.!=0]*2,
                  ux=(x,y,z,t)->v*t
                  )
@save "/home/ccp/Downloads/brain.phantom" phantom #Hacer ]up para version Koma v0.3.6