text:
using Test
using Printf
using FileIO

# ScatterPatch depends on this
include("../src/image.jl")
include("../src/synthResult.jl")

# the tested code
include("../src/scatterPatch.jl")

include("tinyImage.jl")

@testset "ScatterPatch" begin

    # This mimics the startup of Resynthesizer
    testImage = createTestMaskedImage()
    testSynthResult = initialSynthResult(testImage)
    testSortedOffsets = createSortedOffsets(testImage.image)

    @debug "Sorted offset" testSortedOffsets

    # scatter patch distinguished point (origin) is upper left of synth
    testDistinguishedPoint = PointInMaskedImage( CartesianIndex(3,2), testImage)

    testScatterPatch = prepareScatterPatch(testDistinguishedPoint, testSynthResult, testSortedOffsets)

    # Original C algorithm.
    # The patchPoint origin is in the ScatterPatch (at offset 0,0)
    # We also know it should be the first element, since ordered by offset
    # @test any(neighbor.offset == CartesianIndex(0,0) for neighbor in testScatterPatch))
    #@test testScatterPatch[1].offset == CartesianIndex(0,0)

    # The patchPoint origin *not* in the ScatterPatch (at offset 0,0)
    @test ! any( neighbor.offset == CartesianIndex(0,0) for neighbor in testScatterPatch )

    # The point above  the origin is in the ScatterPatch ( at offset -1,0)
    # We know its the third element.
    # More generally, use any()
    @test testScatterPatch[2].offset == CartesianIndex(-1,0)

    # The point below the origin is not in the scatterPatch
    # because it is in the synth and has no value (not synthesized)
    @test ! any( neighbor.offset == CartesianIndex(1,0) for neighbor in testScatterPatch )

    # The point left of the origin is in the scatterPatch
    # because it is in the context and has a value
    @test any( neighbor.offset == CartesianIndex(0, -1) for neighbor in testScatterPatch )

    # The points are distinct
    # Specifically, the origin is not repeated
    # TODO

end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
using TreeKnit
using Documenter

Documenter.makedocs(;
	source = "src",
	clean = true,
	doctest = true,
	modules = Module[TreeKnit],
	repo = "",
	highlightsig = true,
	sitename = "TreeKnit documentation",
	expandfirst = [],
	pages = [
		"Index" => "index.md",
		"Usage" => [
			"Overview" => "overview.md",
			"MCCs" => "mccs.md",
			"Important options" => "options.md",
			# "Flu example" => "flu_example.md",
		],
		"Under the hood" => [
			"`opttrees`" => "opttrees.md",
			"`runopt`" => "runopt.md",
			"Resolving" => "resolving.md",
		],
		"Library" => [
			"OptArgs" => "types.md",
			"Functions" => "functions.md",
		],
	]
)

deploydocs(
    repo = "github.com/PierreBarrat/TreeKnit.jl.git";
    versions = nothing,
)

ext: jl

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
@testset "imfill" begin
    #2 Dimensional case
    image = Bool[
        0 0 0 0 0 0 1 0 0 0
        0 1 1 1 1 1 0 0 0 0
        0 1 0 0 1 1 0 0 0 0
        0 1 1 1 0 1 0 0 0 0
        0 1 1 1 1 1 0 0 0 0
        0 0 0 0 0 0 0 1 1 1
        0 0 0 0 0 0 0 1 0 1
        0 0 0 0 0 0 0 1 1 1
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
    ]
    answer = Bool[
        0 0 0 0 0 0 1 0 0 0
        0 1 1 1 1 1 0 0 0 0
        0 1 1 1 1 1 0 0 0 0
        0 1 1 1 1 1 0 0 0 0
        0 1 1 1 1 1 0 0 0 0
        0 0 0 0 0 0 0 1 1 1
        0 0 0 0 0 0 0 1 1 1
        0 0 0 0 0 0 0 1 1 1
        0 0 0 0 0 0 0 0 0 0
        0 0 0 0 0 0 0 0 0 0
    ]
    @test imfill(.!image, (0, 2)) == (.!answer)

    #3 Dimensional case
    A = zeros(Bool, 3, size(image)[1], size(image)[2])
    A[1, :, :] = image
    A[2, :, :] = image
    A[3, :, :] = image
    B = zeros(Bool, 3, size(answer)[1], size(answer)[2])
    B[1, :, :] = answer
    B[2, :, :] = answer
    B[3, :, :] = answer
    @test imfill(.!A, (0, 8)) == (.!B)

    # Complete white image and complete black image case
    img = zeros(Bool, 10, 10)
    @test imfill(img, (0, 10)) == img
    @test imfill(.!img, (0, 99)) == .!img
    # Miscellaneous case
    img[5, 5] = true
    @test imfill(.!img, (0, 10)) == .!img
    @test imfill(.!img, (0, 99)) == zeros(Bool, 10, 10)
end


text:
# adapted from https://www.youtube.com/watch?v=vAp6nUMrKYg
struct Dual <: Number
    x::Float64
    y::Float64
end

import Base: +, -, *, /, convert, promote_rule, show

+(a::Dual, b::Dual) = Dual(a.x+b.x, a.y+b.y)
-(a::Dual, b::Dual) = Dual(a.x-b.x, a.y-b.y)

*(a::Dual, b::Dual) = Dual(a.x*b.x, a.y*b.x + a.x*b.y)
/(a::Dual, b::Dual) = Dual(a.x/b.x, (a.y*b.x - a.x*b.y)/(b.x^2))

convert(::Type{Dual}, x::Real) = Dual(x, zero(x))
promote_rule(::Type{Dual}, ::Type{<:Number}) = Dual
show(io::IO, d::Dual) = print(io, d.x, " + ", d.y," ϵ")

f(x) = x^2
f(1)
f(1.0)
f(Dual(1,1))

"""
Compute the square root of x using the Babylonian algorithm
"""
function babylonian(x; nmax = 10)
    t = (1+x)/2
    for i in 2:nmax
        t = (t + x/t)/2
    end
    return t
end

babylonian(π), √π
babylonian(2), √2

# What happens when we pass a Dual number to babylonian?
x = Dual(2,1)
# We automagically get the derivative!
babylonian(x)
√2, 0.5/√2

max_stars_repo_name: csimal/Julia-Unamur

max_issues_repo_name: csimal/Julia-Unamur

max_forks_repo_name: csimal/Julia-Unamur


text:
module Relations

using HyperRectangles

export before, meets, overlaps, starts, during, finishes

# http://en.wikipedia.org/wiki/Allen%27s_interval_algebra

