text:
module SSGE

using LinearAlgebra
using ForwardDiff
using Distances
using Statistics

include("CovKer.jl")
include("Nystroem.jl")

"""
    SSGEstimator(x_samples, num_eig_func, cov_kernel)

Stores the eigenfunctions `ψ` and the projection matrix `β`, with which
the Stein gradient can be approximated. These are constructed from the points
sampled `x_samples` with `num_eig_func` eigenfunctions regarding the
the covariance kernel `cov_kernel`.
The procedure is described in

Shi, Jiaxin and Sun, Shengyang and Zhu, Jun,
*A Spectral Approach to Gradient Estimation for Implicit Distributions*,
Proceedings of the 35th International Conference on Machine Learning,
4651--4660,
2018

See also: [`(ssge::SSGEstimator)(x)`](@ref)
"""
struct SSGEstimator{F <: Real, V <: AbstractVector{F}, M <: AbstractMatrix{F}, C <: CovKer}
    ψ::Nystroem{F, V, M, C}
    β::M
end

function SSGEstimator(x_samples::M, num_eig_func::Integer, cov_kernel::C) where {
        F <: Real, M <: AbstractMatrix{F}, C <: CovKer}
    ψ = Nystroem(x_samples, num_eig_func, cov_kernel)

    ∇ψ(x) = ForwardDiff.jacobian(ψ, x)

    β =  ∇ψ(x_samples[:,1])
    for i in 2:ψ.num_samples
       β = β + ∇ψ(x_samples[:,i])
    end
    β = -1/ψ.num_samples .* β

    SSGEstimator{F, typeof(ψ.eig_val_inv), M, C}(ψ, β)
end

function SSGEstimator(x_samples::M, r_bar::F) where {
    F <: Real, M <: AbstractMatrix{F}}
    ψ = Nystroem(x_samples, r_bar)

    ∇ψ(x) = ForwardDiff.jacobian(ψ, x)

    β =  ∇ψ(x_samples[:,1])
    for i in 2:ψ.num_samples
    β = β + ∇ψ(x_samples[:,i])
    end
    β = -1/ψ.num_samples .* β

    SSGEstimator{F, typeof(ψ.eig_val_inv), M, typeof(ψ.cov_kernel)}(ψ, β)
end