function before{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    for i = 1:N
        b1.max[i] < b2.min[i] || return false
    end
    true
end

function meets{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    b1.max == b2.min
end

function overlaps{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    for i = 1:N
        b1.max[i] > b2.min[i] && b1.min[i] > b2.min[i] || return false
    end
    true
end

function starts{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    b1.min == b2.min
end

function during{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    for i = 1:N
        b1.max[i] < b2.max[i] && b1.min[i] > b2.min[i] || return false
    end
    true
end

function finishes{T1, T2, N}(b1::HyperRectangle{T1,N}, b2::HyperRectangle{T2, N})
    b1.max == b2.max
end

end # module


text:
using Crayons
using Crayons.Box
using Test

withenv("FORCE_COLOR" => true) do

Crayons.print_logo()

# Styles
@test string(Crayon(bold = true)) == "\e[1m"
@test string(Crayon(bold = false)) == "\e[22m"
@test string(Crayon(underline = true)) == "\e[4m"
@test string(Crayon(underline = false)) == "\e[24m"
@test string(Crayon(bold = true, underline = false)) == "\e[1;24m"
@test string(Crayon(italics = false)) == "\e[23m"
@test string(Crayon(italics = true, underline = false)) == "\e[3;24m"
@test string(Crayon(faint = true)) == "\e[2m"
@test string(Crayon(conceal = true)) == "\e[8m"
@test string(Crayon(faint = false)) == "\e[22m"
@test string(Crayon(conceal = false)) == "\e[28m"
@test string(Crayon(strikethrough = true)) == "\e[9m"
@test string(Crayon(blink = true)) == "\e[5m"
@test string(Crayon(strikethrough = false)) == "\e[29m"
@test string(Crayon(blink = false)) == "\e[25m"
@test string(Crayon(negative = true)) == "\e[7m"
@test string(Crayon(negative = false)) == "\e[27m"
@test string(Crayon(reset = true)) == "\e[0m"
@test string(Crayon(reset = false)) == ""
@test string(Crayon()) == ""

# 16 colors
@test string(Crayon(foreground = :red)) == "\e[31m"
@test string(Crayon(background = :green)) == "\e[42m"
@test string(Crayon(foreground = :red, background = :green)) == "\e[31;42m"
@test string(Crayon(foreground = :red, bold = true)) == "\e[31;1m"
@test string(Crayon(foreground = :default, background = :default)) == "\e[39;49m"

# 256 colors
@test string(Crayon(foreground = 32)) == "\e[38;5;32m"
@test string(Crayon(foreground = 32, bold = true)) == "\e[38;5;32;1m"
@test string(Crayon(foreground = 32, background = 15, bold = true)) == "\e[38;5;32;48;5;15;1m"

# 24 bit colors
@test string(Crayon(foreground = (100,150,200))) == "\e[38;2;100;150;200m"
@test string(Crayon(foreground = (100,150,200), background = (200,210,220))) == "\e[38;2;100;150;200;48;2;200;210;220m"
@test string(Crayon(foreground = (100,150,200), background = (200,210,220), bold = true)) == "\e[38;2;100;150;200;48;2;200;210;220;1m"

# Mixing color modes
@test string(Crayon(foreground = :red, background = (10,20,30))) == "\e[31;48;2;10;20;30m"
@test string(Crayon(foreground = :60, background = (10,20,30))) == "\e[38;5;60;48;2;10;20;30m"

# Hex colors
@test string(Crayon(foreground = 0x000000)) ==  string(Crayon(foreground = (0,0,0)))
@test string(Crayon(foreground = 0xffffff)) ==  string(Crayon(foreground = (255,255,255)))
@test string(Crayon(foreground = 0xffaadd)) ==  string(Crayon(foreground = (255,170,221)))

# String macro
@test string(crayon"0xffffff") == string(Crayon(foreground = 0xffffff))
@test string(crayon"ffffff") == string(Crayon(foreground = 0xffffff))
@test string(crayon"#ffffff") == string(Crayon(foreground = 0xffffff))
@test string(crayon"fg:0xffffff") == string(Crayon(foreground = 0xffffff))
@test string(crayon"bg:0xff00ff fg:0xffffff") == string(Crayon(foreground = 0xffffff, background = 0xff00ff))
@test string(crayon"bg:red bold !underline") == string(Crayon(background = :red, bold = true, underline = false))
@test string(crayon"bg:(1,2,3) fg:(2,1,5)") == string(Crayon(background = (1,2,3), foreground = (2,1,5)))

#
withenv("FORCE_256_COLORS" => true) do
    @test string(crayon"(0,0,255)") == string(Crayon(foreground = 21))
    @test string(crayon"fg:(0,0,255) bg:(255,0,255)") == string(Crayon(foreground = 21, background = 201))
end

# CrayonStack
cs = CrayonStack()
@test string(cs) == string(Crayon(foreground = :default, background = :default, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
push!(cs, Crayon(foreground = :red))
@test string(cs) == string(Crayon(foreground = :red, background = :default, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
push!(cs, Crayon(foreground = :green))
@test string(cs) == string(Crayon(foreground = :green, background = :default, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
push!(cs, Crayon(bold = true, background = 125))
@test string(cs) == string(Crayon(foreground = :green, background = 125, bold = true, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
push!(cs, Crayon(bold = false))
@test string(cs) == string(Crayon(foreground = :green, background = 125, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
pop!(cs) # Popping the bold = false
@test string(cs) == string(Crayon(foreground = :green, background = 125, bold = true, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
pop!(cs) # Popping the bold = true + background = 125
@test string(cs) == string(Crayon(foreground = :green, background = :default, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
pop!(cs) # Popping the foreground = :green
@test string(cs) == string(Crayon(foreground = :red, background = :default, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
pop!(cs) # Popping the foreground = :red
@test string(cs) == string(Crayon(foreground = :default, background = :default, bold = false, italics = false, underline = false, strikethrough = false, blink = false, conceal = false, negative = false, faint = false))
@test_throws ArgumentError pop!(cs)

# Incremental mode CrayonStack
cs = CrayonStack(incremental = true) # State:
push!(cs, Crayon(foreground = :red)) # State: fg = red, change: fg = red
@test string(cs) == string(Crayon(foreground = :red))
push!(cs, Crayon(foreground = :red)) # State: fg = red, change:
@test string(cs) == ""
push!(cs, Crayon(foreground = :blue)) # State: fg = blue, change: fg = blue
@test string(cs) == string(Crayon(foreground = :blue))
push!(cs, Crayon(bold = true)) # State fg = blue, bold = true, change: bold = true
@test string(cs) == string(Crayon(bold = true))
pop!(cs) # State fg = blue, change: bold = false
@test string(cs) == string(Crayon(bold = false))
pop!(cs) # State fg = red, change: fg = red
@test string(cs) == string(Crayon(foreground = :red))
pop!(cs) # State fg = red, change:
@test string(cs) == ""
pop!(cs) # State change: fg = default
@test string(cs) == string(Crayon(foreground = :default))

# Merge
@test string(merge(Crayon(foreground = :blue, background = :red))) == string(Crayon(foreground = :blue, background = :red))
@test string(merge(Crayon(foreground = :blue), Crayon(background = :red)))  == string(Crayon(foreground = :blue, background = :red))
@test string(merge(Crayon(foreground = :blue), Crayon(background = :red), Crayon(bold = true)))  == string(Crayon(foreground = :blue, background = :red, bold = true))
@test string(merge(Crayon(foreground = :blue), Crayon(background = :red), Crayon(bold = true)))  == string(Crayon(foreground = :blue, background = :red, bold = true))
@test string(merge(Crayon(foreground = :red), Crayon(foreground = :blue))) == string(Crayon(foreground = :blue))
@test string(merge(Crayon(foreground = :red), Crayon(negative = true))) == string(Crayon(foreground = :red, negative = true))

string(BLACK_BG * WHITE_FG * BOLD) == string(Crayon(foreground = :white, background = :black, bold = true))

# Call overloading
@test string(Crayon()("hello")) == "hello"
@test string(Crayon()(split("hello world")[1])) == "hello" # test substrings
@test string(Crayon(bold=true)("hello")) == string(BOLD, "hello", inv(BOLD))
@test string(Crayon(bold=true, foreground = :red)("hello")) == string(Crayon(foreground=:red, bold=true), "hello", Crayon(foreground=:default, bold=false))

@test string("normal", BOLD("bold", RED_FG("red bold", UNDERLINE("red_underline"), "red bold"), "bold"), "normal" ) ==
    string("normal", BOLD, "bold", RED_FG, "red bold", UNDERLINE, "red_underline", inv(UNDERLINE), "red bold", inv(RED_FG), "bold", inv(BOLD), "normal")

@test string("normal", BOLD*UNDERLINE("bold_underline", ITALICS*RED_FG("everything"), "bold_underline"), "normal") ==
    string("normal", Crayon(bold=true, underline=true),"bold_underline", Crayon(italics=true, foreground=:red), "everything", Crayon(italics=false, foreground=:default),
     "bold_underline", Crayon(bold=false, underline=false), "normal")

Crayons.test_system_colors(IOBuffer())
Crayons.test_256_colors(IOBuffer())
Crayons.test_24bit_colors(IOBuffer())

@test Crayons.to_256_colors(crayon"0000ff") == Crayon(foreground = 21)
@test Crayons.to_256_colors(crayon"00ff00") == Crayon(foreground = 46)
@test Crayons.to_256_colors(crayon"ff0000") == Crayon(foreground = 196)

@test Crayons.to_system_colors(crayon"0000ff") == Crayon(foreground = :light_blue)
@test Crayons.to_system_colors(crayon"00ff00") == Crayon(foreground = :light_green)
@test Crayons.to_system_colors(crayon"ff0000") == Crayon(foreground = :light_red)

end # withenv


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the wedge(x, y) function
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using LinearAlgebra: norm
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

@testset "wedge(M::Multivector, N::Multivector)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "wedge(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "wedge(B::Blade, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::Blade, C::Blade)" begin
    # --- Preparations

    B_vectors = hcat([1; 1; 0; 0; 0],
                     [0; 2; 0; 0; 0])
    B = Blade(B_vectors)

    C_vectors = hcat([0; 1; 3; 0; 0],
                     [0; 0; 0; 4; 0])
    C = Blade(C_vectors)

    # --- Tests

    # ------ dim(B) == dim(C)

    # grade(B) + grade(C) <= dim(B)
    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(hcat(B_vectors, C_vectors))
    @test B ∧ C == B_wedge_C

    # grade(B) + grade(C) > dim(B)
    C_vectors = hcat([1; 0; 0; 0; 0],
                     [0; 1; 0; 0; 0],
                     [0; 0; 1; 0; 0],
                     [0; 0; 0; 1; 0])
    C = Blade(C_vectors)
    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)

    # ------ dim(B) != dim(C)

    C = Blade(rand(dim(B) + 1, 2))
    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

@testset "wedge(B::Blade, C::Pseudoscalar)" begin
    # --- Preparations

    test_dim = 10

    test_vectors = rand(test_dim, 3)

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Test

    # dim(B) == dim(C)
    B = Blade(test_vectors)
    C = Pseudoscalar(test_dim, test_value)

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)

    # dim(B) != dim(C)
    B = Blade(test_vectors)
    C = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

@testset "wedge(B::Blade, C::Scalar)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Scalar(test_value)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(basis(B), volume=value(C) * volume(B))
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Blade, C::One)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))

    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Blade, C::Zero)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))

    C = Zero()

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Blade, C::Real)" begin
    test_dim = 10
    B = Blade(rand(test_dim, 3))

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(basis(B), volume=C * volume(B))
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Blade, C::Vector)" begin
    # --- Preparations

    B_vectors = hcat([1; 0; 0; 0; 0],
                     [0; 2; 0; 0; 0])
    B = Blade(B_vectors)

    C = [0; 1; 3; 0; 0]

    # --- Tests

    # dim(B) == dim(C)
    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(hcat(B_vectors, C))
    @test B ∧ C == B_wedge_C

    # dim(B) != dim(C)
    C = Vector(rand(dim(B) + 1))
    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

# ------ B::Pseudoscalar

@testset "wedge(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::Pseudoscalar, C::Blade)" begin
    # --- Preparations

    test_dim = 10

    test_vectors = rand(test_dim, 3)

    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Test

    # dim(B) == dim(C)
    B = Pseudoscalar(test_dim, test_value)
    C = Blade(test_vectors)

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)

    # dim(B) != dim(C)
    B = Pseudoscalar(test_dim + 1, test_value)
    C = Blade(test_vectors)
    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

@testset "wedge(B::Pseudoscalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_dim = 10

    test_value_1 = get_random_value(1)  # add 1 to keep value away from 0
    test_value_2 = get_random_value(1)  # add 1 to keep value away from 0

    # --- Tests

    # dim(B) == dim(C)
    B = Pseudoscalar(test_dim, test_value_1)
    C = Pseudoscalar(test_dim, test_value_2)

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)

    # dim(B) != dim(C)
    B = Pseudoscalar(test_dim, test_value_1)
    C = Pseudoscalar(test_dim + 1, test_value_2)
    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

@testset "wedge(B::Pseudoscalar, C::Scalar)" begin
    # --- Preparations

    test_dim = 10
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C == Pseudoscalar(test_dim, test_value_1 * test_value_2)
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Pseudoscalar, C::One)" begin
    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Pseudoscalar, C::Zero)" begin
    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Pseudoscalar(test_dim, test_value)

    C = Zero()

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Pseudoscalar, C::Real)" begin
    # --- Preparations

    test_dim = 10
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Pseudoscalar(test_dim, test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    B_wedge_C = wedge(B, C)
    @test B_wedge_C == Pseudoscalar(test_dim, test_value_1 * test_value_2)
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Pseudoscalar, C::Vector)" begin
    # --- Preparations

    test_dim = 10

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Tests

    # dim(B) == dim(C)
    for test_dim in 5:8
        B = Pseudoscalar(test_dim, test_value)
        C = Vector(rand(test_dim))

        @test iszero(wedge(B, C))
        @test iszero(B ∧ C)
    end

    # dim(B) != dim(C)
    B = Pseudoscalar(test_dim + 1, test_value)
    C = Vector(rand(test_dim))

    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

# ------ B::Scalar

@testset "wedge(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::Scalar, C::Blade)" begin
    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Scalar(test_value)

    test_dim = 10
    C = Blade(rand(test_dim, 3))

    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(basis(C), volume=value(B) * volume(C))
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Scalar, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Scalar(test_value_1)

    test_dim = 10
    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C == Pseudoscalar(test_dim, test_value_1 * test_value_2)
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Scalar, C::Scalar)" begin
    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value_1)

    test_value_2 = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value_2)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C isa Scalar
    @test B_wedge_C == test_value_1 * test_value_2
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::Scalar, C::One)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::Scalar, C::Zero)" begin
    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    B = Scalar(test_value)

    C = Zero()

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Scalar, C::Real)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = Scalar(test_value_1)

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = test_value_2

    B_wedge_C = wedge(B, C)
    @test B_wedge_C isa AbstractScalar
    @test B_wedge_C == test_value_1 * test_value_2
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Scalar, C::Vector)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = Scalar(test_value)

    C = Vector(rand(5))

    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(C, volume=norm(C) * value(B))
    @test B ∧ C == B_wedge_C
end

# ------ B::One

@testset "wedge(B::One, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::One::, C::Blade)" begin
    B = One()

    test_dim = 10
    C = Blade(rand(test_dim, 3))

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === C
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::One::, C::Pseudoscalar)" begin
    B = One()

    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === C
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::One, C::Scalar)" begin
    B = One()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === C
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::One, C::One)" begin
    B = One()
    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::One, C::Zero)" begin
    B = One()
    C = Zero()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === C
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::One, C::Real)" begin
    B = One()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    B_wedge_C = wedge(B, C)
    @test B_wedge_C isa AbstractScalar
    @test wedge(B, C) == C
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::One, C::Vector)" begin
    B = One()
    C = Vector(rand(5))

    B_wedge_C = wedge(B, C)
    @test B_wedge_C == Blade(C)
    @test B ∧ C == B_wedge_C
end

# ------ B::Zero

@testset "wedge(B::Zero, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::Zero, C::Blade)" begin
    B = Zero()

    test_dim = 10
    C = Blade(rand(test_dim, 3))

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Zero, C::Pseudoscalar)" begin
    B = Zero()

    test_dim = 10
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Pseudoscalar(test_dim, test_value)

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Zero, C::Scalar)" begin
    B = Zero()

    test_value = get_random_value(2)  # add 2 to keep value away from 0 and 1
    C = Scalar(test_value)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::Zero, C::One)" begin
    B = Zero()
    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::Zero, C::Zero)" begin
    B = Zero()
    C = Zero()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C === B
    @test (B ∧ C) == B_wedge_C
end

@testset "wedge(B::Zero, C::Real)" begin
    B = Zero()

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = test_value

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Zero, C::Vector)" begin
    B = Zero()
    C = Vector(rand(5))

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

# ------ B::Real

@testset "wedge(B::Real, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::Real, C::Blade)" begin
    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    test_dim = 10
    C = Blade(rand(test_dim, 3))

    # Exercise functionality and check results
    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(basis(C), volume=B * volume(C))
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Real, C::Pseudoscalar)" begin
    # --- Preparations

    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = test_value_1

    test_dim = 10
    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Pseudoscalar(test_dim, test_value_2)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C == Pseudoscalar(test_dim, test_value_1 * test_value_2)
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Real, C::Scalar)" begin
    test_value_1 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_1 = rand() > 0.5 ? test_value_1 : -test_value_1
    B = test_value_1

    test_value_2 = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value_2 = rand() > 0.5 ? test_value_2 : -test_value_2
    C = Scalar(test_value_2)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C isa AbstractScalar
    @test B_wedge_C == test_value_1 * test_value_2
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Real, C::One)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C isa AbstractScalar
    @test B_wedge_C == test_value
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Real, C::Zero)" begin
    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    B = test_value

    C = Zero()

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

# ------ B::Vector

@testset "wedge(B::Vector, M::Multivector)" begin
    @test_skip 1
end

@testset "wedge(B::Vector, C::Blade)" begin
    # --- Preparations

    B = [0; 1; 3; 0; 0]

    C_vectors = hcat([1; 0; 0; 0; 0],
                     [0; 2; 0; 0; 0])
    C = Blade(C_vectors)

    # --- Tests

    # dim(B) == dim(C)
    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(hcat(B, C_vectors))
    @test B ∧ C == B_wedge_C

    # dim(B) != dim(C)
    B = Vector(rand(dim(C) + 1))
    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

@testset "wedge(B::Vector, C::Pseudoscalar)" begin
    # --- Preparations

    test_dim = 10

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Tests

    # dim(B) == dim(C)
    for test_dim in 5:8
        B = Vector(rand(test_dim))
        C = Pseudoscalar(test_dim, test_value)

        @test iszero(wedge(B, C))
        @test iszero(B ∧ C)
    end

    # dim(B) != dim(C)
    B = Vector(rand(test_dim))
    C = Pseudoscalar(test_dim + 1, test_value)

    @test_throws DimensionMismatch wedge(B, C)
    @test_throws DimensionMismatch B ∧ C
end

@testset "wedge(B::Vector, C::Scalar)" begin
    B = Vector(rand(5))

    test_value = rand() + 2  # add 2 to keep value away from 0 and 1
    test_value = rand() > 0.5 ? test_value : -test_value
    C = Scalar(test_value)

    B_wedge_C = wedge(B, C)
    @test B_wedge_C ≈ Blade(B, volume=norm(B) * value(C))
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Vector, C::One)" begin
    B = Vector(rand(5))
    C = One()

    B_wedge_C = wedge(B, C)
    @test B_wedge_C == Blade(B)
    @test B ∧ C == B_wedge_C
end

@testset "wedge(B::Vector, C::Zero)" begin
    B = Vector(rand(5))
    C = Zero()

    @test iszero(wedge(B, C))
    @test iszero(B ∧ C)
end

@testset "wedge(B::Vector, C::Vector)" begin
    # --- Preparations

    B_vector = Vector{Float64}([1, 2, 3, 0, 0])
    B = Blade(B_vector)

    C_vector = [0; 1; 3; 0; 0]
    C = Blade(C_vector)

    # --- Tests

    # dim(B) == dim(C)
    B_wedge_C = wedge(B_vector, C_vector)
    @test B_wedge_C ≈ Blade(hcat(B_vector, C_vector))
    @test B_vector ∧ C_vector == B_wedge_C

    # dim(B) != dim(C)
    C_vector = Vector(rand(dim(B) + 1))
    @test_throws DimensionMismatch wedge(B_vector, C_vector)
    @test_throws DimensionMismatch B_vector ∧ C_vector
end


text:
p = ℝ³(rand(3))
q = Quaternion(rand() * 2pi, normalize(ℝ³(rand(3))))
q₀ = Quaternion(0, normalize(ℝ³(rand(3))))
q2 = normalize(Quaternion(rand(4)))
points = [ℝ³(rand(3)) for i in 1:5]
i = ℝ³(1, 0, 0)
j = ℝ³(0, 1, 0)
k = ℝ³(0, 0, 1)
g = getrotation(i, j)

@test typeof(rotate(p, q)) == typeof(p)
@test isapprox(rotate(p, q₀), p) # rotations with angle zero
@test typeof(rotate(points, q)) == typeof(points) # Array input
@test isapprox(g, Quaternion(pi / 2, k))

q = Biquaternion(Quaternion(rand() * 2pi - pi, ℝ³(rand(3))), ℝ³(rand(3)))
transformed = applyconfig(points, q)

@test typeof(transformed) == typeof(points)
@test size(transformed) == size(points)
@test isapprox(transformed, points) == false


text:
module CellFieldsTests

end # module


text:
abstract AbstractCoordinate
import Base: +

immutable HexGridCoordinate{T} <: AbstractCoordinate
  x::T
  y::T
end

Base.zero{T}(::Type{HexGridCoordinate{T}}) = HexGridCoordinate(zero(T), zero(T))
+{T}(a::HexGridCoordinate{T}, b::HexGridCoordinate{T}) = HexGridCoordinate(a.x + b.x, a.y + b.y)
Base.print(io::IO, c::HexGridCoordinate) = print(io, c.x, ", ", c.y)

# typealias for the almost exclusively used integer grid positions
typealias Position HexGridCoordinate{Int}

function Position(d::Dict)
  return Position(Int(d["x"]), Int(d["y"]))
end

to_dict(p::HexGridCoordinate) = Dict("x" => p.x, "y" => p.y)
# this seems a little strange, but it means that we can now easily build functions
# that takes positional objects or positions as parameters and just call this function to
# get the actual position
position(p::AbstractCoordinate) = p
distance(a::HexGridCoordinate, b::HexGridCoordinate) = distance(CubeGridCoordinate(a), CubeGridCoordinate(b))

# ######################################################
# cube coordinate position helper type for internal
# computations on the hexagonal grid
########################################################
immutable CubeGridCoordinate{T} <: AbstractCoordinate
  x::T
  y::T
  z::T
end

CubeGridCoordinate(p::HexGridCoordinate) = CubeGridCoordinate(p.x, -p.x - p.y, p.y)
HexGridCoordinate(p::CubeGridCoordinate) = HexGridCoordinate(p.x, p.z)
function distance(a::CubeGridCoordinate, b::CubeGridCoordinate)
  return max(abs(a.x-b.x), abs(a.y-b.y), abs(a.z-b.z))
end

Base.zero{T}(::Type{CubeGridCoordinate{T}}) = HexGridCoordinate(zero(T), zero(T), zero(T))

########################################################
# type for circular areas in hex grid
########################################################

abstract AbstractArea
abstract CircularArea <: AbstractArea

# circular area that is centered around the origin
immutable OriginCircularArea <: CircularArea
  radius::Int
end

radius(a::CircularArea) = a.radius
center(a::OriginCircularArea) = Position(0, 0)
circle(r::Int) = OriginCircularArea(r)

Base.in(p::Position, a::CircularArea) = distance(p, center(a)) <= radius(a)

Base.eltype(::Type{AbstractArea}) = Position
Base.start(a::CircularArea) = Position(-radius(a), 0)
Base.length(iter::CircularArea) = 1 + 6div(radius(iter) * (radius(iter) + 1), 2)

function Base.next(iter::CircularArea, state::Position)
  rad = radius(iter)
  if state.y < min(rad, -state.x+rad)
    nxt = Position(state.x, state.y + 1)
  else
    nxt = Position(state.x + 1, max(-rad, -1-state.x-rad))
  end
  return state + center(iter), nxt
end

function Base.done(iter::CircularArea, state::Position)
  return state.x > radius(iter)
end

# this type defines a circular area in a hexagonal grid.
# it is equipped with an iterator interface
immutable GeneralCircularArea <: CircularArea
  center::Position
  radius::Int
end

radius(a::GeneralCircularArea) = a.radius
center(a::GeneralCircularArea) = a.center
circle(o::Position, r::Int) = GeneralCircularArea(o, r)
circle(origin, radius::Int) = circle(position(origin), radius)

# shift circles
+(a::CircularArea, b::Position) = GeneralCircularArea(center(a)+b, radius(a))

# type for intersection of spherical shape and a complete map
immutable IntersectionArea{T1<:AbstractArea, T2 <: AbstractArea} <: AbstractArea
  A::T1
  B::T2
end

immutable IntIterState
  value::Position
  state::Position
  done::Bool
end

function Base.start(iter::IntersectionArea)
  state = start(iter.A)
  value, state = next(iter.A, state)
  while value ∉ iter.B && !done(iter.A, state)
    value, state = next(iter.A, state)
  end
  return IntIterState(value, state, value ∉ iter.B)
end
# cannot predict length