"""
    (ssge::SSGEstimator)(x)

Estimates the corresponding Stein gradient of `ssge` at point `x`.
"""
(ssge::SSGEstimator)(x::V) where{V <: AbstractVector} =
    (ssge.ψ(x)' * ssge.β)'

end # module


text:
#=

  Decomposition of global constraint alldifferent_except_0 in Julia + ConstraintSolver.

  From Global constraint catalogue:
  http://www.emn.fr/x-info/sdemasse/gccat/Calldifferent_except_0.html
  """
  Enforce all variables of the collection VARIABLES to take distinct
  values, except those variables that are assigned to 0.

  Example
     (<5, 0, 1, 9, 0, 3>)

  The alldifferent_except_0 constraint holds since all the values
  (that are different from 0) 5, 1, 9 and 3 are distinct.
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc
const CS = ConstraintSolver
include("constraints_utils.jl")

# Note: The all_different_except_c is defined in constraints_utils.jl

function all_different_except_0_test(n=5)
    println("n:$n")
    model = Model(optimizer_with_attributes(CS.Optimizer,
                                            "all_solutions"=>true,
                                            # "all_solutions"=>false,
                                            "logging"=>[],
                                            "time_limit"=>10,
                                            )
                                            )
    @variable(model, 0 <=  x[1:n] <= n, Int)

    bs = all_different_except_c(model,x,0)
    increasing(model, x)

    println("solve")
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL

        num_sols = MOI.get(model, MOI.ResultCount())
        println("\nnum_sols:$num_sols\n")

        for sol in 1:num_sols
            println("solution #$sol")
            xx = convert.(Integer,JuMP.value.(x,result=sol))
            # bss = convert.(Integer,JuMP.value.(bs,result=sol))
            println("x:$xx")
            # println("bs:$bss\n")

        end

    end

end

all_different_except_0_test()

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
################################################################################
# Mutable simplex where vertices are represented by some type of abstract vector.
################################################################################

mutable struct MutableSimplex{D, T} <: AbstractMutableSimplex{D, T}
    vertices::Vector{Vector{T}}

    """
    MutableSimplex(pts::Vector{<:AbstractVector{T}})

    Construct a simplex from a vector of vectors.
    """
    function MutableSimplex(pts::Vector{<:AbstractVector{T}}) where {T}
        if !(length(pts) == length(pts[1]) + 1)
            err = """ The input cannot be converted to a simplex.
            Vertices need to have `dim` elements, and there needs to be `dim + 1` vertices.
            """
            throw(DomainError(pts, err))    end
        dim = length(pts[1])
        new{dim, T}([pts[i] for i = 1:length(pts)])
    end

    """
    MutableSimplex(pts::AbstractArray{T, 2}) where T

    Construct a simplex from an array of size `(dim + 1)-by-(dim)` or
    `(dim)-by-(dim + 1)` (faster).
    """
    function MutableSimplex(pts::AbstractArray{T, 2}) where {T}
        s = size(pts)

        if (maximum(s) - minimum(s)) != 1
            err = """ The input cannot be converted to a simplex.
                    size(pts) must be (dim, dim + 1) or (dim + 1, dim).
                """
                throw(DomainError(pts, err))
        end

        if s[1] > s[2] # Rows are points
            dim = s[2]
            return new{dim, T}([pts[i, :] for i = 1:maximum(s)])
        end

        if s[2] > s[1] # Columns are points
            dim = s[1]
            return new{dim, T}([pts[:, i] for i = 1:maximum(s)])
        end
    end
end

# Overwriting the i-th vertex
function Base.setindex!(simplex::MutableSimplex, v, i)
    simplex[i] .= v
end

# Overwriting elements of the i-th vertex
function Base.setindex!(simplex::MutableSimplex, v, i::Int, j)
    if length(v) != length(j)
        err = """
            Trying to overwrite elements $j of vertex $i with $v, which does not
            have the same number of elements as the target.
        """
        throw(ArgumentError(err))
    end
    simplex[i][j] = v
end

export MutableSimplex


text:
#args
dataset = "3dmnist"
root = "./$(dataset)/"
K = 10 # k nearest-neighbors
batch_size = 32
num_classes = 10 #possible values {10,40}
npoints = 1024
epochs = 50

max_stars_repo_name: nirmal-suthar/DGCNN

max_issues_repo_name: nirmal-suthar/DGCNN

max_forks_repo_name: nirmal-suthar/DGCNN


text:
using Clang

error("$(@__FILE__) is broken until Clang.jl supports wrapping C++")

# HIP_HEADERS are those headers to be wrapped.
const HIP_INCLUDE = joinpath(@__DIR__, "include") |> normpath

hip_headers_base = joinpath(@__DIR__, "include", "hip")
const HIP_HEADERS = String[]
push!(HIP_HEADERS, joinpath(hip_headers_base, "julia_hip_runtime.h"))

wc = init(; headers = HIP_HEADERS,
            output_file = joinpath(@__DIR__, "hip_interface.jl"),
            common_file = joinpath(@__DIR__, "hip_types.jl"),
            clang_includes = vcat(HIP_INCLUDE, CLANG_INCLUDE),
            #clang_args = ["-I", joinpath(HIP_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"libhip-hcc",
            clang_diagnostics = true,
            )

run(wc)

# Replace stuff
#=
for (src,dst) in [("libhip-hcc","\"libhip-hcc\"")]
  for file in ("hip_interface.jl", "hip_types.jl")
    run(Cmd(["sed", "-i", "s/$src/$dst/g", file]))
  end
end
=#
# Remove useless HIP_API macros
#run(Cmd(["sed", "-i", "/^.*HIP_API.*\$/d", "hip_types.jl"]))


text:
isspace(x, y, z = 1364) = iseven(count_ones(x^2 + 3x + 2x * y + y + y^2 + z))

neighbours(c::CartesianIndex) = (c + n for n in (
    CartesianIndex(0, 1),
    CartesianIndex(0, -1),
    CartesianIndex(1, 0),
    CartesianIndex(-1, 0)
))

function visitavailable(magic = 1364)
    queue = [CartesianIndex(1, 1) => 0]
    visited = Dict{CartesianIndex,Int}()
    while !isempty(queue)
        coord, visits = popfirst!(queue)
        if coord[1] >= 0 &&
           coord[2] >= 0 &&
           isspace(coord[1], coord[2], magic) &&
           (!haskey(visited, coord) || visited[coord] > visits)
            visited[coord] = visits
            push!(queue, (n => visits + 1 for n in neighbours(coord))...)
        end
    end
    visited
end
part1(x, y, magic = 1364) = visitavailable(magic)[CartesianIndex(x, y)]

@assert part1(7, 4, 10) == 11
part1(31, 39)

part2() = count(<=(50), values(visitavailable(1364)))
part2()

function prettyprint(magic::Integer = 1364)
    visited = visitavailable(magic)
    maxi = maximum(keys(visited))
    for j in 0:maxi[2]
        for i in 0:maxi[1]
            coord = CartesianIndex(i, j)
            if coord == CartesianIndex(1, 1)
                printstyled('X', color = :blue)
            elseif haskey(visited, coord)
                printstyled('O', color = visited[coord] <= 50 ? :green : :light_red)
            elseif isspace(i, j, magic)
                printstyled('.', color = :light_black)
            else
                print('#')
            end
        end
        println()
    end
end
# prettyprint()
# prettyprint(10)

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
module FilePaths

using MacroTools
using Reexport
using Requires

@reexport using FilePathsBase

include("compat.jl")

function __init__()
    @require Glob="c27321d9-0574-5035-807b-f59d2c89b15c" include("glob.jl")
    @require URIParser="30578b45-9adc-5946-b283-645ec420af67" include("uri.jl")
end

end # end of module


text:
using MLBasePlotting
using Base.Test

# plotperf
labels = round(Int, bitrand(1000))
values = randn(1000)
plotperf(labels, values)
plotperf(labels, values, name="ROC Test", curveType="roc")
plotperf(Dict(
	"model1" => (labels, values),
	"model2" => (labels, randn(1000))
))

# area_under_pr
x = rand([true,false], 100)
@test_approx_eq area_under_pr(x, x) 1.0
@test area_under_pr(x, -x) < 0.5

# area_under_roc
# @test_approx_eq area_under_roc(x, x) 1.0 TODO: this is broken and I don't know why yet
@test area_under_roc(x, -x) < 0.5


text:
# # Comparing several models for water waves
#
#md # [`notebook`](@__NBVIEWER_ROOT_URL__notebooks/WWvsXX.ipynb)
#
using WaterWaves1D;
using Test
using TimerOutputs
import Base.Threads: @threads, @sync, @async, @spawn, nthreads, threadid

function run_models()

    reset_timer!()
    #---- parameters
    param = ( μ  = .1,
    			ϵ  = 1,
            	N  = 2^11, 	# number of collocation points
                L  = 10,	# size of the mesh (-L,L)
                T  = 5,		# final time of computation
                dt = 0.001, # timestep
    				);
    #---- initial data
    g(x) = exp.(-abs.(x).^4);
    z(x) = 0*exp.(-x.^2);
    init = Init(g,z);

    #---- models to compare
    models=AbstractModel[]
    push!(models,WaterWaves(param))
    push!(models,PseudoSpectral(param;order=2,dealias=1,lowpass=1/100))
    push!(models,PseudoSpectral(param;order=3,dealias=1,lowpass=1/100))

    #---- Initialize
    function Problems(models)
    	problems = Problem[]
    	for model in models
    		push!(problems, Problem(model, init, param) )
    	end
    	return problems
    end

    #---- Method 1
    @timeit "Method 1 with threads" begin
        problems=Problems(models)
        solve!(problems)
        p1=(last(problems[1].data.U), last(problems[2].data.U) ,last(problems[3].data.U))
    end

    #---- Method 2
    @timeit "Method 2 with threads" begin
        problems=Problems(models)
        @sync for problem in problems
        	@spawn solve!( problem )
        end
        p2=(last(problems[1].data.U), last(problems[2].data.U) ,last(problems[3].data.U))
    end

    #---- Method 3
    @timeit "Consecutive solving" begin
        problems=Problems(models)
        for problem in problems
        	solve!( problem )
        end
        p3=(last(problems[1].data.U), last(problems[2].data.U) ,last(problems[3].data.U))
    end

    #---- Tests
    @testset "Final" begin
    	@test p1==p2
    	@test p2==p3
    end

    print_timer()

end

run_models()


text:
module PuMac

    include("funciones.jl")

    export unir

end


text:
struct Storage{T,N}
    x::Vector{Vector{SVector{3,T}}}
    q::Vector{Vector{UnitQuaternion{T}}}
    v::Vector{Vector{SVector{3,T}}}
    ω::Vector{Vector{SVector{3,T}}}

    function Storage{T}(steps, nbodies) where T
        x = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        q = [[one(UnitQuaternion{T}) for i = steps] for j = 1:nbodies]
        v = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        ω = [[szeros(T, 3) for i = steps] for j = 1:nbodies]
        new{T,length(steps)}(x, q, v, ω)
    end

    Storage{T}() where T = Storage{T}(Base.OneTo(0),0)
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, storage::Storage{T,N}) where {T,N}
    summary(io, storage)
end


text:
module MP3
include("AudioDisplay.jl")

using Plots
using FFTW
using WAV

export AB1A1aPlay
export AB1Aufgabe1a
export AB1A1bPlay
export AB1Aufgabe1b
export AB1A1cPlay
export AB1Aufgabe1c
export AB1A1dPlay
export AB1Aufgabe1d
export AB1A2a1Play
export AB1A2a2Play
export AB1Aufgabe2a
export AB1A2bPlay
export AB1Aufgabe2b

function PlaySinusSound(frequency, amplitude, name)
    #Generate and play a sinus sound.
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = amplitude * sin(2*pi*frequency*x[i]);
    end
    #display(y)
    wavwrite(y, string(name, ".wav"), Fs=44100);
    #y, fs = wavread("SinustonAB1A1a.wav");
    #wavplay(y, fs)
	inline_audioplayer(string(name, ".wav"))
    plot(y)
end

function AB1A1aPlay()
    PlaySinusSound(200, 1, "SinustonAB1A1a")

    # Display 200 Hz sin sound.
    amplitude = 1
    frequenz = 200

    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = amplitude*sin(2*pi*frequenz*x[i]);
    end
    labels = ["Sinuston"];
    plot(x,y, label = labels)
end

function AB1Aufgabe1a(amplitude)
    # Check the amplitude was provided correctly
    if amplitude!=1
        @warn "Amplitude falsch!"
    else
        println("Richtig")
    end
end

function AB1A1bPlay(amplitude)
    PlaySinusSound(200, amplitude, "SinustonAB1A1b")
end

function AB1Aufgabe1b(amplitude)
    if amplitude == 1
        println("Probier doch auch mal eine andere Amplitude.")
    end
    # Display 200 Hz sin sound.
    frequenz = 200

    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = amplitude*sin(2*pi*frequenz*x[i]);
    end
    labels = ["Sinuston"];
    plot(x,y, label = labels)
end

function AB1A1cPlay(A)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    #display(ton)
    wavwrite(y, "SinustonAB1A1c1.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A1c1.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A1c1.wav")
end

function AB1Aufgabe1c(A)
    if A(0) != 0
		@warn "Die Amplitude hat zur Zeit t = 0 nicht den Wert 0!"
    end
    m1 = A(1)-A(0);
    m2 = A(2)-A(1);
    if m1 != m2
		@warn "Die Funktion A besitzt kein lineares Wachstum!"
    end
    labels = ["dein Ton"];
    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    plot(x,y, label = labels)
end

function AB1A1dPlay(A)
    eps = 1/44100
    x = range(0.0,stop=1,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A1d.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A1d.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A1d.wav")
end

function AB1Aufgabe1d(A)
    if A(6) != 0
		@warn "Bei t = 1 hat die Amplitude nicht den Wert 0!"
    end
    m1 = A(1)-A(0);
    m2 = A(2)-A(1);
    if abs(m1 - m2) > 0.0000000001
		@warn "Die Funktion A besitzt kein lineares Abklingverhalten!"
    end
    labels = ["dein Ton"];
    x = 0:0.0001:6;
    y = zeros(length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    plot(x,y, label = labels)
end

function AB1A2a1Play()
    PlaySinusSound(500, 1, "SinustonAB1A2a")

      # Display 500 Hz sin sound.
      amplitude = 1
      frequenz = 500

      x = 0:0.0001:0.1;
      y = zeros(1001);
      for i = 1:1001
          y[i] = amplitude*sin(2*pi*frequenz*x[i]);
      end
      labels = ["Sinuston"];
      plot(x,y, label = labels)
end

function AB1A2a2Play(F_hoeher,F_tiefer)
    eps = 1/44100
    x = range(0.0,stop=2,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(2*pi*200*x[i]);
    end
    z = zeros(Float64,length(x));
    for i = 1:length(z)
        z[i] = sin(F_hoeher*2*pi*200*x[i]);
    end
    w = zeros(Float64,length(x));
    for i = 1:length(w)
        w[i] = sin(F_tiefer*2*pi*200*x[i]);
    end
    ton = [y; z; w];
    #display(ton)
    wavwrite(ton, "SinustonAB1A2a.wav", Fs=44100);
    #ton, fs = wavread("SinustonAB1A2a.wav");
    #wavplay(ton, fs)
	inline_audioplayer("SinustonAB1A2a.wav")
end

function AB1Aufgabe2a(F_hoeher,F_tiefer)
    if F_hoeher >= 6 || F_hoeher <= 0
		@warn "F_hoeher liegt nicht zwischen 0 und 6!"
    end
    if F_tiefer >= 6 || F_tiefer <= 0
		@warn "F_tiefer liegt nicht zwischen 0 und 6!"
    end
    if F_hoeher <= 1
		@warn "Der erste Ton ist nicht höher als der Referenzton!"
    end
    if F_tiefer >= 1
		@warn "Der zweite Ton ist nicht tiefer als der Referenzton!"
    end
    x = 0:0.00001:0.01;
    y = zeros(length(x),3);
    for i = 1:length(x)
        y[i,1] = sin(2*pi*200*x[i]);
        y[i,2] = sin(F_hoeher*2*pi*200*x[i]);
        y[i,3] = sin(F_tiefer*2*pi*200*x[i]);
    end
    labels = ["Referenz" "hoeher" "tiefer"];
    plot(x,y, label = labels)
end

function AB1A2bPlay(F)
    eps = 1/44100
    x = range(0.0,stop=1,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A2b.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A2b.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A2b.wav")
end

function AB1Aufgabe2b(F)
    if F(0) != 1
		@warn "Zum Zeitpunkt t = 0 hat die Frequenz von g nicht den Wert 200 Hz!"
    end
    if F(0.1) != 2
		@warn "Zum Zeitpunkt t = 0.1 hat die Frequenz von g nicht den Wert 400Hz!"
    end
    m1 = F(1)-F(0);
    m2 = F(2)-F(1);
    if abs(m1-m2) > 0.00000001
		@warn "Die Funktion F besitzt kein lineares Wachstum!"
    end
    labels = ["dein Signal"];
    x = 0:0.0001:0.1;
    y = zeros(1001);
    for i = 1:1001
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    plot(x,y, label = labels)
end

function AB1ZAPlay(A,a,b)
    eps = 1/44100
    a_ton = a - 1;
    b_ton = b + 1;
    x = range(a_ton,stop=b_ton,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1A2b.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1A2b.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1A2b.wav")
end

function AB1ZusatzaufgabeA(A4,t_min,a,b)
    if t_min <= a || t_min >= b
		@warn " t_min liegt nicht zwischen a und b!"
    end
    min = 1;
    for i = -10:10
        if A4(t_min + i/100) < A4(t_min)
            min = 0;
        end
    end
    if min == 0
		@warn "Bei t = t_min ist die Amplitude nicht minimal!"
    end
    labels = ["dein Ton"];
    if t_min > a && t_min < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = A4(x[i])*sin(2*pi*200*x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB1ZBPlay(A,a,b)
    eps = 1/44100
    x = range(a,stop=b,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = A(x[i])*sin(2*pi*50*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1ZB.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1ZB.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1ZB.wav")
end

function AB1ZusatzaufgabeB(A5,t_max,a,b)
    if t_max <= a || t_max >= b
		@warn "t_max liegt nicht zwischen a und b!"
    end
    if A5(a) != 0
		@warn "Bei t = a hat die Amplitude nicht den Wert 0!"
    end
    if A5(b) != 0
		@warn "Bei t = b hat die Amplitude nicht den Wert 0!"
    end
    max = 1;
    for i = -10:10
        if A5(t_max + i/100) > A5(t_max)
            max = 0;
        end
    end
    if max == 0
		@warn "Bei t = t_max ist die Amplitude nicht maximal!"
    end
    labels = ["dein Ton"];
    if t_max > a && t_max < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = A5(x[i])*sin(2*pi*50x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB1ZCPlay(a,F)
    eps = 1/44100
    b = a + 1;
    x = range(a,stop=b,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1ZC.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1ZC.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1ZC.wav")
end

function AB1ZusatzaufgabeC(a,b,F)
    if b <= a
		@warn "b muss größer sein als a!"
    end
    increase = 1;
    for i = 1:100
        if F(a+(b-a)/100*i) < F(a+(b-a)/100*(i+1))
            increase = 0;
        end
    end
    if increase == 0
		@warn "Die Frequenz fällt nicht auf dem Intervall [a,b]!"
    end
    labels = ["dein Signal"];
    if a < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = sin(F(x[i])*2*pi*200*x[i]);
        end
        plot(x,y, label = labels)
    end
end

function AB1ZDPlay(a,b,F)
    eps = 1/44100
    x = range(a,stop=b,step=eps);
    y = zeros(Float64,length(x));
    for i = 1:length(y)
        y[i] = sin(F(x[i])*2*pi*200*x[i]);
    end
    #display(y)
    wavwrite(y, "SinustonAB1ZD.wav", Fs=44100);
    #y, fs = wavread("SinustonAB1ZD.wav");
    #wavplay(y, fs)
	inline_audioplayer("SinustonAB1ZD.wav")
end

function AB1ZusatzaufgabeD(a,b,F)
    if b <= a
		@warn "b muss größer sein als a!"
    end
    labels = ["dein Signal"];
    if a < b
        x = a:(b-a)/10000:b;
        y = zeros(length(x));
        for i = 1:length(y)
            y[i] = sin(F(x[i])*2*pi*200*x[i]);
        end
        plot(x,y, label = labels)
    end
end

end


text:
module TestTensorNetworks
using Test

using LinearAlgebra: tr

using Catlab.CategoricalAlgebra.CSets, Catlab.WiringDiagrams
using Catlab.LinearAlgebra.TensorNetworks

# Parsing and code generation
#############################

# Parsing
#--------

parsed = @tensor_network (i,j,k,ℓ) D[i,j,k] = A[i,ℓ] * B[j,ℓ] * C[k,ℓ]
d = RelationDiagram{Symbol}(3)
add_box!(d, 2, name=:A); add_box!(d, 2, name=:B); add_box!(d, 2, name=:C)
add_junctions!(d, 4, variable=[:i,:j,:k,:ℓ])
set_junction!(d, [1,4,2,4,3,4])
set_junction!(d, 1:3, outer=true)
@test parsed == d

parsed = @tensor_network D[i,j,k] = A[i,ℓ] * B[j,ℓ] * C[k,ℓ]
@test parsed == d
parsed = @tensor_network D[i,j,k] := A[i,ℓ] * B[j,ℓ] * C[k,ℓ]
@test parsed == d

# Degenerate case: single term.
parsed = @tensor_network B[i,j,k] = A[i,j,k]
d = singleton_diagram(RelationDiagram{Symbol}, 3, name=:A)
set_subpart!(d, :variable, [:i,:j,:k])
@test parsed == d

# Degenerate case: no terms.
parsed = @tensor_network A[i,j] = 1
d = RelationDiagram{Symbol}(2)
add_junctions!(d, 2, variable=[:i,:j])
set_junction!(d, 1:2, outer=true)
@test parsed == d

# Round trip: parse then compile
#-------------------------------

macro roundtrip_tensor(tensor)
  quote
    compiled = gen_tensor_notation(@tensor_network($tensor),
                                   assign_op=:(=), assign_name=:out)
    @test compiled == $(QuoteNode(tensor))
  end
end

@roundtrip_tensor out[i,k] = A[i,j] * B[j,k]
@roundtrip_tensor out[i,j,k] = A[i,ℓ] * B[j,ℓ] * C[k,ℓ]
@roundtrip_tensor out = u[i] * A[i,j] * v[j]
@roundtrip_tensor out[j] = α * a[i] * B[i,j]

# Evaluation
############

# Binary matrix multiply.
A, B = randn((3,4)), randn((4,5))
network = @tensor_network C[i,k] = A[i,j] * B[j,k]
@test oapply(network, [A,B]) ≈ A*B

# Ternary matrix multiply.
C = randn((5,6))
network = @tensor_network D[i,ℓ] = A[i,j] * B[j,k] * C[k,ℓ]
@test oapply(network, [A,B,C]) ≈ A*B*C

# Trace.
A = randn((5,5))
network = @tensor_network out[] = A[i,i]
@test oapply(network, [A])[] ≈ tr(A)

# Tensor product.
A, B = randn((3,4)), randn((5,6))
network = @tensor_network C[i,j,k,ℓ] = A[i,j] * B[k,ℓ]
@test oapply(network, [A,B]) ≈ (reshape(A, (3,4,1,1)) .* reshape(B, (1,1,5,6)))

end


text:
using GeometricEquations
using GeometricEquations: symplectic_matrix
using Test

include("functions.jl")
include("initial_conditions.jl")

@testset "$(rpad("Ordinary Differential Equations (ODE)",80))" begin

    ode  = ODE(ode_v, NullInvariants(), NullParameters(), NullPeriodicity())
    ode1 = ODE(ode_v)
    ode2 = ODE(ode_v; invariants=NullInvariants())
    ode3 = ODE(ode_v; parameters=NullParameters())
    ode4 = ODE(ode_v; periodicity=NullPeriodicity())

    @test ode == ode1
    @test ode == ode2
    @test ode == ode3
    @test ode == ode4

    @test hash(ode) == hash(ode1)
    @test hash(ode) == hash(ode2)
    @test hash(ode) == hash(ode3)
    @test hash(ode) == hash(ode4)

    @test functions(ode) == NamedTuple{(:v,)}((ode_v,))
    @test solutions(ode) == NamedTuple()

    @test parameters(ode) == NullParameters()
    @test invariants(ode) == NullInvariants()
    @test periodicity(ode) == NullPeriodicity()

    @test hasvectorfield(ode) == true
    @test hassolution(ode) == false
    @test hasprimary(ode) == false
    @test hassecondary(ode) == false

    @test hasinvariants(ode) == false
    @test hasparameters(ode) == false
    @test hasperiodicity(ode) == false

    @test hashamiltonian(ode) == false
    @test haslagrangian(ode) == false

    funcs = functions(ode)

    @test_nowarn funcs.v(t₀, x₀, zero(x₀), NullParameters())

    funcs = functions(ode, NullParameters())

    @test_nowarn funcs.v(t₀, x₀, zero(x₀))

end

@testset "$(rpad("Split Ordinary Differential Equations (SODE)",80))" begin

    sode  = SODE(sode_v, nothing, NullInvariants(), NullParameters(), NullPeriodicity())
    sode1 = SODE(sode_v)
    sode2 = SODE(sode_v; invariants=NullInvariants())
    sode3 = SODE(sode_v; parameters=NullParameters())
    sode4 = SODE(sode_v; periodicity=NullPeriodicity())

    @test sode == sode1
    @test sode == sode2
    @test sode == sode3
    @test sode == sode4

    @test hash(sode) == hash(sode1)
    @test hash(sode) == hash(sode2)
    @test hash(sode) == hash(sode3)
    @test hash(sode) == hash(sode4)

    @test functions(sode) == NamedTuple{(:v,)}((sode_v,))
    @test solutions(sode) == NamedTuple()

    @test parameters(sode) == NullParameters()
    @test invariants(sode) == NullInvariants()
    @test periodicity(sode) == NullPeriodicity()

    @test hasvectorfield(sode) == true
    @test hassolution(sode) == false
    @test hasprimary(sode) == false
    @test hassecondary(sode) == false

    @test hasinvariants(sode) == false
    @test hasparameters(sode) == false
    @test hasperiodicity(sode) == false

    @test hashamiltonian(sode) == false
    @test haslagrangian(sode) == false

    sode  = SODE(sode_v, sode_q, NullInvariants(), NullParameters(), NullPeriodicity())
    sode1 = SODE(sode_v, sode_q)
    sode2 = SODE(sode_v, sode_q; invariants=NullInvariants())
    sode3 = SODE(sode_v, sode_q; parameters=NullParameters())
    sode4 = SODE(sode_v, sode_q; periodicity=NullPeriodicity())

    @test sode == sode1
    @test sode == sode2
    @test sode == sode3
    @test sode == sode4

    @test hash(sode) == hash(sode1)
    @test hash(sode) == hash(sode2)
    @test hash(sode) == hash(sode3)
    @test hash(sode) == hash(sode4)

    @test functions(sode) == NamedTuple{(:v,)}((sode_v,))
    @test solutions(sode) == NamedTuple{(:q,)}((sode_q,))

    @test parameters(sode) == NullParameters()
    @test invariants(sode) == NullInvariants()
    @test periodicity(sode) == NullPeriodicity()

    @test hasvectorfield(sode) == true
    @test hassolution(sode) == true
    @test hasprimary(sode) == false
    @test hassecondary(sode) == false

    @test hasinvariants(sode) == false
    @test hasparameters(sode) == false
    @test hasperiodicity(sode) == false

    @test hashamiltonian(sode) == false
    @test haslagrangian(sode) == false

    sode  = SODE(nothing, sode_q, NullInvariants(), NullParameters(), NullPeriodicity())
    sode1 = SODE(nothing, sode_q)
    sode2 = SODE(nothing, sode_q; invariants=NullInvariants())
    sode3 = SODE(nothing, sode_q; parameters=NullParameters())
    sode4 = SODE(nothing, sode_q; periodicity=NullPeriodicity())

    @test sode == sode1
    @test sode == sode2
    @test sode == sode3
    @test sode == sode4

    @test hash(sode) == hash(sode1)
    @test hash(sode) == hash(sode2)
    @test hash(sode) == hash(sode3)
    @test hash(sode) == hash(sode4)

    @test functions(sode) == NamedTuple()
    @test solutions(sode) == NamedTuple{(:q,)}((sode_q,))

    @test parameters(sode) == NullParameters()
    @test invariants(sode) == NullInvariants()
    @test periodicity(sode) == NullPeriodicity()

    @test hasvectorfield(sode) == false
    @test hassolution(sode) == true
    @test hasprimary(sode) == false
    @test hassecondary(sode) == false

    @test hasinvariants(sode) == false
    @test hasparameters(sode) == false
    @test hasperiodicity(sode) == false

    @test hashamiltonian(sode) == false
    @test haslagrangian(sode) == false

end

@testset "$(rpad("Partitioned Ordinary Differential Equations (PODE)",80))" begin

    pode  = PODE(pode_eqs..., NullInvariants(), NullParameters(), NullPeriodicity())
    pode1 = PODE(pode_eqs...)
    pode2 = PODE(pode_eqs...; invariants=NullInvariants())
    pode3 = PODE(pode_eqs...; parameters=NullParameters())
    pode4 = PODE(pode_eqs...; periodicity=NullPeriodicity())

    @test pode == pode1
    @test pode == pode2
    @test pode == pode3
    @test pode == pode4

    @test hash(pode) == hash(pode1)
    @test hash(pode) == hash(pode2)
    @test hash(pode) == hash(pode3)
    @test hash(pode) == hash(pode4)

    @test functions(pode) == NamedTuple{(:v,:f)}(pode_eqs)
    @test solutions(pode) == NamedTuple()

    @test parameters(pode) == NullParameters()
    @test invariants(pode) == NullInvariants()
    @test periodicity(pode) == NullPeriodicity()

    @test hasvectorfield(pode) == true
    @test hassolution(pode) == false
    @test hasprimary(pode) == false
    @test hassecondary(pode) == false

    @test hasinvariants(pode) == false
    @test hasparameters(pode) == false
    @test hasperiodicity(pode) == false

    @test hashamiltonian(pode) == false
    @test haslagrangian(pode) == false

    funcs = functions(pode)

    @test_nowarn funcs.v(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(p₀), NullParameters())

    funcs = functions(pode, NullParameters())

    @test_nowarn funcs.v(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(p₀))

    # funcs = functions(pode)
    # @test funcs.v == pode_v == pode.v
    # @test funcs.f == pode_f == pode.f

    # rode = ODE(ode_v, t₀, [x₀])
    # code = convert(ODE, pode)
    # v₁ = zero(x₀)
    # v₂ = zero(x₀)
    # rode.v(rode.t₀, rode.q₀[begin], v₁)
    # code.v(code.t₀, code.q₀[begin], v₂)
    # @test v₁ == v₂

    # rode = SODE((v_sode_1, v_sode_2), t₀, [x₀])
    # code = convert(SODE, pode)
    # v₁ = zero(x₀)
    # v₂ = zero(x₀)
    # rode.v[1](rode.t₀, rode.q₀[begin], v₁)
    # code.v[1](code.t₀, code.q₀[begin], v₂)
    # @test v₁ == v₂
    # rode.v[2](rode.t₀, rode.q₀[begin], v₁)
    # code.v[2](code.t₀, code.q₀[begin], v₂)
    # @test v₁ == v₂

end

@testset "$(rpad("Implicit Ordinary Differential Equations (IODE)",80))" begin

    iode  = IODE(iode_eqs..., iode_v, iode_f, NullInvariants(), NullParameters(), NullPeriodicity())

    iode1 = IODE(iode_eqs...; v̄=iode_v)
    iode2 = IODE(iode_eqs...; v̄=iode_v, invariants=NullInvariants())
    iode3 = IODE(iode_eqs...; v̄=iode_v, parameters=NullParameters())
    iode4 = IODE(iode_eqs...; v̄=iode_v, periodicity=NullPeriodicity())
    iode5 = IODE(iode_eqs...; v̄=iode_v, f̄=iode_f)
    iode6 = IODE(iode_eqs...; v̄=iode_v, f̄=iode_f, invariants=NullInvariants())
    iode7 = IODE(iode_eqs...; v̄=iode_v, f̄=iode_f, parameters=NullParameters())
    iode8 = IODE(iode_eqs...; v̄=iode_v, f̄=iode_f, periodicity=NullPeriodicity())

    @test iode == iode1
    @test iode == iode2
    @test iode == iode3
    @test iode == iode4
    @test iode == iode5
    @test iode == iode6
    @test iode == iode7
    @test iode == iode8

    @test hash(iode) == hash(iode1)
    @test hash(iode) == hash(iode2)
    @test hash(iode) == hash(iode3)
    @test hash(iode) == hash(iode4)
    @test hash(iode) == hash(iode5)
    @test hash(iode) == hash(iode6)
    @test hash(iode) == hash(iode7)
    @test hash(iode) == hash(iode8)

    # @test functions(iode) == NamedTuple{(:ϑ,:f,:g,:v̄,:f̄)}(iode_eqs)
    @test solutions(iode) == NamedTuple()

    @test parameters(iode) == NullParameters()
    @test invariants(iode) == NullInvariants()
    @test periodicity(iode) == NullPeriodicity()

    @test hasvectorfield(iode) == true
    @test hassolution(iode) == false
    @test hasprimary(iode) == false
    @test hassecondary(iode) == false

    @test hasinvariants(iode) == false
    @test hasparameters(iode) == false
    @test hasperiodicity(iode) == false

    @test hashamiltonian(iode) == false
    @test haslagrangian(iode) == false

    funcs = functions(iode)

    @test_nowarn funcs.ϑ(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.g(t₀, q₀, p₀, zero(q₀), NullParameters())

    funcs = functions(iode, NullParameters())

    @test_nowarn funcs.ϑ(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.g(t₀, q₀, p₀, zero(q₀))

    # @test periodicity(iode) == zero(q₀)
    # @test initial_conditions(iode) == (t₀, [q₀], [p₀], [λ₀])

    # funcs = functions(iode)
    # @test funcs.ϑ == iode_ϑ == iode.ϑ
    # @test funcs.f == iode_f == iode.f
    # @test funcs.g == iode_g == iode.g
    # @test funcs.v̄ == iode_v == iode.v̄
    # @test funcs.f̄ == iode_f == iode.f̄

    # @test iode == similar(iode, t₀, q₀, p₀, λ₀)
    # @test iode == similar(iode, t₀, q₀, p₀)
    # @test iode == similar(iode, q₀, p₀)

end

@testset "$(rpad("Hamiltonian Ordinary Differential Equations (HODE)",80))" begin

    hode_eqs = (pode_v, pode_f, hode_ω, hode_h)

    hode  = HODE(pode_v, pode_f, hode_ω, hode_h, NullInvariants(), NullParameters(), NullPeriodicity())

    hode1 = HODE(hode_eqs...)
    hode2 = HODE(hode_eqs...; invariants=NullInvariants())
    hode3 = HODE(hode_eqs...; parameters=NullParameters())
    hode4 = HODE(hode_eqs...; periodicity=NullPeriodicity())

    @test hode == hode1
    @test hode == hode2
    @test hode == hode3
    @test hode == hode4

    @test hash(hode) == hash(hode1)
    @test hash(hode) == hash(hode2)
    @test hash(hode) == hash(hode3)
    @test hash(hode) == hash(hode4)

    @test functions(hode) == NamedTuple{(:v,:f,:poisson,:h)}((pode_v, pode_f, hode_ω, hode_h))
    @test solutions(hode) == NamedTuple()

    @test parameters(hode) == NullParameters()
    @test invariants(hode) == NullInvariants()
    @test periodicity(hode) == NullPeriodicity()

    @test hasvectorfield(hode) == true
    @test hassolution(hode) == false
    @test hasprimary(hode) == false
    @test hassecondary(hode) == false

    @test hasinvariants(hode) == false
    @test hasparameters(hode) == false
    @test hasperiodicity(hode) == false

    @test hashamiltonian(hode) == true
    @test haslagrangian(hode) == false

    funcs = functions(hode)

    @test_nowarn funcs.v(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(p₀), NullParameters())
    @test_nowarn funcs.h(t₀, q₀, p₀, NullParameters())
    @test_nowarn funcs.poisson(t₀, q₀, p₀, zeros(2,2), NullParameters())

    funcs = functions(hode, NullParameters())

    @test_nowarn funcs.v(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(p₀))
    @test_nowarn funcs.h(t₀, q₀, p₀)
    @test_nowarn funcs.poisson(t₀, q₀, p₀, zeros(2,2))

    # funcs = functions(hode)
    # @test funcs.v == pode_v == hode.v
    # @test funcs.f == pode_f == hode.f
    # @test funcs.h == pode_h == hode.hamiltonian

    # @test hode == similar(hode, t₀, q₀, p₀)
    # @test hode == similar(hode, q₀, p₀)

    # rode = ODE(ode_v, t₀, [x₀])
    # code = convert(ODE, hode)
    # v₁ = zero(x₀)
    # v₂ = zero(x₀)
    # rode.v(rode.t₀, rode.q₀[begin], v₁)
    # code.v(code.t₀, code.q₀[begin], v₂)
    # @test v₁ == v₂

    # pode = convert(PODE, hode)
    # v₁ = zero(q₀)
    # v₂ = zero(q₀)
    # f₁ = zero(p₀)
    # f₂ = zero(p₀)
    # hode.v(hode.t₀, hode.q₀[begin], hode.p₀[begin], v₁)
    # pode.v(pode.t₀, pode.q₀[begin], pode.p₀[begin], v₂)
    # hode.f(hode.t₀, hode.q₀[begin], hode.p₀[begin], f₁)
    # pode.f(pode.t₀, pode.q₀[begin], pode.p₀[begin], f₂)
    # @test v₁ == v₂
    # @test f₁ == f₂

    # rode = SODE((v_sode_1, v_sode_2), t₀, [x₀])
    # code = convert(SODE, hode)
    # v₁ = zero(x₀)
    # v₂ = zero(x₀)
    # rode.v[1](rode.t₀, rode.q₀[begin], v₁)
    # code.v[1](code.t₀, code.q₀[begin], v₂)
    # @test v₁ == v₂
    # rode.v[2](rode.t₀, rode.q₀[begin], v₁)
    # code.v[2](code.t₀, code.q₀[begin], v₂)
    # @test v₁ == v₂

end

@testset "$(rpad("Lagrangian Ordinary Differential Equations (LODE)",80))" begin

    lode_eqs = (iode_ϑ, iode_f, iode_g, lode_ω, lode_l)

    lode = LODE(iode_ϑ, iode_f, iode_g, lode_ω, iode_v, iode_f, lode_l, NullInvariants(), NullParameters(), NullPeriodicity())

    lode1 = LODE(lode_eqs..., v̄=iode_v)
    lode2 = LODE(lode_eqs..., v̄=iode_v, invariants=NullInvariants())
    lode3 = LODE(lode_eqs..., v̄=iode_v, parameters=NullParameters())
    lode4 = LODE(lode_eqs..., v̄=iode_v, periodicity=NullPeriodicity())

    @test lode == lode1
    @test lode == lode2
    @test lode == lode3
    @test lode == lode4

    @test hash(lode) == hash(lode1)
    @test hash(lode) == hash(lode2)
    @test hash(lode) == hash(lode3)
    @test hash(lode) == hash(lode4)

    @test functions(lode) == NamedTuple{(:ϑ,:f,:g,:ω,:v̄,:f̄,:l)}((iode_ϑ, iode_f, iode_g, lode_ω, iode_v, iode_f, lode_l))
    @test solutions(lode) == NamedTuple()

    @test parameters(lode) == NullParameters()
    @test invariants(lode) == NullInvariants()
    @test periodicity(lode) == NullPeriodicity()

    @test hasvectorfield(lode) == true
    @test hassolution(lode) == false
    @test hasprimary(lode) == false
    @test hassecondary(lode) == false

    @test hasinvariants(lode) == false
    @test hasparameters(lode) == false
    @test hasperiodicity(lode) == false

    @test hashamiltonian(lode) == false
    @test haslagrangian(lode) == true

    funcs = functions(lode)

    @test_nowarn funcs.ϑ(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.g(t₀, q₀, p₀, zero(q₀), NullParameters())
    @test_nowarn funcs.l(t₀, q₀, p₀, NullParameters())
    @test_nowarn funcs.ω(t₀, q₀, p₀, zeros(2,2), NullParameters())

    funcs = functions(lode, NullParameters())

    @test_nowarn funcs.ϑ(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.f(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.g(t₀, q₀, p₀, zero(q₀))
    @test_nowarn funcs.l(t₀, q₀, p₀)
    @test_nowarn funcs.ω(t₀, q₀, p₀, zeros(2,2))

    # funcs = functions(lode)
    # @test funcs.ϑ == iode_ϑ == lode.ϑ
    # @test funcs.f == iode_f == lode.f
    # @test funcs.g == iode_g == lode.g
    # @test funcs.ω == lode_ω == lode.ω
    # @test funcs.v̄ == iode_v == lode.v̄
    # @test funcs.f̄ == iode_f == lode.f̄
    # @test funcs.l == lode_l == lode.lagrangian

    # @test lode == similar(lode, t₀, q₀, p₀, λ₀)
    # @test lode == similar(lode, t₀, q₀, p₀)
    # @test lode == similar(lode, q₀, p₀)

    # iode = convert(IODE, lode)
    # v₁ = zero(q₀)
    # v₂ = zero(q₀)
    # ϑ₁ = zero(q₀)
    # ϑ₂ = zero(q₀)
    # f₁ = zero(p₀)
    # f₂ = zero(p₀)
    # g₁ = zero(p₀)
    # g₂ = zero(p₀)
    # lode.v̄(lode.t₀, lode.q₀[begin], lode.p₀[begin], v₁)
    # iode.v̄(iode.t₀, iode.q₀[begin], iode.p₀[begin], v₂)
    # lode.ϑ(lode.t₀, lode.q₀[begin], v₁, ϑ₁)
    # iode.ϑ(iode.t₀, iode.q₀[begin], v₂, ϑ₂)
    # lode.f(lode.t₀, lode.q₀[begin], v₁, f₁)
    # iode.f(iode.t₀, iode.q₀[begin], v₂, f₂)
    # lode.g(lode.t₀, lode.q₀[begin], v₁, g₁)
    # iode.g(iode.t₀, iode.q₀[begin], v₂, g₂)
    # @test v₁ == v₂
    # @test ϑ₁ == ϑ₂
    # @test f₁ == f₂
    # @test g₁ == g₂

end


text:
using Test, ReferenceTests, BSON
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil

n = 32 # grid size

# Reference test using ReferenceTests.jl
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2) for (v1,v2) in zip(values(d1), values(d2))])


text:
#
# Merge Sort
#

function mergeSort(A)

  Anew = copy(A)
  mergesort_array(Anew,  A,  1,  length(A)+1)

end

function mergesort_array(A, result, left, right)

  if right - left < 2
    return
  end
  if right - left == 2
    if result[left] > result[left+1]
      result[left],  result[left+1] = result[left+1], result[left]
    end
    return
  end

  mid = Int(round((right + left) / 2))
  mergesort_array(result, A, left,  mid)
  mergesort_array(result, A, mid, right)

  i = left
  j = mid
  idx = left
  while idx < right
    if j >= right || (i < mid && A[i] < A[j])
      result[idx] = A[i]
      i = i + 1
    else
      result[idx] = A[j]
      j = j + 1
    end
    idx = idx + 1
  end

end

function main()

  # list to be sorted
  rng = MersenneTwister(1234);
  A = shuffle(rng, Vector(1:100))
  Aorigin = copy(A)

  # sorting
  @time mergeSort(A)

  print("The original list: ")
  println(Aorigin)

  print("The sorted list: ")
  println(A)

end

if contains(@__FILE__, PROGRAM_FILE)
    main()
end

max_stars_repo_name: szmlb/algorithmsInJulia

max_issues_repo_name: szmlb/algorithmsInJulia

max_forks_repo_name: szmlb/algorithmsInJulia


text:
struct TreeNode{T}
  node::T
  showid::Bool
end

similar_tree_node(n::TreeNode,node) = TreeNode(node,n.showid)

AbstractTrees.children(a::TreeNode) = get_children(a,a.node)

AbstractTrees.printnode(io::IO,a::TreeNode) = print_node(io,a,a.node)

get_children(n::TreeNode,a) = ()

function print_node(io,n::TreeNode,a)
  if length(get_children(n,a)) == 0
    print(io,typeof(a))
    if n.showid
      print(io," (")
      show(io,objectid(a))
      print(io,")")
    end
  else
    print(io,nameof(typeof(a)))
    if n.showid
      print(io," (")
      show(io,objectid(a))
      print(io,")")
    end
  end
end

function print_op_tree(a,args...;kwargs...)
  print_op_tree(stdin,a,args...;kwargs...)
end

function print_op_tree(io::IO,a,args...;showid=false,kwargs...)
  tree = TreeNode(a,showid)
  AbstractTrees.print_tree(io,tree,args...;kwargs...)
end

# Specialization

function get_children(n::TreeNode, a::Fill)
   (similar_tree_node(n,a.value),)
end

function get_children(n::TreeNode, a::CompressedArray)
  (similar_tree_node(n,a.values),similar_tree_node(n,a.ptrs))
end

function get_children(n::TreeNode, a::LazyArray)
  (similar_tree_node(n,a.maps),map(i->similar_tree_node(n,i),a.args)...)
end

function get_children(n::TreeNode, a::SubVector)
  (similar_tree_node(n,a.vector),)
end

function get_children(n::TreeNode, a::Reindex)
  (similar_tree_node(n,a.values),)
end

function get_children(n::TreeNode, a::PosNegReindex)
  (similar_tree_node(n,a.values_pos),similar_tree_node(n,a.values_neg))
end

function get_children(n::TreeNode, a::PosNegPartition)
  (similar_tree_node(n,a.i_to_iposneg),similar_tree_node(n,a.ipos_to_i),similar_tree_node(n,a.ineg_to_i))
end

function get_children(n::TreeNode, a::AppendedArray)
  (similar_tree_node(n,a.a),similar_tree_node(n,a.b))
end

function get_children(n::TreeNode, a::Table)
  (similar_tree_node(n,a.data),similar_tree_node(n,a.ptrs))
end

function get_children(n::TreeNode, a::VectorWithEntryRemoved)
  (similar_tree_node(n,a.a),)
end

function get_children(n::TreeNode, a::VectorWithEntryInserted)
  (similar_tree_node(n,a.a),)
end


text:
abstract type GeometricEntity end

struct Segment <: GeometricEntity
    l::Symbol
    r::Symbol
end
show(io::IO, AB::Segment) = print(io, AB.l, AB.r)
const segment = Segment

struct Circle <: GeometricEntity
    O::Symbol
    A::Symbol
end
show(io::IO, γ::Circle) = print(io, "circle center ", γ.O, " radius ", γ.O, γ.A)
const circle = Circle


text:
struct bbox
    xmin::Float64
    ymin::Float64
    xmax::Float64
    ymax::Float64
end

bbox(;xmin, ymin, xmax, ymax) = bbox(xmin, ymin, xmax, ymax)

"""
    GeoArray{T <: Union{Real,Union{Missing,Real}}} <: AbstractArray{T,3}

A GeoArray is an AbstractArray, an AffineMap for calculating coordinates based
on the axes and a CRS definition to interpret these coordinates into in the real
world. It's three dimensional and can be seen as a stack (3D) of 2D geospatial
rasters (bands), the dimensions are :x, :y, and :bands. The AffineMap and CRS
(coordinates) only operate on the :x and :y dimensions.
"""
mutable struct GeoArray{T<:Union{Real, Union{Missing, Real}}} <: AbstractArray{T, 3}
    A::AbstractArray{T, 3}
    f::AffineMap
    crs::WellKnownText{GeoFormatTypes.CRS, <:String}
end


text:
# using Revise # uncomment for development

using Test
using Random
using StatsFuns
using Optim
using StatsBase: countmap
using LinearAlgebra
using ForwardDiff
using DataFrames
using Distributions
using BenchmarkTools: @btime

using rcmnl

@testset "Gauss Hermite integration" begin
    @test rcmnl.integrate_wrt_normal_pdf(x -> x, 1.23, 2.34; npts=15) ≈ 1.23
    @test sqrt(rcmnl.integrate_wrt_normal_pdf(x -> (x-1.23)^2, 1.23, 2.34; npts=15)) ≈ 2.34

    μ = [1.0, 0.5]
    L = [1.0  0.0; 0.7 0.3]

    @test rcmnl.integrate_wrt_normal_pdf(x -> (x-μ)*(x-μ)', μ, L; npts=7) ≈ L*L'
    @test rcmnl.integrate_wrt_normal_pdf(x -> x,            μ, L; npts=7) ≈ μ
end

Random.seed!(1234)

# dim of data
nindv    = 1_000
nt       = 20
nquadpts = 7

# coefs
β = [ 1.0  0.3;
     -2.0  0.0;
      0.0  1.0;
      0.2  0.5]
k, nchoice = size(β)
@assert nchoice == 2

# cholesky decomposition of Σ
Σchol = [1.0 0.0;
         0.5 0.5]
θtrue = vcat(vec(β), Σchol[:,1], Σchol[end])

# generate X and V
X = randn(nt, nindv, k)
V = randn(nindv,2)*Σchol' # need Σchol b/c usu. assume U is a 2-dim vector

# utilities
U = zeros(3, nt, nindv)
mul!(reshape(view(U,2:3,:,:), 2,:), β', reshape(X, :, k)')
U[2,:,:] .+= V[:,1]'
U[3,:,:] .+= V[:,2]'

# probabilities
PrU = mapslices(softmax, U; dims=1)
cumPrU = cumsum(PrU; dims=1)

# instead of simulating random type-1 extreme values, we just
# use a uniform variable and the CDF
y = Array{Int64}(undef, nt, nindv)
for idx in CartesianIndices(y)
    vw = view(cumPrU, :, idx)
    y[idx] = searchsortedfirst(vw, rand())
end

# how many of each choice?
countmap(vec(y))

# run function once
simlogL(y, X, θtrue; npts=nquadpts)
println("It runs!")

# wrapper for optimizer
f(θ) = simlogL(y, X, θ; npts=nquadpts)

# starting value
theta0 = vcat(vec(β), Σchol[:,1], Σchol[end]) .* 2

# time function
@btime f(theta0)

# profile function in IDE
# @profileview f(theta0)

# optimize
res = optimize(f, theta0, BFGS(), Optim.Options(;show_trace=true), autodiff=:forward)

# get results
thetahat = res.minimizer
H = ForwardDiff.hessian(f, thetahat)
se = sqrt.(diag(inv(H)))
pval = 2*normcdf.(-abs.(thetahat./se))

# nice table
results_to_print = DataFrame(
    θtrue = θtrue,
    θhat = thetahat,
    se = se,
    pval = pval
)
@show results_to_print

# Wald test
alpha_reject = 0.001
wald_test_stat = (thetahat - θtrue)' * (H \ (thetahat - θtrue))
wald_test_p = ccdf(Chisq(length(θtrue)), wald_test_stat)
if wald_test_p < alpha_reject
    println("Reject that we got the right answer")
else
    println("We fail to reject our estimate!!")
end

@test wald_test_p > alpha_reject


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module Unknown

export HereSomewhere, Consumed, Nonsense, Shifted

include("unknowns.jl")

end


text:
function f1(x, y)
  m = 0
  n = 2

  if x > 10
    m = x - 10
    println("OK")
  else
    m = 3 * x
  end

  while y < 3
    y += 1
    n = 2 * y
  end

  d = m + n
  return d
end

res = f1(1, 1) # should be 9
if res < 10
  println("Here1")
else
  println("Here2")
end

s = [1,2,3,4,5,6]

for num in s
  k = num + 1
  println(k)
end

max_stars_repo_name: judy-vscode/judy

max_issues_repo_name: judy-vscode/judy

max_forks_repo_name: judy-vscode/judy


text:
# Dirichlet distribution

export Dirichlet

@parameterized Dirichlet(α)

TV.as(d::Dirichlet{(:α,)}) = TV.UnitSimplex(length(d.α))

function basemeasure(μ::Dirichlet{(:α,)})
    t = as(μ)
    d = dimension(t)
    return Pushforward(t, Lebesgue(ℝ)^d, false)
end

@kwstruct Dirichlet(α)

function logdensity(d::Dirichlet{(:α,)}, x)
    α = d.α
    s = 0.0
    for j in eachindex(x)
        s += α[j] * log(x[j])
    end
    return s
end

Base.rand(rng::AbstractRNG, T::Type, μ::Dirichlet) = rand(rng, Dists.Dirichlet(μ.α))

distproxy(d::Dirichlet{(:α,)}) = Dists.Dirichlet(d.α)

function testvalue(d::Dirichlet{(:α,)})
    n = length(d.α)
    Fill(1/n, n)
end


text:
module Wind

using ..Ahorn, Maple

placements = Maple.Wind

function Ahorn.canFgBg(effect::Maple.Wind)
    return true, true
end

end


text:
# Example script to test the branch replacement feature in multi-period optimization of demand flexibility, AC & DC lines and storage investments

## Import relevant packages

import FlexPlan; const _FP = FlexPlan
import PowerModelsACDC; const _PMACDC = PowerModelsACDC # For DC grid
import PowerModels; const _PM = PowerModels # For AC grid and common functions

# Include sample data from FlexPlan repository; you can of course also use your own data
include("../io/create_profile.jl")

# Add solver packages
# > Note: solver packages are needed to handle communication between the solver and JuMP;
# > the commercial ones do not include the solver itself.
import Cbc
#import Ipopt
#import Juniper
#import SCS
#import Mosek
#import Gurobi

# Solver configuration
cbc = _FP.optimizer_with_attributes(Cbc.Optimizer, "logLevel"=>0)
#ipopt = _FP.optimizer_with_attributes(Ipopt.Optimizer, "tol"=>1e-4, "print_level"=>0)
#juniper = _FP.optimizer_with_attributes(Juniper.Optimizer, "nl_solver"=>ipopt, "mip_solver"=>cbc, "time_limit"=>7200)
#scs = _FP.optimizer_with_attributes(SCS.Optimizer, "max_iters"=>100000)
#gurobi = Gurobi.Optimizer
#mosek = Mosek.Optimizer
optimizer = cbc

## Input parameters

number_of_hours = 24 # Number of time points
planning_horizon = 10 # Years to scale generation costs
file = "./test/data/case6/case6_replacement.m" # Input case, in Matpower m-file format: here 6-bus case with candidate AC, DC lines and candidate storage
scenario_properties = Dict(
    1 => Dict{String,Any}("probability"=>0.5, "start"=>1514764800000), # 1514764800000 is 2018-01-01T00:00, needed by `create_profile_data_italy!` when `mc=false`
    2 => Dict{String,Any}("probability"=>0.5, "start"=>1546300800000), # 1546300800000 is 2019-01-01T00:00, needed by `create_profile_data_italy!` when `mc=false`
)
scenario_metadata = Dict{String,Any}("mc"=>false) # Needed by `create_profile_data_italy!`

## Load test case

data = _FP.parse_file(file) # Parse input file to obtain data dictionary
_FP.add_dimension!(data, :hour, number_of_hours) # Add dimension, e.g. hours
_FP.add_dimension!(data, :scenario, scenario_properties; metadata=scenario_metadata) # Add dimension, e.g. scenarios
_FP.add_dimension!(data, :year, 1; metadata = Dict{String,Any}("scale_factor"=>planning_horizon)) # Add_dimension, e.g. years
_FP.scale_data!(data) # Scale investment & operational cost data based on planning years & hours
data, loadprofile, genprofile = create_profile_data_italy!(data) # Load time series data based demand and RES profiles of the six market zones in Italy from the data folder
time_series = create_profile_data(number_of_hours*_FP.dim_length(data,:scenario), data, loadprofile, genprofile) # Create time series data to be passed to the data dictionay
mn_data = _FP.make_multinetwork(data, time_series) # Create the multinetwork data dictionary

## Solve the planning problem

# PowerModels(ACDC) and FlexPlan settings
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => false, "allow_line_replacement" => true)

# Build optimisation model, solve it and write solution dictionary:
# This is the "problem file" which needs to be constructed individually depending on application
# In this case: multi-period optimisation of demand flexibility, AC & DC lines and storage investments
println("Solving planning problem...")
result = _FP.simple_stoch_flex_tnep(mn_data, _PM.DCPPowerModel, optimizer; setting = s)

println("Test completed")


text:
deck = Cint[4,4,4,4,4,4,4,4,4,16]
global d = 0

for i = 1:10

    # Dealer showing
    deck[i] -= 1
    p = 0
    for j = 1:10
        deck[j] -= 1
        n = ccall((:partitions, "./libpartitions.so"), Cint,
                  (Ptr{Cint}, Cint), deck, j)
        deck[j] += 1
        p += n
    end
    println("Dealer showing ", i-1," partitions = ",p)
    d += p
    deck[i] += 1
end

println("Total partitions = ",d)

max_stars_repo_name: octonion/recursion

max_issues_repo_name: octonion/recursion

max_forks_repo_name: octonion/recursion


text:
function get_table_config(opts)
    fields = [:ONodes,:CLevel,:Incumbent,:BestBound,:Gap,:Time]
    field_chars = [9,8,28,28,7,8]
    opbs = opts.branch_strategy

    if (opbs == :StrongPseudoCost || opbs == :Reliability) && opts.strong_restart
        push!(fields,:Restarts)
        push!(field_chars,10)
    end

    if opts.processors > 1
        pushfirst!(fields,:p)
        pushfirst!(field_chars, 3)
    end

    if opbs == :StrongPseudoCost || opbs == :PseudoCost || opbs == :Reliability
        push!(fields, :GainGap)
        push!(field_chars, 10)
    end
    return fields, field_chars
end

function print_table_header(fields, field_chars)
    println("")
    println(get_table_header_line(fields, field_chars))
    println(repeat("=", sum(field_chars)))
end

function get_table_header_line(fields, field_chars)
    ln = ""
    i = 1
    for f in fields
        f = string(f)
        padding = field_chars[i]-length(f)
        ln *= repeat(" ",trunc(Int, floor(padding/2)))
        ln *= f
        ln *= repeat(" ",trunc(Int, ceil(padding/2)))
        i += 1
    end
    return ln
end

function is_table_diff(fields, last_arr, new_arr)
    if length(last_arr) != length(new_arr)
        return true
    end

    time_idx = findfirst(fields .== :Time)
    if time_idx !== nothing
        last_arr = vcat(last_arr[1:time_idx-1],last_arr[time_idx+1:end])
        new_arr = vcat(new_arr[1:time_idx-1],new_arr[time_idx+1:end])
    end
    for i=1:length(last_arr)
        last_arr[i] != new_arr[i] && return true
    end
    return false
end

function print_table(p, tree, node, step_obj, start_time, fields, field_chars; last_arr=[])
    table_line, table_arr = get_table_line(p, tree, node, step_obj, start_time, fields, field_chars;
                                           last_arr=[])
    is_table_diff(fields, last_arr,table_arr) && println(table_line)
    return table_arr
end

function get_table_line(p, tree, node, step_obj, start_time, fields, field_chars; last_arr=[])
    gain_gap = step_obj.gain_gap
    opbs = tree.options.branch_strategy
    counter = step_obj.counter
    if (opbs != :StrongPseudoCost || counter > tree.options.strong_branching_nsteps) && step_obj.upd_gains != :GainsToTree
        step_obj.nrestarts = -1 # will be displayed as -
    end
    if counter <= tree.options.strong_branching_nsteps || step_obj.upd_gains == :GainsToTree
        gain_gap = -1.0 # will be displayed as -
    end

    arr = []
    nrestarts = step_obj.nrestarts
    i = 1
    ln = ""

    # check if incumbent and best bound should be printed and get the shorted field_chars
    found_both = 0
    min_field_chars = Inf
    j = 0
    for f in fields
        j += 1
        if f == :Incumbent || f ==:BestBound
            found_both += 1
            if field_chars[j] < min_field_chars
                min_field_chars = field_chars[j]
            end
        end
    end
    precision = 2
    if isdefined(tree,:incumbent) && found_both == 2
        if tree.incumbent.objval != tree.best_bound
            while round(tree.incumbent.objval; digits=precision) == round(tree.best_bound; digits=precision) && precision < 10
                if length(string(tree.incumbent.objval,precision+1)) > min_field_chars-2
                    break
                end
                if length(string(tree.best_bound,precision+1)) > min_field_chars-2
                    break
                end
                precision += 1
            end
        end
    end

    for f in fields
        val = ""
        if f == :p
            val = string(p)
        elseif f == :ONodes
            val = string(length(tree.branch_nodes))
        elseif f == :Incumbent
            if isdefined(tree,:incumbent)
                incu = tree.incumbent
                val = string(round(incu.objval; digits=precision))
            else
                val = "-"
            end
        elseif f == :BestBound
            val = string(round(tree.best_bound; digits=precision))
        elseif f == :Gap
            if isdefined(tree,:incumbent)
                b = tree.best_bound
                o = tree.incumbent.objval
                val = round(abs(b-o)/abs(o)*100; digits=2)
                if isnan(val)
                    if o == 0 && b == o
                        val = "0%"
                    else
                        val = "NaN%"
                    end
                else
                    if length(string(val)) > field_chars[i]
                        if val > 0 && val < tree.options.mip_gap*100
                            val = "< "*string(tree.options.mip_gap*100)*"%"
                        elseif val > 1000
                            val = ">>"
                        else
                            val = string(val)*"%"
                        end
                    else
                        val = string(val)*"%"
                    end
                end
            else
                val = "-"
            end
        elseif f == :Time
            val = string(round(time()-start_time; digits=1))
        elseif f == :Restarts
            if nrestarts == -1
                val = "-"
            else
                val = string(nrestarts)
            end
        elseif f == :CLevel
            val = string(node.level+1)
        elseif f == :GainGap
            if gain_gap == Inf
                val = "∞"
            elseif gain_gap == -1.0
                val = "-"
            else
                val = string(round(gain_gap*100; digits=1))*"%"
            end
            if length(val) > field_chars[i]
                val = ">>"
            end
        end
        if length(val) > field_chars[i]
            # too long to display shouldn't happen normally but is better than error
            # if it happens
            val = "t.l."
        end

        padding = field_chars[i]-length(val)
        ln *= repeat(" ",trunc(Int, floor(padding/2)))
        ln *= val
        ln *= repeat(" ",trunc(Int, ceil(padding/2)))
        push!(arr,val)
        i += 1
    end
    return ln, arr
end


text:
"""
Imagine we have several dataframes, each have columns "t" and "f(t)", and we try to get a long dataframe of columns "t", "f(t)" and "sheet" which is specified by `sheetnames`; in this case, `dfvstack(name_df_pairs, t_name, ft_name)` can help!

# Example
```julia
julia> name_df_pairs = sheetnames .=> dfs

Dict{Any, Any} with 3 entries:
  "s763.mat" => 4001×4 DataFrame…
  "s761.mat" => 1499×4 DataFrame…
  "s915.mat" => 824×4 DataFrame…

julia> t_name = "t";
julia> ft_name = "f(t)"
julia> longdf = dfvstack(name_df_pairs, t_name, ft_name)
```
"""
function dfvstack(name_df_pairs, t_name, ft_name; newcol=:sheet)

    dfs = []
    for (k, dfi0) in name_df_pairs
        dfi = copy(dfi0)
        select!(dfi, t_name, ft_name)
        insertcols!(dfi, newcol => k)
        push!(dfs, dfi)
    end
    return vcat(dfs...)
end


text:
using Pkg, Pkg.TOML

root = joinpath(@__DIR__, "..")

meta = TOML.parsefile(joinpath(@__DIR__, "Notebooks.toml"))
meta = meta[ARGS[1]]

path = meta["path"]
deps = get(meta, "deps", [])
deps = deps isa String ? [deps] : deps

for d in ["Project.toml", "Manifest.toml", ".gitignore"]
  isfile(joinpath(root, path, d)) && push!(deps, d)
end

mkpath(joinpath(root, "notebooks", path))
for dep in deps
  cp(joinpath(root, path, dep), joinpath(root, "notebooks", path, dep), force = true)
end

pushfirst!(LOAD_PATH, @__DIR__)
Pkg.activate(joinpath(root, "notebooks", path))

using Literate

function init_nb(content)
  content = "using Pkg; Pkg.activate(\"$root\"); Pkg.status();\n\n" * content
  return content
end

scripts = meta["notebook"]
scripts isa String && (scripts = [scripts])

for script in scripts
  Literate.notebook(joinpath(root, path, script),
                    joinpath(root, "notebooks", path),
                    credit = false, preprocess = init_nb)
end

scripts = map(x -> x[1:end - 3] * ".ipynb", scripts)
nbs = filter(x -> endswith(x, ".ipynb"), readdir(joinpath(root, path)))
keep = union(deps, scripts, nbs)
files = readdir(joinpath(root, "notebooks", path))

for r in files
  r in keep || rm(joinpath(root, "notebooks", path, r), force = true)
end

max_stars_repo_name: jstrube/model-zoo

max_issues_repo_name: jstrube/model-zoo

max_forks_repo_name: jstrube/model-zoo


text:
function arnoldi(a::AbstractArray{T,2},b::AbstractArray{T,1},m::S = length(b)) where {T <: AbstractFloat, S <: Integer} # Exact

    n = length(b)
    if m > n
        error("'m' must be no larger than the length of 'b' ")
    end

    q = zeros(n,m)
    h = zeros(m+1,m)

    β = norm(b)
    qk = copy(b)

    @views for k = 1:m
        q[:,k] .= qk/β
        qk .= a*q[:,k]
        for j = 1:k
            h[j,k] = qk'q[:,j]
            qk .-= h[j,k]*q[:,j]
        end
        β = norm(qk)
        h[k+1,k] = β
        if β < eps(T)
            return q[:,1:k], h[1:k,1:k], k
        end
    end

    return q, h[1:m,1:m], m

end

function arnoldi(f::Function,x::Array{T,1},r::Array{T,1},m::S=length(x)) where {T <: AbstractFloat, S <: Integer} #Exact

    n = length(x)
    if m > n
        error("'m' must be no larger than the length of 'b' ")
    end

    q = zeros(n,m)
    h = zeros(m+1,m)

    β = norm(r)
    qk = copy(r)

    @views for k = 1:m
        q[:,k] = qk/β
        jacvec!(qk,f,x,q[:,k])
        for j = 1:k
            h[j,k] = qk'q[:,j]
            qk -= h[j,k]*q[:,j]
        end
        β = norm(qk)
        h[k+1,k] = β
        if β < eps(T)
            return q[:,1:k], h[1:k,1:k], k
        end
    end

    return q, h[1:m,:], m

end

function arnoldi(a::AbstractArray{T,2},b::AbstractArray{T,1},tol::T,m::S=length(b)) where {T <: AbstractFloat, S <: Integer} #Inexact

    n = length(b)
    if m > n
        error("'m' must be no larger than the length of 'b' ")
    end

    q = zeros(n,m)
    h = zeros(m+1,m)

    β = norm(b)
    qk = copy(b)

    @views for k = 1:m
        q[:,k] .= qk/β
        qk .= a*q[:,k]
        for j = 1:k
            h[j,k] = qk'q[:,j]
            qk .-= h[j,k]*q[:,j]
        end
        β = norm(qk)
        h[k+1,k] = β
        if β < tol
            return q[1:n,1:k], h[1:k,1:k], k
        end
    end

    return q[1:n,1:m], h[1:m,1:m], m

end

function arnoldi(f::Function,x::Array{T,1},r::Array{T,1},tol::T,m::S=length(x)) where {T <: AbstractFloat, S <: Integer} # Inexact

    n = length(x)
    if m > n
        error("'m' must be no larger than the length of 'b' ")
    end

    q = zeros(n,m)
    h = zeros(m+1,m)

    β = norm(r)
    qk = copy(r)

    @views for k = 1:m
        q[:,k] = qk/β
        jacvec!(qk,f,x,q[:,k])
        for j = 1:k
            h[j,k] = qk'q[:,j]
            qk -= h[j,k]*q[:,j]
        end
        β = norm(qk)
        h[k+1,k] = β
        if β < tol
            return q[1:n,1:k], h[1:k,1:k], k
        end
    end

    return q[1:n,1:m], h[1:m,1:m], m

end

function arnoldi_inplace(f::Function,x::Array{T,1},r::Array{T,1},tol::T,m::S=length(x)) where {T <: AbstractFloat, S <: Integer} # Inexact

    n = length(x)
    if m > n
        error("'m' must be no larger than the length of 'b' ")
    end

    q = zeros(n,m)
    h = zeros(m+1,m)

    β = norm(r)
    qk = copy(r)

    @views for k = 1:m
        q[:,k] = qk/β
        jacvec_inplace!(qk,f,x,q[:,k])
        for j = 1:k
            h[j,k] = qk'q[:,j]
            qk -= h[j,k]*q[:,j]
        end
        β = norm(qk)
        h[k+1,k] = β
        if β < tol
            return q[1:n,1:k], h[1:k,1:k], k
        end
    end

    return q[1:n,1:m], h[1:m,1:m], m

end


text:
# Mu' = f
struct SteadyStateProblem{uType,isinplace,P,F} <: AbstractSteadyStateProblem{uType,isinplace}
  f::F
  u0::uType
  p::P
  @add_kwonly function SteadyStateProblem{iip}(f::AbstractODEFunction{iip},
                                               u0,p=nothing) where {iip}
    new{typeof(u0),isinplace(f),typeof(p),typeof(f)}(f,u0,p)
  end

  function SteadyStateProblem{iip}(f,u0,p=nothing) where iip
    SteadyStateProblem(ODEFunction{iip}(f),u0,p)
  end
end

function SteadyStateProblem(f::AbstractODEFunction,u0,p=nothing;kwargs...)
  SteadyStateProblem{isinplace(f)}(f,u0,p;kwargs...)
end

function SteadyStateProblem(f,u0,p=nothing;kwargs...)
  SteadyStateProblem(ODEFunction(f),u0,p;kwargs...)
end

SteadyStateProblem(prob::AbstractODEProblem) =
      SteadyStateProblem{isinplace(prob)}(prob.f,prob.u0,prob.p)


text:
"""
    angle_of_incidence(lat,lon,jd,time)

Calculate the cosine of the sun incident angle for a given location (lat,lon),
julian day (jd) and time in hours and fraction of hours.
"""
function angle_of_incidence(lat, lon, jd, time)
    # inclination angle
    delta = (asin(-sin(23.45 / 180 * pi) .* cos(360 / 365.25 * (jd + 10) / 180 * pi)))
    # lat to radians
    lamda = lat / 180 * pi
    # time to radians
    omega = time / 24 * 2 * pi + lon / 360 * 2 * pi + pi

    return cos(delta) * cos(lamda) * cos(omega) + sin(delta) * sin(lamda)
end

"""
    delta(jd)

Calculates the Sun's declination angle as a function of Julian day (jd)
"""
delta(n) = (asin(-sin(23.45 / 180 * pi) .* cos(360 / 365.25 * (n + 10) / 180 * pi)))

"""
    heaviside(x)

Heaviside function : 0 when x<0 and 1 when x>=0
"""
function heaviside(x)
    y = zeros(size(x))
    y[x.>=0] .= 1
    return y
end


text:
# ## Type Structure

# The storage of power system data has traditionally been done with general use tables based on the power flow problem. However, tables are not inherently designed to store data with mixed data representations or hierarchal structures. This limitation was addressed in the early 1990's with the advent of automation, spurred by increasingly complex data needs for power systems operations. The industry required standardized models to exchange larger and more complex sets of information, and so resorted to an object-oriented data model. The CIM was developed and later made a standard maintained by the IEC - Technical Committee 57 Working Group 13. The aim was to provide a standard definition for power system components that could be used in automated EMS and asset-management databases.

# From the computer engineering point-of-view, using tables would require to re-write a lot of code. In most languages with strongly type hierarchies, using a tree structure it is advantageous. This is particularly relevant for code re-use and data encapsulation. Collecting all potential devices is not possible; neither is it desirable. However, it is necessary to have the capability to extend the data model such that different user can store their data in an organized ontology.

# Each device in PowerSystems.jl is defined using a Julia structure embedded in a type hierarchy. This implementation enables the categorization of the devices by their operational characteristics and the data required to model them. The type-based and method-forward paradigm in PowerSystems.jl incentivizes code re-use in modeling, enhances the reproducibility of code, and enables modularization for specific data types.


text:
function get_norm_rc!(v, A_colptr, A_rowval, A_nzval, n, ax)
  T = eltype(v)
  v .= zero(T)
  for j = 1:n
    @inbounds @simd for i = A_colptr[j]:(A_colptr[j + 1] - 1)
      k = ax == :row ? A_rowval[i] : j
      if abs(A_nzval[i]) > v[k]
        v[k] = abs(A_nzval[i])
      end
    end
  end

  v .= sqrt.(v)
  @inbounds @simd for i = 1:length(v)
    if v[i] == zero(T)
      v[i] = one(T)
    end
  end
end

function mul_A_D1_D2!(A_colptr, A_rowval, A_nzval, d1, d2, r, c, uplo)
  for j = 1:length(c)
    @inbounds @simd for i = A_colptr[j]:(A_colptr[j + 1] - 1)
      A_nzval[i] /= r[A_rowval[i]] * c[j]
    end
  end
  if uplo == :U
    d1 ./= c
    d2 ./= r
  else
    d1 ./= r
    d2 ./= c
  end
end

function mul_A_D3!(A_colptr, A_rowval, A_nzval, n, d3, uplo)
  for j = 1:n
    @inbounds @simd for i = A_colptr[j]:(A_colptr[j + 1] - 1)
      A_nzval[i] *= (uplo == :U) ? d3[A_rowval[i]] : d3[j]
    end
  end
end

function mul_Q_D!(Q_colptr, Q_rowval, Q_nzval, d, c)
  for j = 1:length(d)
    @inbounds @simd for i = Q_colptr[j]:(Q_colptr[j + 1] - 1)
      Q_nzval[i] /= c[Q_rowval[i]] * c[j]
    end
  end
  d ./= c
end

function mul_Q_D2!(Q_colptr, Q_rowval, Q_nzval, d2)
  for j = 1:length(d2)
    @inbounds @simd for i = Q_colptr[j]:(Q_colptr[j + 1] - 1)
      Q_nzval[i] *= d2[Q_rowval[i]] * d2[j]
    end
  end
end

function scaling_Ruiz!(
  fd_T0::QM_FloatData{T},
  id::QM_IntData,
  sd::ScaleData{T},
  ϵ::T;
  max_iter::Int = 100,
) where {T <: Real}
  d1, d2, d3, r_k, c_k = sd.d1, sd.d2, sd.d3, sd.r_k, sd.c_k
  # scaling Q (symmetric)
  if length(fd_T0.Q.data.rowval) > 0
    if fd_T0.uplo == :U
      get_norm_rc!(r_k, fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, id.nvar, :row)
      convergence = maximum(abs.(one(T) .- r_k)) <= ϵ
      mul_Q_D!(fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, d3, r_k)
    else
      get_norm_rc!(c_k, fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, id.nvar, :col)
      convergence = maximum(abs.(one(T) .- c_k)) <= ϵ
      mul_Q_D!(fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, d3, c_k)
    end
    k = 1
    while !convergence && k < max_iter
      if fd_T0.uplo == :U
        get_norm_rc!(
          r_k,
          fd_T0.Q.data.colptr,
          fd_T0.Q.data.rowval,
          fd_T0.Q.data.nzval,
          id.nvar,
          :row,
        )
        convergence = maximum(abs.(one(T) .- r_k)) <= ϵ
        mul_Q_D!(fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, d3, r_k)
      else
        get_norm_rc!(
          c_k,
          fd_T0.Q.data.colptr,
          fd_T0.Q.data.rowval,
          fd_T0.Q.data.nzval,
          id.nvar,
          :col,
        )
        convergence = maximum(abs.(one(T) .- c_k)) <= ϵ
        mul_Q_D!(fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, d3, c_k)
      end
      k += 1
    end
    mul_A_D3!(fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, fd_T0.A.n, d3, fd_T0.uplo)
    fd_T0.c .*= d3
    fd_T0.lvar ./= d3
    fd_T0.uvar ./= d3
  end

  # r (resp. c) norm of rows of AT (resp. cols)
  # scaling: D2 * AT * D1
  get_norm_rc!(r_k, fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, fd_T0.A.n, :row)
  get_norm_rc!(c_k, fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, fd_T0.A.n, :col)
  convergence = maximum(abs.(one(T) .- r_k)) <= ϵ && maximum(abs.(one(T) .- c_k)) <= ϵ
  mul_A_D1_D2!(fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, d1, d2, r_k, c_k, fd_T0.uplo)
  k = 1
  while !convergence && k < max_iter
    get_norm_rc!(r_k, fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, fd_T0.A.n, :row)
    get_norm_rc!(c_k, fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, fd_T0.A.n, :col)
    convergence = maximum(abs.(one(T) .- r_k)) <= ϵ && maximum(abs.(one(T) .- c_k)) <= ϵ
    mul_A_D1_D2!(fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, d1, d2, r_k, c_k, fd_T0.uplo)
    k += 1
  end
  length(fd_T0.Q.data.rowval) > 0 &&
    mul_Q_D2!(fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, d2)
  fd_T0.b .*= d1
  fd_T0.c .*= d2
  fd_T0.lvar ./= d2
  fd_T0.uvar ./= d2
end

function div_D2D3_Q_D3D2!(Q_colptr, Q_rowval, Q_nzval, d2, d3, n)
  for j = 1:n
    @inbounds @simd for i = Q_colptr[j]:(Q_colptr[j + 1] - 1)
      Q_nzval[i] /= d2[Q_rowval[i]] * d2[j] * d3[Q_rowval[i]] * d3[j]
    end
  end
end

function div_D1_A_D2D3!(A_colptr, A_rowval, A_nzval, d1, d2, d3, n, uplo)
  for j = 1:n
    @inbounds @simd for i = A_colptr[j]:(A_colptr[j + 1] - 1)
      A_nzval[i] /=
        (uplo == :U) ? d1[j] * d2[A_rowval[i]] * d3[A_rowval[i]] : d1[A_rowval[i]] * d2[j] * d3[j]
    end
  end
end

function post_scale!(
  d1::Vector{T},
  d2::Vector{T},
  d3::Vector{T},
  pt::Point{T},
  res::AbstractResiduals{T},
  fd_T0::QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
) where {T <: Real}
  pt.x .*= d2 .* d3
  div_D2D3_Q_D3D2!(fd_T0.Q.data.colptr, fd_T0.Q.data.rowval, fd_T0.Q.data.nzval, d2, d3, id.nvar)
  mul!(itd.Qx, fd_T0.Q, pt.x)
  itd.xTQx_2 = dot(pt.x, itd.Qx) / 2
  div_D1_A_D2D3!(fd_T0.A.colptr, fd_T0.A.rowval, fd_T0.A.nzval, d1, d2, d3, fd_T0.A.n, fd_T0.uplo)
  pt.y .*= d1
  fd_T0.uplo == :U ? mul!(itd.ATy, fd_T0.A, pt.y) : mul!(itd.ATy, fd_T0.A', pt.y)
  fd_T0.uplo == :U ? mul!(itd.Ax, fd_T0.A', pt.x) : mul!(itd.Ax, fd_T0.A, pt.x)
  fd_T0.b ./= d1
  fd_T0.c ./= d2 .* d3
  itd.cTx = dot(fd_T0.c, pt.x)
  itd.pri_obj = itd.xTQx_2 + itd.cTx + fd_T0.c0
  fd_T0.lvar .*= d2 .* d3
  fd_T0.uvar .*= d2 .* d3
  pt.s_l ./= @views d2[id.ilow] .* d3[id.ilow]
  pt.s_u ./= @views d2[id.iupp] .* d3[id.iupp]
  itd.dual_obj =
    dot(fd_T0.b, pt.y) - itd.xTQx_2 + dot(pt.s_l, view(fd_T0.lvar, id.ilow)) -
    dot(pt.s_u, view(fd_T0.uvar, id.iupp)) + fd_T0.c0
  res.rb .= itd.Ax .- fd_T0.b
  res.rc .= itd.ATy .- itd.Qx .- fd_T0.c
  res.rc[id.ilow] .+= pt.s_l
  res.rc[id.iupp] .-= pt.s_u
  #         rcNorm, rbNorm = norm(rc), norm(rb)
  res.rcNorm, res.rbNorm = norm(res.rc, Inf), norm(res.rb, Inf)
end


text:
using Documenter, IterableTables

makedocs(
	modules = [IterableTables],
	analytics="UA-132838790-1",
	sitename = "IterableTables.jl",
	pages = [
		"Introduction" => "index.md",
        "User Guide" => "userguide.md",
        "Integration Guide" => "integrationguide.md",
        "Developer Guide" => "developerguide.md"]
)

deploydocs(
    repo = "github.com/queryverse/IterableTables.jl.git"
)


text:
using Criterion

function test1_svd(k::Integer)
    for _ in k
        svd(zeros(20, 20))
    end
end

function test2_svd(k::Integer)
    for _ in k
        svd(ones(20, 20))
    end
end

function test3_svd(k::Integer)
    for _ in k
        svd(ones(10, 10))
    end
end

env = Criterion.measure_environment()

bench1 = Criterion.Benchmark("test1_svd", 1, test1_svd)
times, sample_analysis, outliers = Criterion.run_analyze_one(env, "test1_svd", bench1)
r1 = Criterion.Report(1, "test1/svd", times, sample_analysis, outliers)

bench2 = Criterion.Benchmark("test2_svd", 2, test2_svd)
times, sample_analysis, outliers = Criterion.run_analyze_one(env, "test2_svd", bench2)
r2 = Criterion.Report(2, "test2/svd", times, sample_analysis, outliers)

bench3 = Criterion.Benchmark("test3_svd", 2, test3_svd)
times, sample_analysis, outliers = Criterion.run_analyze_one(env, "test3_svd", bench3)
r3 = Criterion.Report(3, "test3/svd", times, sample_analysis, outliers)

Criterion.report([r1, r2, r3])

#using Winston
#p = FramedPlot()
#setattr(p, "title", "Benchmark Times (sec)")
#setattr(p, "ylabel", "counts")
#setattr(p, "xlabel", "times (sec)")
#add(p, Histogram(hist(times, 100)...))
#add(p, LineX(sa.mean.point, "color", "red"))
#add(p, LineX(sa.mean.lbound, "type", "dot"))
#add(p, LineX(sa.mean.ubound, "type", "dot"))
#Winston.display(p)


text:
#using Base.Threads
import StatsBase: BaseRadixSortSafeTypes
import Base: Forward, ForwardOrdering, Reverse, ReverseOrdering, sortperm, Ordering,
            setindex!, getindex, similar, sort!, Reverse

import Base.Threads: nthreads, threadid

# create bits types for easy loading of bytes of lengths up to 15
primitive type Bits24 24 end
primitive type Bits40 40 end
primitive type Bits48 48 end
primitive type Bits56 56 end
primitive type Bits72 72 end
primitive type Bits80 80 end
primitive type Bits88 88 end
primitive type Bits96 96 end
primitive type Bits104 104 end
primitive type Bits112 112 end
primitive type Bits120 120 end

# winner from benchmarks/which_is_the_fastest_UInt_histogram.jl
"""
    uint_hist(bits, [RADIX_SIZE = 16, RADIX_MASK = 0xffff])
Computes a histogram (counts) for the vector RADIX_SIZE bits at a time. E.g. if eltype(bits) is UInt64 and RADIX_SIZE is 16
then 4 histograms are created for each of the 16 bit chunks.
"""
function uint_hist(bits::Vector{T}, RADIX_SIZE = 16, RADIX_MASK = 0xffff) where T
    #println("hello")
    iter = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    #hist = zeros(UInt32, 2^RADIX_SIZE, iter, nthreads())
    hist = zeros(UInt32, 2^RADIX_SIZE, iter)
    #@threads for j = 1:length(bits)
    for j = 1:length(bits)
        for i = 0:iter-1
            @inbounds hist[1+Int((bits[j] >> (i * RADIX_SIZE)) & RADIX_MASK), i+1] += 1
        end
    end
    #@threads for j in 1:iter
    # for j in 1:iter
    #     for i = 2:nthreads()
    #         @inbounds hist[:, j, 1] .+= hist[:, j, i]
    #     end
    # end
    hist
end

# sort it by using sorttwo! on the pointer
function sort_spointer!(svec::Vector{String}, lo::Int, hi::Int, o::O) where O <: Union{ForwardOrdering, ReverseOrdering}
    if lo >= hi;  return svec;  end
    # the length subarray to sort
    l = hi - lo + 1

    # find the maximum string length
    lens = maximum(sizeof, svec)
    skipbytes = lens

    if lens > 0
        while lens > 4
            skipbytes = max(0, skipbytes - 8)
            bits64 = zeros(UInt64, l)
            if o == Reverse
                bits64[lo:hi] .= .~load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits64, svec, lo, hi)
            else
                bits64[lo:hi] .= load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits64, svec, lo, hi)
            end
            lens -= 8
        end
        if lens > 0
            skipbytes = max(0, skipbytes - 4)
            bits32 = zeros(UInt32, l)
            if o == Reverse
                bits32[lo:hi] .= .~load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits32, svec, lo, hi)
            else
                bits32[lo:hi] .= load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                SortingLab.sorttwo!(bits32, svec, lo, hi)
            end
            lens -= 4
        end
    end
    # unsafe_pointer_to_objref.(ptrs - 8)
    svec
end

"""
    sorttwo!(vs, index)

Sort both the `vs` and reorder `index` at the same. This allows for faster sortperm
for radix sort.
"""
function sorttwo!(vs::Vector{T}, index, lo::Int = 1, hi::Int=length(vs), RADIX_SIZE = 16, RADIX_MASK = 0xffff) where T <:Union{BaseRadixSortSafeTypes}
    # Input checking
    if lo >= hi;  return (vs, index);  end
    #println(vs)
    o = Forward

    # Init
    iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
    # number of buckets in the counting step
    nbuckets = 2^RADIX_SIZE

    # Histogram for each element, radix
    bin = uint_hist(vs, RADIX_SIZE, RADIX_MASK)

    # bin = zeros(UInt32, nbuckets, iters)
    # if lo > 1;  bin[1,:] = lo-1;  end

    # Sort!
    swaps = 0
    len = hi-lo+1

    index1 = similar(index)
    ts=similar(vs)
    for j = 1:iters
        # Unroll first data iteration, check for degenerate case
        v = uint_mapping(o, vs[hi])
        idx = Int((v >> ((j-1)*RADIX_SIZE)) & RADIX_MASK) + 1

        # are all values the same at this radix?
        if bin[idx,j] == len;  continue;  end

        # cbin = cumsum(bin[:,j])
        # tries to achieve the above one-liner with more efficiency
        cbin = zeros(UInt32, nbuckets)
        cbin[1] = bin[1,j]
        for i in 2:nbuckets
            cbin[i] = cbin[i-1] + bin[i,j]
        end

        ci = cbin[idx]
        #println((ci, hi))
        ts[ci] = vs[hi]
        index1[ci] = index[hi]
        # println(cbin[idx])
        cbin[idx] -= 1

        # Finish the loop...
        #@inbounds
        for i in hi-1:-1:lo
            v = uint_mapping(o, vs[i])
            idx = Int((v >> ((j-1)*RADIX_SIZE)) & RADIX_MASK) + 1
            ci = cbin[idx]
            #println(ci)
            ts[ci] = vs[i]
            index1[ci] = index[i]
            cbin[idx] -= 1
        end
        vs,ts = ts,vs
        index, index1 = index1, index
        swaps += 1
    end

    if isodd(swaps)
        vs,ts = ts,vs
        index, index1 = index1, index
        for i = lo:hi
            @inbounds vs[i] = ts[i]
            @inbounds index[i] = index1[i]
        end
    end
    (vs, index)
end

# loads `remaining_bytes_to_load` bytes from `ptrs` which is a C-style pointer to a string
# these functions assumes that remaining_bytes_to_load > 0
function load_bits_with_padding(::Type{UInt128}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt128
    nbits_to_shift_away = 8(sizeof(UInt128) - remaining_bytes_to_load)

    # the below checks if the string is less than 16 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xff0
        if  remaining_bytes_to_load == 15
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits120}(ptrs))))
        elseif  remaining_bytes_to_load == 14
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits112}(ptrs))))
        elseif  remaining_bytes_to_load == 13
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits104}(ptrs))))
        elseif  remaining_bytes_to_load == 12
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits96}(ptrs))))
        elseif  remaining_bytes_to_load == 11
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits88}(ptrs))))
        elseif  remaining_bytes_to_load == 10
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits80}(ptrs))))
        elseif  remaining_bytes_to_load == 9
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits72}(ptrs))))
        elseif  remaining_bytes_to_load == 8
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt64}(ptrs))))
        elseif  remaining_bytes_to_load == 7
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits56}(ptrs))))
        elseif  remaining_bytes_to_load == 6
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits48}(ptrs))))
        elseif  remaining_bytes_to_load == 5
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits40}(ptrs))))
        elseif  remaining_bytes_to_load == 4
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt32}(ptrs))))
        elseif  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt128, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt128}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

function load_bits_with_padding(::Type{UInt64}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt64
    nbits_to_shift_away = 8(sizeof(UInt64) - remaining_bytes_to_load)

    # the below checks if the string is less than 8 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xff8
        if  remaining_bytes_to_load == 7
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits56}(ptrs))))
        elseif  remaining_bytes_to_load == 6
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits48}(ptrs))))
        elseif  remaining_bytes_to_load == 5
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits40}(ptrs))))
        elseif  remaining_bytes_to_load == 4
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt32}(ptrs))))
        elseif  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt64, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt64}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

function load_bits_with_padding(::Type{UInt32}, ptrs::Ptr{UInt8}, remaining_bytes_to_load)::UInt32
    nbits_to_shift_away = 8(sizeof(UInt32) - remaining_bytes_to_load)

    # the below checks if the string is less than 4 bytes away from the page
    # boundary assuming the page size is 4kb
    # see https://discourse.julialang.org/t/is-there-a-way-to-check-how-far-away-a-pointer-is-from-a-page-boundary/8147/11?u=xiaodai
    if (UInt(ptrs) & 0xfff) > 0xffc
        if  remaining_bytes_to_load == 3
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{Bits24}(ptrs))))
        elseif  remaining_bytes_to_load == 2
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{UInt16}(ptrs))))
        else
            return ntoh(Base.zext_int(UInt32, unsafe_load(Ptr{UInt8}(ptrs))))
        end
    else
        return ntoh(unsafe_load(Ptr{UInt32}(ptrs))) >> nbits_to_shift_away << nbits_to_shift_away
    end
end

"""
    load_bits([type,] s, skipbytes)

Load the underlying bits of a string `s` into a `type` of the user's choosing.
The default is `UInt`, so on a 64 bit machine it loads 64 bits (8 bytes) at a time.
If the `String` is shorter than 8 bytes then it's padded with 0.

- `type`:       any bits type that has `>>`, `<<`, and `&` operations defined
- `s`:          a `String`
- `skipbytes`:  how many bytes to skip e.g. load_bits("abc", 1) will load "bc" as bits
"""
# Some part of the return result should be padded with 0s.
# To prevent any possibility of segfault we load the bits using
# successively smaller types
# it is assumed that the type you are trying to load into needs padding
# i.e. `remaining_bytes_to_load > 0`
function load_bits(::Type{T}, s::String, skipbytes = 0)::T where T <: Unsigned
    n = sizeof(s)
    if n < skipbytes
        res = zero(T)
    elseif n - skipbytes >= sizeof(T)
        res = ntoh(unsafe_load(Ptr{T}(pointer(s, skipbytes+1))))
    else
        ptrs  = pointer(s) + skipbytes
        remaining_bytes_to_load = n - skipbytes
        res = load_bits_with_padding(T, ptrs, remaining_bytes_to_load)
    end
    return res
end

# Radix sort for strings
"""
    radixsort!(svec)

Applies radix sort to the string vector, svec, and sort it in place.
"""
function radixsort!(svec::AbstractVector{String}, lo::Int, hi::Int, o::O; RADIX_SIZE = 16, RADIX_MASK = 0xffff) where O <: Union{ForwardOrdering, ReverseOrdering}
    if lo >= hi;  return svec;  end
    # the length subarray to sort
    l = hi - lo + 1

    # find the maximum string length
    lens = maximum(sizeof, svec)
    skipbytes = lens
    if lens > 0
        # while lens > 8
        #     skipbytes = max(0, skipbytes - 16)
        #     bits128 = zeros(UInt128, l)
        #     if o == Reverse
        #         bits128[lo:hi] .= .~load_bits.(UInt128, @view(svec[lo:hi]), skipbytes)
        #         sorttwo!(bits128, svec, lo, hi)
        #     else
        #         bits128[lo:hi] .= load_bits.(UInt128, @view(svec[lo:hi]), skipbytes)
        #         sorttwo!(bits128, svec, lo, hi)
        #     end
        #     lens -= 16
        # end
        while lens > 4
            skipbytes = max(0, skipbytes - 8)
            bits64 = zeros(UInt64, l)
            if o == Reverse
                bits64[lo:hi] .= .~load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits64, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            else
                bits64[lo:hi] .= load_bits.(UInt64, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits64, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            end
            lens -= 8
        end
        if lens > 0
            skipbytes = max(0, skipbytes - 4)
            bits32 = zeros(UInt32, l)
            if o == Reverse
                bits32[lo:hi] .= .~load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits32, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            else
                bits32[lo:hi] .= load_bits.(UInt32, @view(svec[lo:hi]), skipbytes)
                sorttwo!(bits32, svec, lo, hi, RADIX_SIZE, RADIX_MASK)
            end
            lens -= 4
        end
    end
    svec
end

"""
    radixsort(svec, rev = false)

Applies radix sort to the string vector, svec, and sort it in place.
"""
fsort(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) = radixsort(svec, rev , (16, 0xffff))
fsort!(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) = radixsort!(svec, rev , (16, 0xffff))
radixsort(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) = radixsort!(copy(svec), rev, radix_opts)
radixsort!(svec::Vector{String}, rev = false, radix_opts = (16, 0xffff)) =  radixsort!(svec, 1, length(svec), rev ? Base.Reverse : Base.Forward, RADIX_SIZE = radix_opts[1], RADIX_MASK = radix_opts[2])


text:
"""
Internal profile checker which check if the config file exist or not.

# Argument
- `path` : File path of the profiles
"""
function profiles_exist(path::String) :: Bool
    if isfile(path)
        return true
    end
    return false
end

"""
Split profile strings in to property to add in to dict collection.

# Argument
- property string with signature key=pair

# Example
```
split_string("repo=my_beautiful_repo")
```
"""
function split_string(s::String)::Pair
    property = split(s, "=")
    return string(property[1]) => string(property[2])
end

"""
Invalidate the loaded properties from profiles config file.

# Argument
- `p` profiles Dict will all profiles from config file
"""
function validate_profiles(p::Dict)
    all_keys = [
        "project_root", "projects_dir", "data_dir",
        "insights_dir", "insights_viz_dir",
        "insights_data_dir",
    ]
    missing_keys = []
    for (key, value) in p
        for property in all_keys
            if !haskey(value, property)
                push!(missing_keys, property)
            end
        end
    end
    if !isempty(missing_keys)
        msg = join(missing_keys, ", ")
        err = ErrorException("$msg keys does not exist in profile.")
        throw(err)
    end
end

"""
Load the profiles from the users profile setup file.

# Argument
- Path of profile file.

# Example

```
const ROOT_DIR = joinpath(homedir(), ".projectflow/profiles")
profiles = load_profiles(ROOT_DIR)
```
"""
function load_profiles(path::String)::Dict
    path = joinpath(path, "profiles")
    all_profiles = Dict{String, Dict{String, String}}()
    if profiles_exist(path)
        is_header = false
        global header = ""
        for line in eachline(path)
            global header
            if occursin(r"\[\w+\]", line)
                is_header = true
                header = replace(line, r"[\[\]]" => "")
                all_profiles[header] = Dict{String, String}()
            else
                is_header = false
            end
            if !is_header && !isempty(strip(line))
                property = split_string(line)
                push!(all_profiles[header], property)
            end
        end
    else
        err = ErrorException("$path file does not exist")
        throw(err)
    end
    try
        validate_profiles(all_profiles)
    catch ex
        rethrow(ex)
    end
    return all_profiles
end

"""
function select_property(path::String, name::String)

Select property of a profile from a given profile value.

# Argument
- path: Path for Profile
- name: Name of the Profile

```
Example:
const ROOT_DIR = joinpath(homedir(), ".projectflow/")
select_property(ROOT_DIR, "default")
```
"""
function select_property(path::String, name::String)
    try
        profiles = load_profiles(path)
        property = profiles[name]
        root = property["project_root"]
        property["projects_dir"] = joinpath(root, property["projects_dir"])
        property["data_dir"] = joinpath(root, property["data_dir"])
        property["insights_dir"] = joinpath(root, property["insights_dir"])
        return property
    catch ex
        rethrow(ex)
    end
end


text:
@testset "Base.show methods" begin
    txt = "sample text"
    buf = IOBuffer()
    # Test Document types
    for dt in [StringDocument, NGramDocument, TokenDocument]
        doc = dt(txt)
        Test.@test try show(buf, doc); true
                   catch; false end
        Test.@test try print(buf, summary(doc)); true
                   catch; false end
    end
    # Metadata
    doc = StringDocument(txt)
    Test.@test try show(buf, doc.metadata); true
               catch; false end
    # Corpus
    crps = Corpus([doc])
    Test.@test try show(buf, crps); true
               catch; false end
    Test.@test try print(buf, summary(crps)); true
               catch; false end
    #Document Term Matrix
    update_lexicon!(crps)
    update_inverse_index!(crps)
    m = DocumentTermMatrix(crps);
    Test.@test try show(buf,m); true
               catch; false end
    c = CooMatrix(crps);
    Test.@test try show(buf,c); true
               catch; false end
    # Stemmer
    s = Stemmer("english")
    Test.@test try show(buf,s); true
               catch; false end

end


text:
s = """
ex = :(1+x)
f = eval(:(x->(\$ex)))
"""
msgs = lintstr(s)
@test isempty(msgs)


text:
module RulesNormalMeanPrecisionPrecisionTest

using Test
using ReactiveMP
using Random

import ReactiveMP: @test_rules

@testset "rules:NormalMeanPrecision:precision" begin

    @testset "Variational: (q_out::Any, q_μ::Any)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out = PointMass(-1.0), q_μ = PointMass(2.0)), output = Gamma(1.5, 2.0 / 9.0)),
            (input = (q_out = PointMass(1.0), q_μ = PointMass(2.0)),  output = Gamma(1.5, 2.0)),
            (input = (q_out = PointMass(2.0), q_μ = PointMass(1.0)),  output = Gamma(1.5, 2.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out = NormalMeanVariance(-1.0, 2.0), q_μ = PointMass(2.0)),         output = Gamma(1.5, 2.0 / 11.0)),
            (input = (q_out = NormalMeanPrecision(1.0, 4.0), q_μ = PointMass(3.0)),         output = Gamma(1.5, 2.0 / 4.25)),
            (input = (q_out = NormalWeightedMeanPrecision(2.0, 4.0), q_μ = PointMass(1.0)), output = Gamma(1.5, 4.0))
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out = PointMass(2.0), q_μ = NormalMeanVariance(-1.0, 2.0)),         output = Gamma(1.5, 2.0 / 11.0)),
            (input = (q_out = PointMass(3.0), q_μ = NormalMeanPrecision(1.0, 4.0)),         output = Gamma(1.5, 2.0 / 4.25)),
            (input = (q_out = PointMass(1.0), q_μ = NormalWeightedMeanPrecision(2.0, 4.0)), output = Gamma(1.5, 4.0))
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out = NormalMeanVariance(2.0, 4.0), q_μ = NormalMeanVariance(-1.0, 2.0)),                  output = Gamma(1.5, 2.0 / 15.0)),
            (input = (q_out = NormalMeanPrecision(3.0, 3.0), q_μ = NormalMeanPrecision(1.0, 4.0)),                 output = Gamma(1.5, 24.0 / 55.0)),
            (input = (q_out = NormalWeightedMeanPrecision(1.0, 1.0), q_μ = NormalWeightedMeanPrecision(2.0, 4.0)), output = Gamma(1.5, 2.0 / 1.5))
        ]

    end

    @testset "Variational: (q_out_μ::Any)" begin

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out_μ = MvNormalMeanCovariance([ 1.0, 1.0 ], [ 1.0 0.0; 0.0 1.0 ]), ),   output = Gamma(1.5, 1.0)),
            (input = (q_out_μ = MvNormalMeanCovariance([ 2.0, 3.0 ], [ 2.0 -0.1; -0.1 3.0 ]), ), output = Gamma(1.5, 2.0 / 6.2)),
            (input = (q_out_μ = MvNormalMeanCovariance([ 4.0, 1.0 ], [ 4.0 1.0; 1.0 9.0 ]), ),   output = Gamma(1.5, 2.0 / 20.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out_μ = MvNormalMeanPrecision([ 1.0, 1.0 ], [ 1.0 0.0; 0.0 1.0 ]), ),   output = Gamma(1.5, 1.0)),
            (input = (q_out_μ = MvNormalMeanPrecision([ 2.0, 3.0 ], [ 2.0 -0.1; -0.1 3.0 ]), ), output = Gamma(1.5, 1198.0 / 1079.0)),
            (input = (q_out_μ = MvNormalMeanPrecision([ 4.0, 1.0 ], [ 4.0 1.0; 1.0 9.0 ]), ),   output = Gamma(1.5, 70.0 / 330.0)),
        ]

        @test_rules [ with_float_conversions = true ] NormalMeanPrecision(:τ, Marginalisation) [
            (input = (q_out_μ = MvNormalWeightedMeanPrecision([ 1.0, 1.0 ], [ 1.0 0.0; 0.0 1.0 ]), ),  output = Gamma(1.5, 1.0)),
            (input = (q_out_μ = MvNormalWeightedMeanPrecision([ 4.0, 1.0 ], [ 4.0 1.0; 1.0 9.0 ]), ),  output = Gamma(1.5, 1.4)),
        ]

    end

end

end


text:
@testset "Truth Tables" begin
    a, b, c = LogicalSymbol.([:a, :b, :c])

    # a ∧ b  (a and b)
    st = a ∧ b
    and_tt = truthtable(st)
    @test all(and_tt[!, :a] .== [false, true, false, true])
    @test all(and_tt[!, :b] .== [false, false, true, true])
    @test all(and_tt[!, string(st)] .== [false, false, false, true])

    # b ∨ a  (b or a)
    st = b ∨ a
    or_tt = truthtable(st)
    @test all(or_tt[!, :b] .== [false, false, true, true])
    @test all(or_tt[!, :a] .== [false, true, false, true])
    @test all(or_tt[!, string(st)] .== [false, true, true, true])

    # material implication
    st = a → b
    mi_tt = truthtable(st)
    @test all(mi_tt[!, string(st)] .== [true, false, true, true])

    # slightly more complex statement
    st = (a → b) ∧ c
    mi_tt = truthtable(st)
    @test all(mi_tt[!, string(st)] .== [false, false, false, false, true, false, true, true])
end