function Base.next(iter::IntersectionArea, state::IntIterState)
  # return the last value now
  if done(iter.A, state.state)
    return state.value, IntIterState(state.value, state.state, true)
  end

  v, n = next(iter.A, state.state)
  while v ∉ iter.B && !done(iter.A, n)
    v, n = next(iter.A, n)
  end
  return state.value, IntIterState(v, n, v ∉ iter.B)
end

Base.done(iter::IntersectionArea, state::IntIterState) = state.done
Base.intersect(a::AbstractArea, b::AbstractArea) = IntersectionArea(a, b)

###################################################
#        convenience functions for bots
###################################################
# centered areas
view_area(config::Config)    = circle(config.see)
radar_area(config::Config)   = circle(config.radar)
move_area(config::Config)    = circle(config.move)
damage_area(config::Config)  = circle(config.cannon)
map_area(config::Config)     = circle(config.field_radius)

# areas around a specific point
view_area(origin,  config::Config)   = intersect(circle(origin, config.see),    map_area(config))
radar_area(origin, config::Config)   = intersect(circle(origin, config.radar),  map_area(config))
move_area(origin,  config::Config)   = intersect(circle(origin, config.move),   map_area(config))
damage_area(origin,  config::Config) = intersect(circle(origin, config.cannon), map_area(config))

max_stars_repo_name: Condimentum/spacetyckiting

max_issues_repo_name: Condimentum/spacetyckiting

max_forks_repo_name: Condimentum/spacetyckiting


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

mutable struct IoK8sApiCoreV1CSIPersistentVolumeSource <: SwaggerModel
    controllerExpandSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: controllerExpandSecretRef
    controllerPublishSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: controllerPublishSecretRef
    driver::Any # spec type: Union{ Nothing, String } # spec name: driver
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    nodePublishSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: nodePublishSecretRef
    nodeStageSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretReference } # spec name: nodeStageSecretRef
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    volumeAttributes::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: volumeAttributes
    volumeHandle::Any # spec type: Union{ Nothing, String } # spec name: volumeHandle

    function IoK8sApiCoreV1CSIPersistentVolumeSource(;controllerExpandSecretRef=nothing, controllerPublishSecretRef=nothing, driver=nothing, fsType=nothing, nodePublishSecretRef=nothing, nodeStageSecretRef=nothing, readOnly=nothing, volumeAttributes=nothing, volumeHandle=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("controllerExpandSecretRef"), controllerExpandSecretRef)
        setfield!(o, Symbol("controllerExpandSecretRef"), controllerExpandSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("controllerPublishSecretRef"), controllerPublishSecretRef)
        setfield!(o, Symbol("controllerPublishSecretRef"), controllerPublishSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("driver"), driver)
        setfield!(o, Symbol("driver"), driver)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("nodePublishSecretRef"), nodePublishSecretRef)
        setfield!(o, Symbol("nodePublishSecretRef"), nodePublishSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("nodeStageSecretRef"), nodeStageSecretRef)
        setfield!(o, Symbol("nodeStageSecretRef"), nodeStageSecretRef)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("volumeAttributes"), volumeAttributes)
        setfield!(o, Symbol("volumeAttributes"), volumeAttributes)
        validate_property(IoK8sApiCoreV1CSIPersistentVolumeSource, Symbol("volumeHandle"), volumeHandle)
        setfield!(o, Symbol("volumeHandle"), volumeHandle)
        o
    end
end # type IoK8sApiCoreV1CSIPersistentVolumeSource

const _property_map_IoK8sApiCoreV1CSIPersistentVolumeSource = Dict{Symbol,Symbol}(Symbol("controllerExpandSecretRef")=>Symbol("controllerExpandSecretRef"), Symbol("controllerPublishSecretRef")=>Symbol("controllerPublishSecretRef"), Symbol("driver")=>Symbol("driver"), Symbol("fsType")=>Symbol("fsType"), Symbol("nodePublishSecretRef")=>Symbol("nodePublishSecretRef"), Symbol("nodeStageSecretRef")=>Symbol("nodeStageSecretRef"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("volumeAttributes")=>Symbol("volumeAttributes"), Symbol("volumeHandle")=>Symbol("volumeHandle"))
const _property_types_IoK8sApiCoreV1CSIPersistentVolumeSource = Dict{Symbol,String}(Symbol("controllerExpandSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("controllerPublishSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("driver")=>"String", Symbol("fsType")=>"String", Symbol("nodePublishSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("nodeStageSecretRef")=>"IoK8sApiCoreV1SecretReference", Symbol("readOnly")=>"Bool", Symbol("volumeAttributes")=>"Dict{String, String}", Symbol("volumeHandle")=>"String")
Base.propertynames(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1CSIPersistentVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1CSIPersistentVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1CSIPersistentVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1CSIPersistentVolumeSource)
    (getproperty(o, Symbol("driver")) === nothing) && (return false)
    (getproperty(o, Symbol("volumeHandle")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1CSIPersistentVolumeSource }, name::Symbol, val)
end


text:
# This file was generated, do not modify it. # hide
using DelimitedFiles

function cnt(ReZ, ImZ, ReC, ImC)
    k::UInt8 = 0
    while (k < 255)
        new_ReZ = ReZ * ReZ - ImZ * ImZ + ReC
        new_ImZ = 2 * ReZ * ImZ + ImC
        if new_ReZ * new_ReZ + new_ImZ * new_ImZ > 4.
            return k
        end
        ReZ = new_ReZ
        ImZ = new_ImZ
        k += 1
    end
    return k
end

function mand(M,N)
    grid = Array{UInt8}(undef, M, N)
    init_Re = 0.0
    init_Im = 0.0
    for j = 1:N
        for i = 1:M
            ReC = 4.0 / (M - 1) * (j - 1) - 2.0
            ImC = 4.0 / (N - 1) * (i - 1) - 2.0
            ReZ = init_Re
            ImZ = init_Im
            grid[i, j] = cnt(ReZ, ImZ, ReC, ImC)
        end
    end
    return grid
end

function main()
    s = time()
    M = N = 4200
    grid = mand(M,N)
    t = time()
    @show(t - s)
    writedlm("/tmp/resultJ.txt",grid,",")
end

main()

using BenchmarkTools
@btime mand(5000, 5000);


text:
using QuantumOptics
using MultimodeCavity
const mm = MultimodeCavity
using Optim

# Particles
const E0 = 0.2
const Nparticles = 1
const Nlevels = 10
const potential = mm.quantum.BoxPotential
const particles = mm.quantum.Bosons(Nparticles, Nlevels, E0, potential)

# Cavity mode
const index1 = 2
const cavitymodes1 = 12
const delta1 = -1.
const eta1 = 4.0
const kappa1 = 1.
const U0_1 = 0
const mode1 = mm.quantum.CavityMode(index1, cavitymodes1, delta1, eta1, kappa1, U0_1)

# Multimode system
const scaling = 1.
const system = mm.quantum.MultimodeSystem(scaling, particles, [mode1])

# Bases
basis_total = mm.quantum.basis(system)
basis_particles = mm.quantum.basis(particles)
basis_mode1 = mm.quantum.basis(mode1)

# Initial state
const ψ₀ = basis_ket(basis_particles, 1) ⊗ basis_ket(basis_mode1, 1)
const ρ₀ = ψ₀ ⊗ dagger(ψ₀)

T = [0.,10.0]

tout, ρt = mm.quantum.timeevolution_master(system, T, ρ₀, reltol=1e-7)
ρp = ptrace(ρt[end], [2])
ρf = ptrace(ρt[end], [1])
ρf /= trace(ρf)

using PyCall
@pyimport matplotlib.pyplot as plt

H = mm.quantum.Hamiltonian(system)
U, S, V = svd(ρt[end].data)
x = Float64[-1:0.01:1;]

for i=1:9
    ϕ = Ket(basis_total, vec(V[:,i]))
    q = qfunc(ptrace(ϕ⊗dagger(ϕ), [1]), [-4:0.1:4;], [-4:0.1:4;])
    plt.figure(1)
    plt.subplot(3,3,i)
    plt.imshow(q)

    n = mm.quantum.particledensity(particles, x, ptrace(ϕ⊗dagger(ϕ), [2]))
    plt.figure(2)
    plt.subplot(3,3,i)
    plt.plot(x, n)
end

H = mm.quantum.Hamiltonian(particles)
U, S, V = svd(ρp.data)

for i=1:9
    ϕ = Ket(basis_particles, vec(V[:,i]))
    n = mm.quantum.particledensity(particles, x, ϕ⊗dagger(ϕ))
    plt.figure(3)
    plt.subplot(3,3,i)
    plt.plot(x, n)
end
plt.show()


text:
## functions that work for AbstractTimedata types and return the same
## output for all types (Timematr, Timenum,...)

####################################
## Timedata information retrieval ##
####################################

function idx(tn::AbstractTimedata)
    return tn.idx
end

import DataFrames.names
function names(tn::AbstractTimedata)
    return names(tn.vals)
end

function core(td::AbstractTimedata)
    ## return all entries of Timedata object
    ## in order to force output other than Array{Any}, use asArr
    return get(td)
end

function core(td::AbstractTimematr)
    ## return Float64
    return convert(Array{Float64}, convert(Array, td.vals))
end

function idxtype(td::AbstractTimedata)
    return typeof(idx(td))
end

#########
## get ##
#########

import Base.get
function get(td::AbstractTimedata, idx1::Int, idx2::Int)
    ## return entry of Timedata object
    return td.vals[idx1, idx2]
end

function get(td::AbstractTimedata)
    ## return all entries of Timedata object
    return [get(td, ii, jj) for ii=1:size(td, 1), jj=1:size(td, 2)]
end

## extract single value
##---------------------

function get(td::AbstractTimedata, x::Colon)
    if size(td) != (1, 1)
        error("TimeData obj must be 1x1 for call with colon")
    end
    return get(td, 1, 1)
end

function get(df::DataFrame, x::Colon)
    if size(df) != (1, 1)
        error("DataFrame obj must be 1x1 for call with colon")
    end
    return df[1, 1]
end

function get(da::DataArray, x::Colon)
    if (size(da) != (1, 1)) | (size(da) != (1,))
        error("DataArray obj must be 1x1 for call with colon")
    end
    return da[1, 1]
end

function get(x::NAtype, y::Colon)
    return NA
end

function get{T}(x::Nullable{T}, y::Colon)
    return x
end

function get(x::Bool, y::Colon)
    return x
end

function get(x::Real, y::Colon)
    return x
end

## function get(x::Any, y::Colon)
##     return x
## end

###################
## Timedata size ##
###################

import Base.size
function size(tn::AbstractTimedata)
    return size(tn.vals)
end

function size(tn::AbstractTimedata, ind::Int)
    return size(tn.vals, ind)
end

import Base.ndims
function ndims(tn::AbstractTimedata)
    return ndims(tn.vals)
end

######################
## isequal function ##
######################

## isequal(NA, NA) -> true
import Base.isequal
function isequal(tn::AbstractTimedata, tn2::AbstractTimedata)
    ## return single true or false
    isequal(typeof(tn), typeof(tn2)) || return false
    isequal(tn.vals, tn2.vals) || return false
    isequal(tn.idx, tn2.idx) || return false
    return true
end

## NA == NA -> NA
import Base.==
function ==(tn::AbstractTimedata, tn2::AbstractTimedata)
    ## return single true or false
    isequal(typeof(tn), typeof(tn2)) || return false
    isequal(tn.idx, tn2.idx) || return false
    return ==(tn.vals, tn2.vals)
end

import Base.(.==)
function .==(tn, tn2)
    if !equMeta(tn, tn2)
        error("Timedata object must be similar for elementwise
isequal")
    end
    (nObs, nVars) = size(tn)
    df = DataFrame()
    for (nam, col) in eachcol(tn)
        df[nam] = tn.vals[nam] .== tn2.vals[nam]
    end
    return Timedata(df, idx(tn))
end

## isequal(NA, NA) -> true
function isequalElw(tn::AbstractTimedata, tn2::AbstractTimedata)
    ## return matrix with true / false elements
    if !equMeta(tn, tn2)
        error("Timedata object must be similar for elementwise
isequal")
    end
    (nObs, nVars) = size(tn)
    vals = Bool[isequal(tn.vals[ii, jj], tn2.vals[ii, jj]) for
        ii=1:nObs,
        jj=1:nVars]
    return Timedata(vals, names(tn), idx(tn))
end

import Base.isapprox
function isapprox(tn::AbstractTimedata, tn2::AbstractTimedata)
    ## return single true or false
    isequal(typeof(tn), typeof(tn2)) || return false

    (nObs, nVars) = size(tn)
    for ii=1:nObs
        for jj=1:nVars
            if !isequal(get(tn, ii, jj), get(tn2, ii, jj))
                if !isapprox(get(tn, ii, jj), get(tn2, ii, jj))
                    return false
                end
            end
        end
    end
    isequal(tn.idx, tn2.idx) || return false
    return true
end

##########
## isna ##
##########

function isnaElw(td::AbstractTimedata)
    ## elementwise test for NA
    df = DataFrame()
    for (nam, col) in eachcol(td.vals)
        df[nam] = [isna(col)...]
    end
    return Timedata(df, idx(td))
end

####################
## complete_cases ##
####################

import DataFrames.complete_cases
function complete_cases(td::AbstractTimedata)
    ## find days without NA
    return complete_cases(td.vals)
end

function complete_cases(tm::AbstractTimematr)
    ## no NAs allowed anyways
    return bool(ones(size(tm, 1)))
end

#####################
## equal structure ##
#####################

function equMeta(td1::AbstractTimedata, td2::AbstractTimedata)
    ## test if column and index metadata are exactly equal
    isequal(typeof(td1), typeof(td2)) || return false
    isequal(names(td1), names(td2)) || return false
    isequal(td1.idx, td2.idx) || return false
end

function similarMeta(td1::AbstractTimedata, td2::AbstractTimedata)
    ## test if either index or column metadata is equal
    ## subtraction for TimeData objects makes sense if either:
    ## - variables are equal: variable minus lagged version of
    ##   variable
    ## - index values are equal: one variable minus another variable

    ## variables must be of equal type
    isequal(super(typeof(td1)), super(typeof(td2))) || return false

    equDates = isequal(td1.idx, td2.idx)
    equColMeta = isequal(names(td1), names(td2))

    (equDates | equColMeta) || return false
end


text:
module Interpolation

include("spline.jl")

export SplineFunction
export ⋮, d2, ∫, ⩪

end # module


text:
function price_digits(symb)
    if symb[4:end] == "JPY"
        3
    else
        5
    end
end


text:
using Qaintessent
using LinearAlgebra
using SparseArrays: sparse
using Memoize

# Simple struct that represents a graph via its edges
struct Graph
    n::Int # number of vertices (indices 1,...,n)
    edges::Set{Set{Int}} # edges, represented as sets of vertices

    function Graph(n::Integer, edges::Vector{Tuple{T, T}}) where T <: Integer
        n >= 1 || throw(DomainError("n must be a positive integer"))

        # Turn into set of sets
        edge_set = Set(Set.(edges))

        # Verify that all edges are valid
        all(edge -> edge ⊆ 1:n && length(edge) == 2, edge_set) || throw(ArgumentError("Some edges have invalid endpoints"))
        new(n, edge_set)
    end
end

"""
    Phase separation gate for Max-κ-colorable subgraph QAOA mapping.
    Represents the objective function which counts the number of invalid
    edges (i.e. between vertices of the same color).
    Implemented the for one-hot encoding.

``U_{P}(\\gamma) = e^{-i \\gamma H_{P}}``
``H_{P} = \\sum_{\\{u, v\\} \\in E} \\sum_{a=1}^{\\kappa} Z_{u, a} Z_{v, a}``

Reference:\n
    Stuart Hadfield, Zhihui Wang, Bryan O'Gorman, Eleanor G. Rieffel, Davide Venturelli and Rupak Biswas\n
    From the Quantum Approximate Optimization Algorithm to a Quantum Alternating Operator Ansatz\n
    Algorithms 12.2 (2019), p.34
"""
struct MaxKColSubgraphPhaseSeparationGate <: AbstractGate
    # use a reference type (array with 1 entry) for compatibility with Flux
    γ::Vector{Float64}
    κ::Int # the number of possible colors
    graph::Graph # the underlying graph which should be colored

    function MaxKColSubgraphPhaseSeparationGate(γ::Float64, κ::Integer, graph::Graph)
        κ > 0 || throw(ArgumentError("Parameter `κ` must be a positive integer."))
        length(graph.edges) > 0 || throw(ArgumentError("Graph `graph` must have at least one edge."))
        new([γ], κ, graph)
    end
end

@memoize function max_k_col_subgraph_phase_separation_hamiltonian(graph::Graph, κ::Int)
    z = matrix(Z)

    # Implementation of Eq. (17)
    # one-hot encoding: n * κ vector. Index (a-1)*n + (b-1) corresponds to vertex a, color b.
    H_P_enc = sum(
        kron((color == a && vertex ∈ edge ? z : I(2) for vertex ∈ 1:graph.n for color ∈ 1:κ)...) # Z_{u,a} Z_{v,a}
        for a ∈ 1:κ for edge ∈ graph.edges # Σ_{(u,v) = edge ∈ E} Σ_{a=1..κ}
    )
    return Diagonal(H_P_enc)
end

function Qaintessent.matrix(g::MaxKColSubgraphPhaseSeparationGate)
    # Calculate the hamiltonian (Eq. 17)
    H_P_enc = max_k_col_subgraph_phase_separation_hamiltonian(g.graph, g.κ)

    # Implementation of one-hot phase seperator, A.3.1, p. 34
    U_P = exp(-im * g.γ[] * H_P_enc)

    U_P
end

Qaintessent.adjoint(g::MaxKColSubgraphPhaseSeparationGate) = MaxKColSubgraphPhaseSeparationGate(-g.γ[], g.κ, g.graph)

Qaintessent.sparse_matrix(g::MaxKColSubgraphPhaseSeparationGate) = sparse(matrix(g))

# Number of wires (κ * n in the one-hot encoding)
Qaintessent.num_wires(g::MaxKColSubgraphPhaseSeparationGate)::Int = g.κ * g.graph.n


text:
# Query incidence matrix in a database table
using PyPlot

# Create a starting set of vertices
Nv0 = 100
v0 = ceil.(Int,10000*rand(Nv0,1))

# Convert to string list.
v0str = "Out|" .* string.(v0)

# Get degrees of vertices.
Edeg = str2num(TedgeDeg[v0str,:])

# Select vertices in an out degree range.
degMin = 5;  degMax = 10;
degMax = max(degMax,minimum(getval(Edeg[:,"Degree,"]))+1)
v1str = getrow((Edeg[:,"Degree,"] > degMin) < degMax )

# Get vertex neighbors.
E = logical(Tedge[getrow(Tedge[:,v1str]),:])

figure()
spy(E)
xlabel("in/out vertex")
ylabel("edge")


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

function isotropic(coef,strain)
    isotropic_ = load_op_and_grad("./build/libIsotropic","isotropic")
    coef,strain = convert_to_tensor([coef,strain], [Float64,Float64])
    isotropic_(coef,strain)
end

# TODO: specify your input parameters
N = 10
strain = rand(N,3)
coef = rand(N, 3)

out = zeros(N,3)
for i = 1:N
    A = [strain[i,1] strain[i,3]/2
        strain[i,3]/2 strain[i,2]]
    S = coef[i,1]*I + coef[i,2] * A + coef[i,3] * A * A
    out[i,:] = [S[1,1];S[2,2];S[1,2]]
end
u = isotropic(coef,strain)
sess = Session(); init(sess)
@show run(sess, u)-out

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(m)
    # return sum(isotropic(coef,strain)^2)
    return sum(isotropic(coef,m)^2)
end

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant([ zeros(N,1) rand(N,1) zeros(N,1)])
# v_ = [ zeros(N,1) rand(N,1) zeros(N,1)]

m_ = constant([ zeros(N,2) rand(N,1)])
v_ = [ zeros(N,2) rand(N,1)]
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")


text:
module TravelingSalesmanExact

using JuMP, UnicodePlots, Logging, LinearAlgebra
import MathOptInterface
const MOI = MathOptInterface
export get_optimal_tour,
       plot_cities,
       simple_parse_tsp,
       set_default_optimizer!

# added in Julia 1.2
if VERSION < v"1.2.0-"
    import Base.>
    >(x) = Base.Fix2(>, x)
end

# added in Julia 1.1
if VERSION < v"1.1.0-"
    isnothing(::Any) = false
    isnothing(::Nothing) = true
end

const default_optimizer = Ref{Any}(nothing)

"""
    set_default_optimizer(O)

Sets the default optimizer. For example,

    using GLPK
    set_default_optimizer(GLPK.Optimizer)
"""
set_default_optimizer!(O) = default_optimizer[] = O

"""
    get_default_optimizer()

Gets the default optimizer, which is set by `set_default_optimizer`.
"""
get_default_optimizer() = default_optimizer[]

reset_default_optimizer!() = default_optimizer[] = nothing

"""
    plot_cities(cities)

Uses `UnicodePlots`'s `lineplot` to make a plot of the tour of the cities in
`cities`, in order (including going from the last city back to the first).
"""
function plot_cities(cities)
    n = length(cities)
    inc(a) = a == n ? one(a) : a + 1
    lineplot([cities[inc(j)][1] for j = 0:n], [cities[inc(j)][2] for j = 0:n])
end

"""
    find_cycle(perm_matrix, starting_ind)

Returns the cycle in the permutation described by `perm_matrix` which includes
`starting_ind`.
"""
function find_cycle(perm_matrix, starting_ind = 1)
    cycle = [starting_ind]
    prev_ind = ind = starting_ind
    while true
        # the comparisons `x > (0.5)` should mean `x == 1`. Due to floating point results returned
        # by the solvers, instead we sometimes have `x ≈ 1.0` instead. Since these are binary
        # values, we might as well just compare to 1/2.
        next_ind = findfirst(>(0.5), @views(perm_matrix[ind, 1:prev_ind-1]))
        another_ind = findfirst(>(0.5), @views(perm_matrix[ind, prev_ind+1:end]))
        if isnothing(next_ind) && isnothing(another_ind)
            cycle = Int[]
            break
        end
        if isnothing(next_ind)
            next_ind = another_ind + prev_ind
        end
        next_ind == starting_ind && break
        push!(cycle, next_ind)
        prev_ind, ind = ind, next_ind
    end
    cycle
end

"""
    get_cycles(perm_matrix)

Returns a list of cycles from the permutation described by `perm_matrix`.
"""
function get_cycles(perm_matrix)
    N = size(perm_matrix, 1)
    remaining_inds = Set(1:N)
    cycles = Vector{Int}[]
    while length(remaining_inds) > 0
        cycle = find_cycle(perm_matrix, first(remaining_inds))
        if length(cycle) > 0
            push!(cycles, cycle)
            setdiff!(remaining_inds, cycle)
        else
            setdiff!(remaining_inds, [first(remaining_inds)])
        end
    end
    cycles
end

"""
    show_tour(cities, perm_matrix)

Show a plot of the tour described by `perm_matrix` of the cities in the vector `cities`.
"""
function plot_tour(cities, perm_matrix)
    cycles = get_cycles(perm_matrix)
    tour = reduce(vcat, cycles)
    plot_cities(cities[tour])
end

"""
    remove_cycles!(model, tour_matrix)

Find the (non-maximal-length) cycles in the current solution `tour_matrix`
and add constraints to the JuMP model to disallow them. Returns the
number of cycles found.
"""
function remove_cycles!(model, tour_matrix; symmetric)
    tour_matrix_val = value.(tour_matrix)
    cycles = get_cycles(tour_matrix_val)
    length(cycles) == 1 && return 1
    for cycle in cycles
        constr = symmetric ? 2 * length(cycle) - 2 : length(cycle) - 1
        @constraint(model, sum(tour_matrix[cycle, cycle]) <= constr)
    end
    return length(cycles)
end

"""
    euclidean_distance(city1, city2)

The usual Euclidean distance measure.
"""
euclidean_distance(city1, city2) = sqrt((city1[1] - city2[1])^2 + (city1[2] - city2[2])^2)

"""
    ATT(city1, city2)

The `ATT` distance measure as specified in TSPLIB:
<https://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/tsp95.pdf>.
"""
function ATT(city1, city2)
    xd = city1[1] - city2[1]
    yd = city1[2] - city2[2]
    r = sqrt((xd^2 + yd^2) / 10.0)
    t = round(Int, r)
    if t < r
        d = t + 1
    else
        d = t
    end
    return d
end

"""
    get_optimal_tour(
        cities::AbstractVector,
        optimizer = get_default_optimizer();
        verbose = false,
        distance = euclidean_distance,
        symmetric = true,
        lazy_constraints = false,
    )

Solves the travelling salesman problem for a list of cities using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path. Optionally specify a distance metric.

The second argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an optimizer, e.g.

    get_optimal_tour(cities, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the distance metric used is symmetric (the default is to assume that it is)
* `lazy_constraints` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    cities::AbstractVector,
    optimizer = get_default_optimizer();
    verbose = false,
    distance = euclidean_distance,
    symmetric = true,
    lazy_constraints = false,
)
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    N = length(cities)
    cost = [distance(cities[i], cities[j]) for i = 1:N, j = 1:N]
    return _get_optimal_tour(cost, optimizer, symmetric, verbose, lazy_constraints, cities)

end

"""
    get_optimal_tour(
        cost::AbstractMatrix,
        optimizer = get_default_optimizer();
        verbose::Bool = false,
        symmetric::Bool = issymmetric(cost),
        lazy::Bool = true,
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

The second argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an
optimizer, e.g.

        get_optimal_tour(cities, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the `cost` matrix is symmetric (the default is to check via `issymmetric`)
* `lazy` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    cost::AbstractMatrix,
    optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(cost),
    lazy_constraints = false,
)
    size(cost, 1) == size(cost, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(cost, optimizer, symmetric, verbose, lazy_constraints)
end

"""
    get_optimal_tour(
        cost::AbstractMatrix,
        optimizer = get_default_optimizer();
        verbose::Bool = false,
        symmetric::Bool = issymmetric(cost),
        lazy::Bool = true,
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

Has additional constraits to ensure that locations are visited inside time windows.
Time windows are given as a 2d array with shape N×2.

The fourth argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an
optimizer, e.g.

        get_optimal_tour(cities, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the `cost` matrix is symmetric (the default is to check via `issymmetric`)
* `lazy` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    cost::AbstractMatrix,
    time_matrix::AbstractMatrix,
    time_windows::AbstractMatrix,
    optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(cost),
    lazy_constraints = false,
)
    size(cost, 1) == size(cost, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(cost, optimizer, false, verbose, lazy_constraints, nothing, time_matrix, time_windows)
end

"""
    get_optimal_tour(
        cost::AbstractMatrix,
        time_matrix::AbstractMatrix,
        time_windows::AbstractMatrix,
        drop_penalties::Dict,
        optimizer = get_default_optimizer();
        verbose = false,
        symmetric = issymmetric(cost),
        lazy_constraints = false,
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

The fifth argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an
optimizer, e.g.

        get_optimal_tour(cities, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the `cost` matrix is symmetric (the default is to check via `issymmetric`)
* `lazy` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    cost::AbstractMatrix,
    time_matrix::AbstractMatrix,
    time_windows::AbstractMatrix,
    drop_penalties::Dict,
    optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(cost),
    lazy_constraints = false,
)
    size(cost, 1) == size(cost, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(cost, optimizer, false, verbose, lazy_constraints, nothing, time_matrix, time_windows, drop_penalties)
end

"""
    get_optimal_tour(
        time_matrix::AbstractMatrix,
        time_windows::AbstractMatrix,
        optimizer = get_default_optimizer();
        verbose = false,
        symmetric = issymmetric(time_matrix),
        lazy_constraints = false,
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

The third argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an
optimizer, e.g.

        get_optimal_tour(cities, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the `cost` matrix is symmetric (the default is to check via `issymmetric`)
* `lazy` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    time_matrix::AbstractMatrix,
    time_windows::AbstractMatrix,
    optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(time_matrix),
    lazy_constraints = false,
)
    size(time_matrix, 1) == size(time_matrix, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(time_matrix, optimizer, false, verbose, lazy_constraints, nothing, time_matrix, time_windows)
end

"""
    get_optimal_tour(
        time_matrix::AbstractMatrix,
        time_windows::AbstractMatrix,
        drop_penalties::Dict,
        optimizer = get_default_optimizer();
        verbose = false,
        symmetric = issymmetric(time_matrix),
        lazy_constraints = false,
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

The fourth argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an
optimizer, e.g.

        get_optimal_tour(times, time_windows, penalties, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the `cost` matrix is symmetric (the default is to check via `issymmetric`)
* `lazy` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    time_matrix::AbstractMatrix,
    time_windows::AbstractMatrix,
    drop_penalties::Dict,
    optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(time_matrix),
    lazy_constraints = false,
)
    size(time_matrix, 1) == size(time_matrix, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(time_matrix, optimizer, false, verbose, lazy_constraints, nothing, time_matrix, time_windows, drop_penalties)
end

"""
    get_optimal_tour(
        cost::AbstractMatrix,
        drop_penalties::Dict,
        optimizer = get_default_optimizer();
        verbose = false,
        symmetric = issymmetric(cost),
        lazy_constraints = false,
    )

Solves the travelling salesman problem for a square cost matrix using JuMP by
formulating a MILP using the Dantzig-Fulkerson-Johnson formulation and
adaptively adding constraints to disallow non-maximal cycles. Returns an optimal
tour and the cost of the optimal path.

The third argument is mandatory if a default optimizer has not been set (via
`set_default_optimizer`). This argument should be a function which creates an
optimizer, e.g.

        get_optimal_tour(cost, penalties, GLPK.Optimizer)

There are three boolean optional keyword arguments:

* `verbose` indicates whether or not to print lots of information as the algorithm proceeds.
* `symmetric` indicates whether or not the `cost` matrix is symmetric (the default is to check via `issymmetric`)
* `lazy` indicates whether lazy constraints should be used (which requires a [compatible solver](https://www.juliaopt.org/JuMP.jl/v0.21/callbacks/#Available-solvers-1)).

"""
function get_optimal_tour(
    cost::AbstractMatrix,
    drop_penalties::Dict,
    optimizer = get_default_optimizer();
    verbose = false,
    symmetric = issymmetric(cost),
    lazy_constraints = false,
)
    size(cost, 1) == size(cost, 2) || throw(ArgumentError("First argument must be a square matrix"))
    isnothing(optimizer) && throw(ArgumentError("An optimizer is required if a default optimizer has not been set."))
    return _get_optimal_tour(cost, optimizer, symmetric, verbose, lazy_constraints, nothing, nothing, nothing, drop_penalties)
end

function build_tour_matrix(model, cost::AbstractMatrix, symmetric::Bool)
    N = size(cost, 1)
    if symmetric
    # `tour_matrix` has tour_matrix[i,j] = 1 iff cities i and j should be connected
        @variable(model, tour_matrix[1:N, 1:N], Symmetric, binary = true)

    # cost of the tour
        @objective(model, Min, sum(tour_matrix[i, j] * cost[i, j] for i = 1:N, j = 1:i))
        for i = 1:N
            @constraint(model, sum(tour_matrix[i, :]) == 2) # degree of each city is 2
            @constraint(model, tour_matrix[i, i] == 0) # rule out cycles of length 1
        end
    else
        # `tour_matrix` will be a permutation matrix
        @variable(model, tour_matrix[1:N, 1:N], binary = true)
        @objective(model, Min, sum(tour_matrix[i, j] * cost[i, j] for i = 1:N, j = 1:N))
        for i = 1:N
            @constraint(model, sum(tour_matrix[i, :]) == 1) # row-sum is 1
            @constraint(model, sum(tour_matrix[:, i]) == 1) # col-sum is 1
            @constraint(model, tour_matrix[i, i] == 0) # rule out cycles of length 1
            for j = 1:N
                @constraint(model, tour_matrix[i, j] + tour_matrix[j, i] <= 1) # rule out cycles of length 2
            end
        end
    end
    return tour_matrix
end

function build_tour_matrix(model, cost::AbstractMatrix, drop_penalties::Dict, symmetric::Bool)
    N = size(cost, 1)
    @variable(model, visit[1:N], binary = true)
    penalty_array = zeros(N)
    for i ∈ 1:N
        if i ∈ keys(drop_penalties)
            penalty_array[i] = drop_penalties[i]
        else
            @constraint(model, visit[i] == 1)
        end
    end

    if symmetric
        # `tour_matrix` has tour_matrix[i,j] = 1 iff cities i and j should be connected
        @variable(model, tour_matrix[1:N, 1:N], Symmetric, binary = true)

       # cost of the tour
        @objective(model, Min, sum(tour_matrix[i, j] * cost[i, j] for i = 1:N, j = 1:i) + sum((1-visit[i])*penalty_array[i] for i = 1:N))
        for i = 1:N
            @constraint(model, sum(tour_matrix[i, :]) == 2*visit[i]) # degree of each city is 2
            @constraint(model, tour_matrix[i, i] == 0) # rule out cycles of length 1
        end
    else
       # `tour_matrix` will be a permutation matrix
        @variable(model, tour_matrix[1:N, 1:N], binary = true)
        @objective(model, Min, sum(tour_matrix[i, j] * cost[i, j] for i = 1:N, j = 1:N) + sum((1-visit[i])*penalty_array[i] for i = 1:N))
        for i = 1:N
            @constraint(model, sum(tour_matrix[i, :]) == visit[i]) # row-sum is 1
            @constraint(model, sum(tour_matrix[:, i]) == visit[i]) # col-sum is 1
            @constraint(model, tour_matrix[i, i] == 0) # rule out cycles of length 1
            for j = 1:N
                @constraint(model, tour_matrix[i, j] + tour_matrix[j, i] <= 1) # rule out cycles of length 2
            end
        end
    end
    return tour_matrix
end

function _get_optimal_tour(
    cost::AbstractMatrix,
    optimizer,
    symmetric,
    verbose,
    lazy_constraints,
    cities = nothing,
    times = nothing,
    time_windows = nothing,
    drop_penalties = Dict{Int,Number}(),
)
    has_cities = !isnothing(cities)

    model = Model(optimizer)
    if length(drop_penalties) > 0
        tour_matrix = build_tour_matrix(model, cost, drop_penalties, symmetric)
    else
        tour_matrix = build_tour_matrix(model, cost, symmetric)
    end
    if !isnothing(times) && !isnothing(time_windows)
        @variable(model, time_windows[i,1] <= time_vars[i = 1:size(time_windows, 1)] <= time_windows[i,2])
        @variable(model, time_windows[1,1] <= end_time <= time_windows[1,2])
        for i in 1:size(cost, 1)
            for j in 2:size(cost, 1)
                @constraint(model, time_vars[i] + times[i,j] <= time_vars[j] + 10*maximum(time_windows)*(1-tour_matrix[i,j]))
            end
            @constraint(model, time_vars[i] + times[i,1] <= end_time + 10*maximum(time_windows)*(1-tour_matrix[i,1]))
        end
    end

    if has_cities && verbose
        @info "Starting optimization." plot_cities(cities)
    elseif verbose
        @info "Starting optimization."
    end

    # counts for logging
    iter = Ref(0)
    tot_cycles = Ref(0)

    if lazy_constraints
        remove_cycles_callback = make_remove_cycles_callback(model, tour_matrix, has_cities, cities, verbose, symmetric, tot_cycles)
        MOI.set(model, MOI.LazyConstraintCallback(), remove_cycles_callback)
    end

    num_cycles = 2 # just something > 1

    while num_cycles > 1
        t = @elapsed optimize!(model)
        status = termination_status(model)
        status == MOI.OPTIMAL || @warn("Problem status not optimal; got status $status")
        has_values(model) || return Int[], Inf
        num_cycles = remove_cycles!(model, tour_matrix; symmetric = symmetric)
        tot_cycles[] += num_cycles
        iter[] += 1
        if verbose
            if num_cycles == 1
                description = "found a full cycle!"
            else
                description = "disallowed $num_cycles cycles."
            end

            if has_cities
                @info "Iteration $(iter[]) took $(round(t, digits = 3))s, $description" plot_tour(
                    cities,
                    value.(tour_matrix),
                )
            else
                @info "Iteration $(iter[]) took $(round(t, digits = 3))s, $description"
            end
        end
    end
    tot_cycles[] -= 1 # remove the true cycle

    status = termination_status(model)
    status == MOI.OPTIMAL || @warn(status)

    cycles = get_cycles(value.(tour_matrix))
    length(cycles) == 1 || error("Something went wrong; did not elimate all subtours. Please file an issue.")

    if verbose
        @info "Optimization finished; adaptively disallowed $(tot_cycles[]) cycles."
        @info "Final path has length $(objective_value(model))."
        @info "Final problem has $(num_constraints(model, VariableRef, MOI.ZeroOne)) binary variables,
            $(num_constraints(model,
            GenericAffExpr{Float64,VariableRef}, MOI.LessThan{Float64})) inequality constraints, and
            $(num_constraints(model, GenericAffExpr{Float64,VariableRef}, MOI.EqualTo{Float64})) equality constraints."
    end
    return first(cycles), objective_value(model)
end

function make_remove_cycles_callback(model, tour_matrix, has_cities, cities, verbose, symmetric, tot_cycles)
    num_triggers = Ref(0)
    return function remove_cycles_callback(cb_data)
        tour_matrix_val = callback_value.(Ref(cb_data), tour_matrix)

        # We only handle integral solutions
        # Could possibly also find cycles in the integer part of a mixed solution.
        any(x -> !(x ≈ round(Int, x)), tour_matrix_val) && return

        num_triggers[] += 1
        cycles = get_cycles(tour_matrix_val)

        if length(cycles) == 1
            if has_cities && verbose
                @info "Lazy constaint triggered ($(num_triggers[])); found a full cycle!" plot_tour(
                    cities,
                    tour_matrix_val,
                )
            elseif verbose
                @info "Lazy constaint triggered ($(num_triggers[])); found a full cycle!"
            end
            return nothing
        end

        for cycle in cycles
            constr = symmetric ? 2 * length(cycle) - 2 : length(cycle) - 1
            cycle_constraint = @build_constraint( sum(tour_matrix[cycle, cycle]) <= constr)
            MOI.submit(model, MOI.LazyConstraint(cb_data), cycle_constraint)
        end

        num_cycles = length(cycles)
        tot_cycles[] += num_cycles
        if has_cities && verbose
            @info "Lazy constaint triggered ($(num_triggers[])); disallowed $num_cycles cycles." plot_tour(
                cities,
                tour_matrix_val,
            )
        elseif verbose
            @info "Lazy constaint triggered ($(num_triggers[])); disallowed $num_cycles cycles."
        end
    end
end

"""
    simple_parse_tsp(filename; verbose = true)

Try to parse the ".tsp" file given by `filename`. Very simple implementation
just to be able to test the optimization; may break on other files. Returns a
list of cities for use in `get_optimal_tour`.
"""
function simple_parse_tsp(filename; verbose = true)
    cities = Vector{Int}[]
    for line in readlines(filename)
        if startswith(line, '1':'9')
            nums = split(line, " ")
            @assert length(nums) == 3
            x = parse(Int, nums[2])
            y = parse(Int, nums[3])
            push!(cities, [x, y])
        elseif verbose
            println(line)
        end
    end
    return cities
end

"""
    get_ATT48_cities() -> Vector{Vector{Int}}

A simple helper function to get the problem data for the ATT48 TSPLIB problem.

# Example

```julia
using TravelingSalesmanExact, GLPK
cities = TravelingSalesmanExact.get_ATT48_cities()
get_optimal_tour(cities, GLPK.Optimizer, distance = TravelingSalesmanExact.ATT)
```
"""
function get_ATT48_cities()
    path = joinpath(@__DIR__, "..", "data", "att48.tsp")
    cities = simple_parse_tsp(path; verbose = false)
    return cities
end

end # module


text:
using ModelSanitizer
using Statistics
using Test

mutable struct LinearModel{T}
    X::Matrix{T}
    y::Vector{T}
    beta::Vector{T}
    function LinearModel{T}()::LinearModel{T} where T
        m::LinearModel{T} = new()
        return m
    end
end

function fit!(m::LinearModel{T}, X::Matrix{T}, y::Vector{T})::LinearModel{T} where T
    m.X = deepcopy(X)
    m.y = deepcopy(y)
    m.beta = beta = (m.X'm.X)\(m.X'm.y)
    return m
end

function predict(m::LinearModel{T}, X::Matrix{T})::Vector{T} where T
    y_hat::Vector{T} = X * m.beta
    return y_hat
end

function predict(m::LinearModel{T})::Vector{T} where T
    X::Matrix{T} = m.X
    y_hat::Vector{T} = predict(m, X)
    return y_hat
end

function mse(y::Vector{T}, y_hat::Vector{T})::T where T
    _mse::T = mean((y .- y_hat).^2)
    return _mse
end

function mse(m::LinearModel{T}, X::Matrix{T}, y::Vector{T})::T where T
    y_hat::Vector{T} = predict(m, X)
    _mse::T = mse(y, y_hat)
    return _mse
end

function mse(m::LinearModel{T})::T where T
    X::Matrix{T} = m.X
    y::Vector{T} = m.y
    _mse::T = mse(m, X, y)
    return _mse
end

rmse(varargs...) = sqrt(mse(varargs...))

function r2(y::Vector{T}, y_hat::Vector{T})::T where T
    y_bar::T = mean(y)
    SS_tot::T = sum((y .- y_bar).^2)
    SS_res::T = sum((y .- y_hat).^2)
    _r2::T = 1 - SS_res/SS_tot
    return _r2
end

function r2(m::LinearModel{T}, X::Matrix{T}, y::Vector{T})::T where T
    y_hat::Vector{T} = predict(m, X)
    _r2::T = r2(y, y_hat)
    return _r2
end

function r2(m::LinearModel{T})::T where T
    X::Matrix{T} = m.X
    y::Vector{T} = m.y
    _r2::T = r2(m, X, y)
    return _r2
end

X = randn(Float64, 200, 14)
y = X * randn(Float64, 14) + randn(200)
m = LinearModel{Float64}()
testing_rows = 1:2:200
training_rows = setdiff(1:200, testing_rows)
fit!(m, X[training_rows, :], y[training_rows])

@test m.X == X[training_rows, :]
@test m.y == y[training_rows]
@test all(m.X .== X[training_rows, :])
@test all(m.y .== y[training_rows])
@test !all(m.X .== 0)
@test !all(m.y .== 0)

# before sanitization, we can make predictions
predict(m, X[testing_rows, :])
predict(m, X[training_rows, :])
@show mse(m, X[training_rows, :], y[training_rows])
@show rmse(m, X[training_rows, :], y[training_rows])
@show r2(m, X[training_rows, :], y[training_rows])
@show mse(m, X[testing_rows, :], y[testing_rows])
@show rmse(m, X[testing_rows, :], y[testing_rows])
@show r2(m, X[testing_rows, :], y[testing_rows])

sanitize!(Model(m), Data(X), Data(y)) # sanitize the model with ModelSanitizer

@test m.X != X[training_rows, :]
@test m.y != y[training_rows]
@test !all(m.X .== X[training_rows, :])
@test !all(m.y .== y[training_rows])
@test all(m.X .== 0)
@test all(m.y .== 0)

# after sanitization, we are still able to make predictions
predict(m, X[testing_rows, :])
predict(m, X[training_rows, :])
@show mse(m, X[training_rows, :], y[training_rows])
@show rmse(m, X[training_rows, :], y[training_rows])
@show r2(m, X[training_rows, :], y[training_rows])
@show mse(m, X[testing_rows, :], y[testing_rows])
@show rmse(m, X[testing_rows, :], y[testing_rows])
@show r2(m, X[testing_rows, :], y[testing_rows])

# if you know exactly where the data are stored inside the model, you can
# directly delete them with ForceSanitize:
sanitize!(ForceSanitize(m.X), ForceSanitize(m.y))

# we can still make predictions even after using ForceSanitize
predict(m, X[testing_rows, :])
predict(m, X[training_rows, :])
@show mse(m, X[training_rows, :], y[training_rows])
@show rmse(m, X[training_rows, :], y[training_rows])
@show r2(m, X[training_rows, :], y[training_rows])
@show mse(m, X[testing_rows, :], y[testing_rows])
@show rmse(m, X[testing_rows, :], y[testing_rows])
@show r2(m, X[testing_rows, :], y[testing_rows])


text:
using LinearAlgebra

begin const
    n = 22
    m = 16
    mr = 3
    md = m - mr

    orig = [ 1 2 3 4 5 6 7 8 8 9 10 11 13 1 2 4 5 7 8 14 2 10];
    dest = [ 4 16 15 5 6 10 16 9 12 10 11 14 15 16 6 8 9 11 13 15 4 13];

    absn = [11 18 38 4 8 15 26 4 10 19 26 7 21 33 33 16];
    ordn = [28 21 8 21 17 17 26 9 13 13 18 4 9 18 12 24];

    r = Float64[100, 10,1000, 100, 100, 10, 1000, 100, 1000, 100,
                1000, 1000, 1000, 10, 10, 100, 100, 1000, 100, 1000, 100, 10];
    pr = Float64[105, 104, 110];
    fd = Float64[0.08, -1.30, 0.13, 0.09, 0.16, 0.14, 0.12, 0.07, 0.17, 0.11,
                0.25, 0.01, 0.13];

    ################################################################################
    A = zeros(m, n);

    for l = 1:n
        A[orig[l], l] = +1
        A[dest[l], l] = -1
    end

    Ar = A[1:mr, :]
    Ad = A[mr+1:m, :];

    AdT = Ad[:, 1:md]
    AdI = inv(AdT);

    AdC = Ad[:, md+1:n];

    # Here I stands for identity matrix (former eye function)
    B = [-AdI * AdC ; I];

    # Vecteur des débits admissibles
    q0 = [AdI * fd; zeros(n-md)];
end

function hydrauliqueP(qc)
    # débit des arcs
    q = q0 + B*qc
    # pertes de charge des arcs
    z = r .* abs.(q) .* q
    # flux aux noeuds
    f = [Ar*q ; fd]
    # pressions aux noeuds
    temp = (Ar'*pr) + z
    p = [pr; -AdI'*temp[1:md]]
    return q, z, f, p
end

function hydrauliqueD(pd)
    # pressions aux noeuds
    p = [pr ; pd]
    # pertes de charge aux arcs
    z = -A' * p;
    # débit des arcs
    q = z ./ sqrt(r .* abs.(z))
    # flux aux noeuds
    f = [Ar*q; fd]
    return q, z, f, p
end

# vérification des équations d'équilibre d'un réseau de
# distribution d'eau
function verification(q, z, f, p)
    # écarts maximaux sur les lois de Kirchhoff
    tol1 = max(abs.(A*q-f))
    tol2 = max(abs.(A'*p+z))
    println("Vérification des équations d'équilibre du réseau")
    println("sur les débits: ", tol1, "\t sur les pressions: ", tol2)
end

max_stars_repo_name: frapac/Knitro-Zoo

max_issues_repo_name: frapac/Knitro-Zoo

max_forks_repo_name: frapac/Knitro-Zoo


text:
module UnivariateProblems

    ### Sources
    ###
    ### [1]  http://infinity77.net/global_optimization/test_functions_1d.html

    struct UnivariateProblem
        name::AbstractString
        f::Function
        bounds::Vector{Float64}
        minimizers::Vector{Float64}
        minima::Vector{Float64}
    end

    examples = Dict{AbstractString, UnivariateProblem}()

    f(x) = 2x^2+3x+1

    examples["Polynomial"] = UnivariateProblem("Polynomial",
                                                  f,
                                                  [-2.0, 1.0],
                                                  [-0.75,],
                                                  [f(-0.75),])

    # Problem 04 from [1]
    p04(x) = -(16x^2-24x+5)exp(-x)

    examples["Problem04"] = UnivariateProblem("Problem04",
                                                p04,
                                                [1.9, 3.9],
                                                [2.868034,],
                                                [p04(2.868034),])

    # Problem 13 from [1]
    p13(x) = -x^(2/3)-(1-x^2)^(1/3)

    examples["Problem13"] = UnivariateProblem("Problem13",
                                                p13,
                                                [0.001, 0.99],
                                                [1 ./ sqrt(2.0),],
                                                [p13(1 ./ sqrt(2.0)),])

    # Problem 18 from [1]
    p18(x) = x <= 3.0 ? (x-2.0)^2 : 2.0*log(x-2.0)+1.0

    examples["Problem18"] = UnivariateProblem("Problem18",
                                                p18,
                                                [0.0, 6.0],
                                                [2.0,],
                                                [p18(2.0),])
end


text:
using GraphTools
using Base.Test

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

@testset "GraphTools" begin
  @testset "Undirected Graph" begin
    include("undirected.jl")
  end
  @testset "Directed Graph" begin
    include("directed.jl")
  end
  @testset "GraphStore" begin
    include("graphstore.jl")
  end
end


text:
"""
`Layout` has following public fields:

* `root` - the root folder of the data
* `subjects` - a vector of `Subject`
* `longitudinal` - `true` if a longitudinal study, i.e. multi-session/multi-visit per
  subject.
* `description` - dictionary parsed from `dataset_description.json`.
* `subjects_details` - `DataFrame` from `subjects.tsv`

The BIDS directory must follow the specification, i.e.
`root/sub-<subject_id>/[ses-<session_id>]`

The `Layout` can be initialized by specifying only `root`. Other behavior can also be
tweaked accordingly with the optional parameters:

* `search::Bool` - defaults to `true`. If `true` this will search over `root` to find
  `Subject`, `Session` and `Files`.
* `load_metadata::Bool` - defaults to `true`. If `true`, this expects a JSON-sidecar to
  be present alongside with every BIDS files.
* `require_modality::Bool` - defaults to `true`. If `true` this expects a modality
  presents in every file name, e.g. sub-subtest_ses-1_run-001_**T1w**.nii.gz
* `longitudinal::Bool` - defaults to `true`. If `true` this expects a session directory
  `ses-` exists under subject directory.
* `strict::Bool` - defaults to `true`. If `true`, `BIDSTools` will throw errors on
  invalid BIDS filenames, this can be turned to `false` to not parse those file names
  and display a warning instead.
* `extract_from_full_path::Bool` - defaults to `true`. If `true`, will try to extract
  `sub` and `ses` from full path and append to `entities` if they can't be found in the
  parsed filename.

# Example

```julia-repl
julia> layout = Layout("/path/to/bids/root/")
Layout:
    root = "/path/to/bids/root/"
    total subject = 49
    total session = 74
    total files = 1161
```
"""
struct Layout
    root::String
    subjects::Vector{Subject}
    # With this turned off, not going to look for ses- in path
    longitudinal::Bool
    description::OrderedDict{String, Any}
    subjects_detail::DataFrame
end

function Layout(
    root::AbstractString;
    search::Bool=true,
    load_metadata::Bool=true,
    require_modality::Bool=true,
    longitudinal::Bool=true,
    strict::Bool=true,
    extract_from_full_path::Bool=true
)
    subjects = Vector{Subject}()
    if search
        # Search for subject, session, and files
        for d in readdir(root)
            if isdir(joinpath(root, d)) && startswith(d, "sub-")
                push!(
                    subjects,
                    Subject(
                        joinpath(root, d),
                        search=search,
                        load_metadata=load_metadata,
                        require_modality=require_modality,
                        longitudinal=longitudinal,
                        strict=strict,
                        extract_from_full_path=extract_from_full_path
                    )
                )
            end
        end
    end
    description = !isfile(joinpath(root, "dataset_description.json")) ?
                      OrderedDict{String,Any}() :
                      JSON.parsefile(
                          joinpath(root, "dataset_description.json"),
                          dicttype=OrderedDict{String,Any}
                      )
    subjects_detail = !isfile(joinpath(root, "subjects.tsv")) ?
                              DataFrame() :
                              CSV.File(joinpath(root, "subjects.tsv"), delim="\t")|>
                                  DataFrame!
    Layout(root, subjects, longitudinal, description, subjects_detail)
end

function Layout(root::AbstractString, longitudinal::Bool)
    Layout(root, [], true)
end

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

"""
    function total_subjects(layout::Layout)

Get number of subject in the layout
"""
function total_subjects(layout::Layout)
    length(layout.subjects)
end

# Docstring in Subject.jl
function total_sessions(layout::Layout)
    sum(total_sessions, layout.subjects)
end

# Docstring in Session.jl
function total_files(layout::Layout)
    sum(total_files, layout.subjects)
end

function Base.show(io::IO, layout::Layout)
    print(io, """
        Layout:
            root = $(layout.root)
            total subject = $(total_subjects(layout))
            total session = $(total_sessions(layout))
            total files = $(total_files(layout))""")
end

"""
    function print_dataset_description(layout::Layout)

Function to pretty-print dataset description
"""
function print_dataset_description(layout::Layout)
    print(json(layout.description, 2))
end

"""
    function list_subject_detail(layout::Layout)

Function to pretty print subject spreadsheet (subjects.tsv) using PrettyTables
"""
function list_subject_detail(layout::Layout)
    pretty_table(layout.subjects_detail, crop= :none)
end

# Docstring of get_files is in File.jl
function get_files(
    layout::Layout; path::Union{String, Regex, Nothing}=nothing, kws...
)
    result = Vector{File}()
    for sub in layout.subjects
        files = get_files(sub; path=path, kws...)
        push!(result, files...)
    end
    return result
end