text:
using Stheno

# Dependencies that are also used in Stheno that are not in the standard library.
# This is a hack to prevent CompatHelper from trying to maintain a version number for these
# packages in both Stheno and test.
using Stheno.AbstractGPs
using Stheno.BlockArrays
using Stheno.KernelFunctions
using Stheno.Zygote

# Dependencies that are test-specific.
using Documenter
using FiniteDifferences
using LinearAlgebra
using Random
using Statistics
using Test
using TimerOutputs

using Stheno:
    mean,
    cov,
    var,
    GPC,
    AV,
    FiniteGP,
    block_diagonal,
    AbstractGP,
    BlockData,
    blocks,
    cross,
    ColVecs,
    Xt_invA_Y,
    Xt_invA_X,
    diag_At_A,
    diag_At_B,
    diag_Xt_invA_X,
    diag_Xt_invA_Y,
    block_diagonal,
    BlockDiagonal,
    blocksizes

using Stheno.AbstractGPs.Distributions: MvNormal
using FiniteDifferences: j′vp

const to = TimerOutput()

macro timedtestset(name, code)
    return esc(:(@timeit to $name @testset $name $code))
end

include("test_util.jl")

@testset "Stheno" begin

    println("util:")
    @timedtestset "util" begin
        include(joinpath("util", "zygote_rules.jl"))
        include(joinpath("util", "covariance_matrices.jl"))
        @testset "block_arrays" begin
            include(joinpath("util", "block_arrays", "test_util.jl"))
            include(joinpath("util", "block_arrays", "dense.jl"))
            include(joinpath("util", "block_arrays", "diagonal.jl"))
        end
        include(joinpath("util", "abstract_data_set.jl"))
    end

    println("gp:")
    @timedtestset "gp" begin
        include(joinpath("gp", "gp.jl"))
    end

    println("composite:")
    @timedtestset "composite" begin
        include(joinpath("composite", "test_util.jl"))
        include(joinpath("composite", "cross.jl"))
        include(joinpath("composite", "product.jl"))
        include(joinpath("composite", "addition.jl"))
        include(joinpath("composite", "compose.jl"))
    end

    println("abstract_gp:")
    @timedtestset "abstract_gp" begin
        include("abstract_gp.jl")
    end

    println("gaussian_process_probabilistic_programme:")
    include("gaussian_process_probabilistic_programme.jl")

    println("sparse_finite_gp:")
    include("sparse_finite_gp.jl")

    println("doctests")
    @timedtestset "doctests" begin
        DocMeta.setdocmeta!(
            Stheno,
            :DocTestSetup,
            :(using Stheno.AbstractGPs, Stheno, Random, Documenter, LinearAlgebra);
            recursive=true,
        )
        doctest(Stheno)
    end
end

display(to)


text:
"""
    networks(dataset::MangalDataset, query::Pair...)

Returns networks that are part of a `MangalDataset`. Allows additional query parameters.
"""
function networks(dataset::MangalDataset, query::Pair...)
    return networks(Pair("dataset_id", dataset.id), query...)
end

"""
    network(name::AbstractString)

Returns a network of a given name.
"""
function network(name::AbstractString)
    q = networks(Pair("name", name))
    return isequal(1)(length(q)) ? only(q) : nothing
end

ext: jl


text:
include("some_functions.jl")
include("some_variables.jl")

function main()
    say_name(n)
    say_age(a)
end

main()


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct SalientPoleMachine <: Machine
        R::Float64
        Td0_p::Float64
        Td0_pp::Float64
        Tq0_pp::Float64
        Xd::Float64
        Xq::Float64
        Xd_p::Float64
        Xd_pp::Float64
        Xl::Float64
        Se::Tuple{Float64, Float64}
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of 3-states salient-pole synchronous machine with quadratic/exponential saturation:
IEEE Std 1110 §5.3.1 (Model 2.1). GENSAL or GENSAE model in PSSE and PSLF.

# Arguments
- `R::Float64`: Armature resistance, validation range: (0, nothing)
- `Td0_p::Float64`: Time constant of transient d-axis voltage, validation range: (0, nothing)
- `Td0_pp::Float64`: Time constant of sub-transient d-axis voltage, validation range: (0, nothing)
- `Tq0_pp::Float64`: Time constant of sub-transient q-axis voltage, validation range: (0, nothing)
- `Xd::Float64`: Reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xq::Float64`: Reactance after EMF in q-axis per unit, validation range: (0, nothing)
- `Xd_p::Float64`: Transient reactance after EMF in d-axis per unit, validation range: (0, nothing)
- `Xd_pp::Float64`: Sub-Transient reactance after EMF in d-axis per unit. Note: Xd_pp = Xq_pp, validation range: (0, nothing)
- `Xl::Float64`: Stator leakage reactance, validation range: (0, nothing)
- `Se::Tuple{Float64, Float64}`: Saturation factor at 1 and 1.2 pu flux: Se(eq_p) = B(eq_p-A)^2
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	eq_p: q-axis generator voltage behind the transient reactance,
	ψ_kd: flux linkage in the first equivalent damping circuit in the d-axis,
	ψq_pp: phasonf of the subtransient flux linkage in the q-axis
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct SalientPoleMachine <: Machine
    "Armature resistance"
    R::Float64
    "Time constant of transient d-axis voltage"
    Td0_p::Float64
    "Time constant of sub-transient d-axis voltage"
    Td0_pp::Float64
    "Time constant of sub-transient q-axis voltage"
    Tq0_pp::Float64
    "Reactance after EMF in d-axis per unit"
    Xd::Float64
    "Reactance after EMF in q-axis per unit"
    Xq::Float64
    "Transient reactance after EMF in d-axis per unit"
    Xd_p::Float64
    "Sub-Transient reactance after EMF in d-axis per unit. Note: Xd_pp = Xq_pp"
    Xd_pp::Float64
    "Stator leakage reactance"
    Xl::Float64
    "Saturation factor at 1 and 1.2 pu flux: Se(eq_p) = B(eq_p-A)^2"
    Se::Tuple{Float64, Float64}
    ext::Dict{String, Any}
    "The states are:
	eq_p: q-axis generator voltage behind the transient reactance,
	ψ_kd: flux linkage in the first equivalent damping circuit in the d-axis,
	ψq_pp: phasonf of the subtransient flux linkage in the q-axis"
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function SalientPoleMachine(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se, ext=Dict{String, Any}(), )
    SalientPoleMachine(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se, ext, [:eq_p, :ψ_kd, :ψq_pp], 3, InfrastructureSystemsInternal(), )
end

function SalientPoleMachine(; R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se, ext=Dict{String, Any}(), )
    SalientPoleMachine(R, Td0_p, Td0_pp, Tq0_pp, Xd, Xq, Xd_p, Xd_pp, Xl, Se, ext, )
end

# Constructor for demo purposes; non-functional.
function SalientPoleMachine(::Nothing)
    SalientPoleMachine(;
        R=0,
        Td0_p=0,
        Td0_pp=0,
        Tq0_pp=0,
        Xd=0,
        Xq=0,
        Xd_p=0,
        Xd_pp=0,
        Xl=0,
        Se=(0.0, 0.0),
        ext=Dict{String, Any}(),
    )
end

"""Get SalientPoleMachine R."""
get_R(value::SalientPoleMachine) = value.R
"""Get SalientPoleMachine Td0_p."""
get_Td0_p(value::SalientPoleMachine) = value.Td0_p
"""Get SalientPoleMachine Td0_pp."""
get_Td0_pp(value::SalientPoleMachine) = value.Td0_pp
"""Get SalientPoleMachine Tq0_pp."""
get_Tq0_pp(value::SalientPoleMachine) = value.Tq0_pp
"""Get SalientPoleMachine Xd."""
get_Xd(value::SalientPoleMachine) = value.Xd
"""Get SalientPoleMachine Xq."""
get_Xq(value::SalientPoleMachine) = value.Xq
"""Get SalientPoleMachine Xd_p."""
get_Xd_p(value::SalientPoleMachine) = value.Xd_p
"""Get SalientPoleMachine Xd_pp."""
get_Xd_pp(value::SalientPoleMachine) = value.Xd_pp
"""Get SalientPoleMachine Xl."""
get_Xl(value::SalientPoleMachine) = value.Xl
"""Get SalientPoleMachine Se."""
get_Se(value::SalientPoleMachine) = value.Se
"""Get SalientPoleMachine ext."""
get_ext(value::SalientPoleMachine) = value.ext
"""Get SalientPoleMachine states."""
get_states(value::SalientPoleMachine) = value.states
"""Get SalientPoleMachine n_states."""
get_n_states(value::SalientPoleMachine) = value.n_states
"""Get SalientPoleMachine internal."""
get_internal(value::SalientPoleMachine) = value.internal

"""Set SalientPoleMachine R."""
set_R!(value::SalientPoleMachine, val::Float64) = value.R = val
"""Set SalientPoleMachine Td0_p."""
set_Td0_p!(value::SalientPoleMachine, val::Float64) = value.Td0_p = val
"""Set SalientPoleMachine Td0_pp."""
set_Td0_pp!(value::SalientPoleMachine, val::Float64) = value.Td0_pp = val
"""Set SalientPoleMachine Tq0_pp."""
set_Tq0_pp!(value::SalientPoleMachine, val::Float64) = value.Tq0_pp = val
"""Set SalientPoleMachine Xd."""
set_Xd!(value::SalientPoleMachine, val::Float64) = value.Xd = val
"""Set SalientPoleMachine Xq."""
set_Xq!(value::SalientPoleMachine, val::Float64) = value.Xq = val
"""Set SalientPoleMachine Xd_p."""
set_Xd_p!(value::SalientPoleMachine, val::Float64) = value.Xd_p = val
"""Set SalientPoleMachine Xd_pp."""
set_Xd_pp!(value::SalientPoleMachine, val::Float64) = value.Xd_pp = val
"""Set SalientPoleMachine Xl."""
set_Xl!(value::SalientPoleMachine, val::Float64) = value.Xl = val
"""Set SalientPoleMachine Se."""
set_Se!(value::SalientPoleMachine, val::Tuple{Float64, Float64}) = value.Se = val
"""Set SalientPoleMachine ext."""
set_ext!(value::SalientPoleMachine, val::Dict{String, Any}) = value.ext = val
"""Set SalientPoleMachine states."""
set_states!(value::SalientPoleMachine, val::Vector{Symbol}) = value.states = val
"""Set SalientPoleMachine n_states."""
set_n_states!(value::SalientPoleMachine, val::Int64) = value.n_states = val
"""Set SalientPoleMachine internal."""
set_internal!(value::SalientPoleMachine, val::InfrastructureSystemsInternal) = value.internal = val


text:
using Test, Tody

include("test.jl")

using Tody: Ex, Eval, define, step!, done, run, lookup, @pseudo_str

# Shortcut
macro p_str(str)
    :(@pseudo_str($str))
end

@testset TodyTestSet "Tody" begin

    # The second round takes basically no time, but makes sure things are sort
    # of idempotent
    for i = 1:2

    @testset TodyTestSet "Orders" begin

        using Tody: Order

        P = Order(Int)

        @test cmp(1, 1) == 0
        @test_throws ErrorException cmp(1, 2, P)

        push!(P, 1 => 2)

        @test cmp(1, 2, P) == -1
        @test cmp(2, 1, P) == +1

        push!(P, 3 => 4)

        @test_throws ErrorException cmp(1, 4, P)

        push!(P, 2 => 3)

        @test cmp(1, 4, P) == -1
        @test_throws ErrorException push!(P, 4 => 2)

    end

    @testset TodyTestSet "Parsing" begin

        using Tody: Ex

        ## Parsing the syntax

        # Corner cases
        @test p" " isa Ex
        @test p"1" isa Ex
        @test p" { } " isa Ex
        @test p"if true {} else {}" isa Ex

        @test string(p"break; return") == "(do (break) (return))"

        @test string(p"-x^2") == "(- (^ x 2))"

        @test_throws ErrorException parse("1 / 2 * 3")

        @test string(p"1 + -2 - +3") == "(+ 1 (- (- 2) (+ 3)))"

        # Reluctant search, from *Pessimal Algorithms and Simplexity Analysis* by
        # Broder and Stolfi.

        ex = pseudo"""

        procedure research(X, i, j) {
            if i > j {
                return -1
            }
            if i == j && X == A {
                return i
            } else {
                return -1
            }
            m = (i + j) / 2
            if X ≤ A[m] {
                k = research(X, m+1, j)
                if k == -1 {
                    return research(X, i, m)
                } else {
                    return k
                }
            } else {
                k = research(X, i, m)
                if k == -1 {
                    return research(X, m+1, j)
                } else {
                    return k
                }
            }
        }

        """

        expected = """
        (procedure
           (research X i j)
           (do (if (> i j) (return (- 1)))
               (if (&& (== i j) (== X A))
                   (return i)
                   else
                   (return (- 1)))
               (= m (/ (+ i j) 2))
               (if (≤ X (at A m))
                   (do (= k (research X (+ m 1) j))
                       (if (== k (- 1))
                           (return (research X i m))
                           else
                           (return k)))
                   else
                   (do (= k (research X i m))
                       (if (== k (- 1))
                           (return (research X (+ m 1) j))
                           else
                           (return k))))))"""

        @test string(ex) == expected

        ## Parsing S-expressions

        @test ex == parse(Ex, expected)

    end

    @testset TodyTestSet "Partial eval" begin

        using Tody: Ex, peek, Eval, done

        input = "(+ (* 2 2) (- 3))"

        ex = parse(Ex, input)
        ev = Eval(ex)

        @test !done(ev)

        @test string(ev) == string(ex) == input

        push!(ev, :+)
        @test string(ev) == input

        @test string(peek(ev)) == "(* 2 2)"
        push!(ev, 4)
        @test string(ev) == "(+ 4 (- 3))"

        @test string(peek(ev)) == "(- 3)"
        push!(ev, -3)
        @test string(ev) == "(+ 4 -3)"

        @test done(ev)

    end

    output(ex::String; trace=false) = output(parse(ex), trace=trace)
    output(ex::Ex; trace=false) = output(run(ex), trace=trace)
    function output(ev::Eval; trace=false)
        trace && println()
        cur = ev
        while cur isa Ex
            trace && println(cur)
            cur = step!(cur)
        end
        String(take!(Tody.lookup(ev, :stdout)))
    end

    @testset TodyTestSet "Short-circuit" begin

        @test output(p"1 + 1 == 2 && print(1)") == "1"
        @test output(p"2 + 2 == 5 && print(1)") == ""
        @test output(p"1 + 1 == 2 || print(1)") == ""
        @test output(p"2 + 2 == 5 || print(1)") == "1"

    end

    @testset TodyTestSet "if/else if/else" begin

        ex = pseudo"""

            if true {
                print(1)
            } else {
                print(0)
            }

            if false {
                print(0)
            } else {
                print(1)
            }

            if false {
                print(0)
            } else if true {
                print(1)
            } else {
                print(0)
            }

            if false {
                print(0)
            } else if false {
                print(0)
            } else {
                print(1)
            }

        """

        @test output(ex) == "1111"

    end

    @testset TodyTestSet "while" begin

        ex = pseudo"""

            a = 93
            b = 42

            while b ≠ 0 {
                tmp = b
                b = a % b
                a = tmp
            }

            print(a)

        """

        @test output(ex) == "3"

    end

    @testset TodyTestSet "break" begin

        ex = p"""

            while true {
                print(1)
                break
                print(0)
            }

        """

        @test output(ex) == "1"

    end

    @testset TodyTestSet "continue" begin

        ex = p"""

            i = 1
            while i ≤ 3 {
                print(i)
                i = i + 1
                if i == 2 {
                    continue
                }
                print(i)
            }

        """

        @test output(ex) == "12334"

    end

    @testset TodyTestSet "for" begin

        ex = p"""

            for i in 0:9 {
                if i == 2 {
                    continue
                }
                print(i)
                if i == 8 {
                    break
                }
            }

        """

        @test output(ex) == "01345678"

        ex = p"""

            for x in A {
                print(x)
            }

        """

        ev = run(ex)

        define(ev, :A, [1, 2, 3])

        @test output(ev) == "123"

    end

    @testset TodyTestSet "definition" begin

        ex = p"""

            i = 1
            j = 1

            while true {
                i := 2
                j = 2
                break
            }

            print(i, j)

        """

        @test output(ex) == "12"

    end

    @testset TodyTestSet "Proc" begin

        using Tody: Proc, Env, BUILTINS

        env = BUILTINS()

        mysum = Proc(:mysum, [:x, :y], p"x + y", Env(env))

        define(env, :mysum, mysum)

        ex = p"print(mysum(2^1, 1+1))"
        ev = Eval(ex, env=env)

        @test output(ev) == "4"

    end

    @testset TodyTestSet "ProcDef" begin

        ex = p"""

        procedure mypow(x, y) {

            res = 1
            i = 0

            while i < y {
                res = res * x
                i = i + 1
            }

            return res

        }

        print(mypow(2, 3))

        """

        @test output(ex) == "8"

    end

    @testset TodyTestSet "return" begin

        ex = p"""

        procedure f() {
            if true {
                return
            }
            print(0)
        }

        print(1)
        f()
        print(2)

        """
        @test output(ex) == "12"

    end

    @testset TodyTestSet "recursion" begin

        ex = p"""

        procedure myfact(n) {
            if n == 0 {
                1
            } else {
                n * myfact(n - 1)
            }
        }

        print(myfact(5))

        """

        @test output(ex) == "120"

    end

    @testset TodyTestSet "closures" begin

        ex = p"""

        procedure outer() {

            var = "inside"

            procedure inner() {
                print(var)
            }

            return inner

        }

        var = "outside"

        proc = outer()

        proc()

        """

        @test output(ex) == "inside"

    end

    @testset TodyTestSet "attributes" begin

        ex = p"""

            print(e.h.π.d)
            print(e.h.d)

        """

        ev = run(ex)

        s = (d=0, π=nothing)
        v = (d=1, π=s)
        e = (t=s, h=v)

        define(ev, :e, e)

        @test output(ev) == "01"

        ex = p"""

            print(a.b.head)
            a.b.head = y
            print(a.b.head)

        """

        ev = run(ex)

        a = (b=Expr(:x),)
        define(ev, :a, a)
        define(ev, :y, :y)
        define(ev, :z, :z)

        @test output(ev) == "xy"

    end

    @testset TodyTestSet "indexing" begin

        ex = p"""

            print(A[3, 2])
            A[1, 1] = 4
            print(A[1, 1])

        """

        ev = run(ex)

        define(ev, :A, collect(reshape(1:16, 4, 4)))

        @test output(ev) == "74"

    end

    end

end


text:
@doc raw"""
ABC equilibrium in (x,y,z) coordinates with covariant components of the vector
potential given by
```math
A (x,y,z) = \big( a \, \sin(z) + c \, \cos(y) , \, b \, \sin(x) + a \, \cos(z) , \, c \, \sin(y) + b \, \cos(x) \big)^T
```
resulting in the magnetic field ``B(x,y,z) = A(x,y,z)``.

Parameters: `a`, `b`, `c`
"""
struct ABC{T <: Number} <: AnalyticEquilibrium
    name::String
    a::T
    b::T
    c::T

    ABC{T}(a::T, b::T, c::T) where T <: Number = new("ABCEquilibrium", a, b, c)
end

ABC(a::T, b::T, c::T) where T <: Number = ABC{T}(a, b, c)

function Base.show(io::IO, equ::ABC)
    print(io, "ABC Equilibrium with\n")
    print(io, "  A = ", equ.a, "\n")
    print(io, "  B = ", equ.b, "\n")
    print(io, "  C = ", equ.c)
end

function X(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    x[1]
end

function Y(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    x[2]
end

function Z(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    x[3]
end

function J(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    one(T)
end

function A₁(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    equ.a * sin(x[3]) + equ.c * cos(x[2])
end

function A₂(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    equ.b * sin(x[1]) + equ.a * cos(x[3])
end

function A₃(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    equ.c * sin(x[2]) + equ.b * cos(x[1])
end

function g₁₁(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    one(T)
end

function g₂₂(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    one(T)
end

function g₃₃(x::AbstractArray{T,1}, equ::ABC) where {T <: Number}
    one(T)
end

macro abc_equilibrium(a, b, c)
    generate_equilibrium_code(ABC(a, b, c); output=false)
end


text:
# Timings from Strong scaling 256 core runs

using PyPlot

# 1 thread, 8 thread, 32 thread
t1 = [6.62, 1.3, 4.0, 0.24, 0.21, 0.69] ./ 5
t2 = [6.85, 1.29, 4.01, 0.24, 0.2, 0.81] ./ 5
t3 = [10.69, 1.34, 4.05, 0.27, 0.88, 3.87] ./ 5

t1_percent = [99.34, 19.6, 60.4, 3.7, 3.1, 10.5]
t2_percent = [99.11, 18.9, 58.5, 3.5, 3.0, 11.8]
t3_percent = [99.53, 12.5, 37.9, 2.5, 8.2, 36.18]

#=
t_all = [6.62, 6.85, 10.69]
t_facevalue = [1.3, 1.29, 1.34]
t_faceflux = [4.0, 4.01, 4.05]
t_source = [0.07, 0.07, 0.08]
t_update = [0.24, 0.24, 0.27]
t_updatecheck = [0.21, 0.2, 0.88]
t_messagepass = [0.69, 0.81, 3.87]
=#

# set width of bar
const barWidth = 0.3

# Set position of bar on X axis
r1 = collect(0:length(t1)-1)
r2 = [x + barWidth for x in r1]
r3 = [x + barWidth for x in r2]

fig = figure("bar_timing",figsize=(14,6))
ax = plt.axes()
#axis("tight")

# Make the plot
a1 = bar(r1, t1, width=barWidth, edgecolor="white", label="1 thread * 256 MPI")
a2 = bar(r2, t2, width=barWidth, edgecolor="white", label="8 threads * 32 MPI")
a3 = bar(r3, t3, width=barWidth, edgecolor="white", label="32 threads * 8 MPI")

# Add xticks on the middle of the group bars
xlabel("module", fontweight=:bold, fontsize=18)
xticks([r + barWidth for r in 0:length(t1)-1],
   ["Advance", "FaceValue", "FaceFlux", "Update", "UpdateCheck", "MessagePass"],
   fontsize=18)
plt.minorticks_on()
ax.tick_params(axis=:x, which=:minor, bottom=false)
ax.yaxis.grid(true, linestyle="--")
#ax.yaxis.grid(true, which="minor")
ylabel("time [s]", fontsize=18)
ax.set_title("Timings per time step for 256 cores, explicit 2nd Order Linde + mc3 limiter",fontsize=18)

# Create legend & Show graphic
legend(fontsize=18)

function autolabel(rects, t_percent)
   """Attach a text label above each bar in *rects*, displaying its height."""
   for (i, rect) in enumerate(rects)
      i == 1 && continue
      height = rect.get_height()
      percent = t_percent[i]
      ax.annotate(string(percent),
         xy=(rect.get_x() + rect.get_width() / 2, height),
         xytext=(0, 3),  # 3 points vertical offset
         textcoords="offset points",
         ha=:center, va=:bottom)
   end
end

autolabel(a1, t1_percent)
autolabel(a2, t2_percent)
autolabel(a3, t3_percent)


text:
import Compat.LinearAlgebra: RealHermSymComplexHerm, normalize
import Base: reverse

convert(::Type{BitSet},r::IntCol) = BitSet(r)
BitSet(elems::Vararg{Int}) = BitSet(elems)

reverse(n::Number) = n

tname(T::Type) = T.name.name

function randomsmooth(n::Integer,m::Integer)
    R = Matrix{Float64}(undef,n,m)
    for i = 1:m
        R[:,i] = randomsmooth(n)
    end
    return R
end

function randomsmooth(n::Integer)
    x = transpose(0:1/(n-1):1)
    N = rand(4:7) # number of sines to use
    f = 2.0 .+ 10 .* rand(N) # frequencies
    A = normalize!((rand(N).-0.5)./f,2) # amplitudes
    ϕ = 2π.*rand(N) # phase
    return vec(sum(A.*sin.(f.*x.+ϕ),1))
end

gaussian(x::Real,σ::Real) = exp(-0.5*(x/σ)^2)
gaussianprime(x::Real,σ::Real) = -exp(-0.5*(x/σ)^2)*x/σ^2

function gaussian_fwhm(x::Real,w)
    # w is the full width half maximum (FWHM)
    return exp(-4ln2*(x/w)^2)
end

function gaussianprime_fwhm(x::Real,w)
    # w is the full width half maximum (FWHM)
    return -x*ln2*exp2(3-4(x/w)^2)/w^2
end

function sqrtfact(n::Integer)
    if 0 <= n < 128
        return sqrtfact_table[n+1]
    elseif n >= 128 # Stirling series for values not stored in the table
        return sqrt(sqrt(2π*n)*(n/ℯ)^n*(1.0+1/12n+1/288n^2))
    else
        throw(ArgumentError("n must be larger than or equal to 0"))
    end
end

# pairwise sum of squared differences (SSD)
function norm2_diff(A::Vector{T},B::Vector{T}) where {T}
    n = length(A)
    n == length(B) || throw(DimensionMismatch())
    n == 0 ? abs2(zero(T)-zero(T)) : norm2_diff(A, B, 1, n)
end

function norm2_diff(A::Vector, B::Vector, i1::Integer, n::Integer)
    if n < 128
        @inbounds s = abs2(A[i1] - B[i1])
        for i in i1+1 : i1+n-1
            @inbounds s += abs2(A[i] - B[i])
        end
        return s
    else
        n2 = n÷2
        return norm2_diff(A, B, i1, n2) + norm2_diff(A, B, i1+n2, n-n2)
    end
end

normalize(z::Complex) = z == 0 ? one(z) : z/abs(z)
normalize(z::Real) = one(z)

function unwrap!(p)
    n,m = size(p,1),size(p,2)
    n < 2 && return p
    for j = 1:m, i = 2:n
        d = p[i,j] - p[i-1,j]
        if abs(d) > π
            p[i,j] -= floor((d+π)/2π) * 2π
        end
    end
    return p
end

expim(H::RealHermSymComplexHerm) = expim!(Matrix{complex(eltype(H))}(H),copy(H))
expim!(H::RealHermSymComplexHerm) = expim!(Matrix{complex(eltype(H))}(H),H)

function expim!(R::Matrix,H::RealHermSymComplexHerm,Λ=Vector{real(eltype(R))}(undef,size(H,1)),U=Matrix(H),B=similar(R))
    # First decompose H into U*Λ*Uᴴ
    #F = eigfact!(H); copy!(Λ,F.values); copy!(U,F.vectors)
    hermfact!(Λ,U,H)
    # Calculate the imaginary exponential of each eigenvalue and multiply
    # each column of U to obtain B = U*exp(iΛ)
    n = length(Λ)
    @inbounds for j = 1:n
        a = cis(Λ[j])
        @simd for i = 1:n
            B[i,j] = a * U[i,j]
        end
    end
    # Finally multiply B by Uᴴ to obtain U*exp(iΛ)*Uᴴ = exp(i*H)
    @static if VERSION < v"0.7.0-"
        A_mul_Bc!(R,B,U)
    else
        mul!(R,B,adjoint(U))
    end
    return R
end

excited_sys(A::QuObject, subs::Vararg{Integer}) = excited_sys(dims(A), sys)
excited_sys(::Union{Dims{N},Val{N}}, subs::Vararg{Integer}) where N =
    ntuple(i->ifelse(i ∈ subs, 1, 0), Val(N))

lengthrange(a::Real,len::Integer) = UnitRange{typeof(a)}(a, oftype(a, a+len-1))

tensored_iterator(dims::Dims) = (revtuple(t) for t ∈ product(lengthrange.(0,revtuple(dims))...))


text:
@testset "zygote_friendly_map" begin
    @testset "$name" for (name, f, x) in [
        ("Vector{Float64}", x -> sin(x) + cos(x) * exp(x), randn(100)),
        ("Fill{Float64}", x -> sin(x) + exp(x) + 5, Fill(randn(), 100)),
        ("Vector{Vector{Float64}}", sum, [randn(25) for _ in 1:33]),
        (
            "zip(Vector{Float64}, Fill{Float64})",
            x -> x[1] * x[2],
            zip(randn(5), Fill(1.0, 5)),
        ),
    ]
        @test TemporalGPs.zygote_friendly_map(f, x) ≈ map(f, x)
        adjoint_test(x -> TemporalGPs.zygote_friendly_map(f, x), (x, ))
    end
end


text:
using DrWatson
@quickactivate "PECUZAL_Julia"

# These tests are only for development and further merging the functionality into
# DelayEmbeddings.jl and are not strictly needed for this scientific project

using DelayEmbeddings
using DynamicalSystemsBase
using Test
using DelimitedFiles
using Random

include("../src/pecuzal_method.jl")
include("../src/data_analysis_functions.jl")

println("\nTesting pecuzal_method.jl for comparison to MATLAB and Python implementations:")
@testset "PECUZAL other algorithm comparison" begin
    s = readdlm("./test/timeseries/lorenz_pecora_uni_x.csv")
    s = vec(s[1:500]) # input timeseries = x component of lorenz
    Tmax = 100
    theiler = 21

    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding_update(s;
                               τs = 0:Tmax , w = theiler)
    @test size(Y,2) == 4
    @test -0.6134 < sum(Ls) < -0.6133
    @test τ_vals[2]  == 21
    @test τ_vals[3]  == 13
    @test τ_vals[4]  == 78

    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding_update(s;
                               τs = 0:Tmax , w = theiler, econ = true)

    @test -0.6078 < sum(Ls) < -0.6077
    @test τ_vals[2]  == 21
    @test τ_vals[3]  == 13
    @test τ_vals[4]  == 78

    tr = readdlm("./test/timeseries/lorenz_pecora_multi.csv")
    tr = Dataset(tr) # input timeseries = x component of lorenz
    w = 15
    Tmax = 100

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding_update(tr[1:500,1:2];
                                         τs = 0:Tmax , w = w)
    @test size(Y,2) == 2
    @test ts_vals[1] == 2
    @test ts_vals[2] == 1
    @test τ_vals[1] == τ_vals[2] == 0
    @test sum(Ls) < -0.5505736

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding_update(tr[1:500,1:2];
                                         τs = 0:Tmax , w = w, econ = true)
    @test size(Y,2) == 2
    @test ts_vals[1] == 2
    @test ts_vals[2] == 1
    @test τ_vals[1] == τ_vals[2] == 0
    @test sum(Ls) < -0.544942
end

println("\nTesting pecuzal_method.jl...")
@testset "PECUZAL" begin
@testset "Univariate example" begin
    # For comparison reasons using Travis CI we carry out the integration on a UNIX
    # OS and load the resulting time series (see /test/timeseries/produce_timeseries.jl)

    s = readdlm("./test/timeseries/lorenz_pecora_uni_x.csv")
    s = vec(s) # input timeseries = x component of lorenz
    w = estimate_delay(s, "mi_min")
    Tmax = 100
    K = 14
    samplesize = 1
    KNN = 3

    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding_update(s[1:5000];
                                        τs = 0:Tmax , w = w, samplesize = samplesize,
                                        K = K, KNN = KNN)
    L1 = sum(Ls)
    L2 = compute_delta_L(s[1:5000], τ_vals, Tmax;
         KNN = KNN, w = w)
    @test L1 == L2
    @test -0.728 < Ls[1] < -0.727
    @test -0.3234 < Ls[2] < -0.3233

    @test τ_vals[2] == 18
    @test τ_vals[3] == 9

    @test length(ts_vals) == 3

    @time Y, τ_vals, ts_vals, Ls , εs = pecuzal_embedding_update(s[1:5000];
                                        τs = 0:Tmax , w = w, samplesize = samplesize,
                                        K = K, KNN = KNN, econ = true)

    L1 = sum(Ls)

    @test L1 < -1.047465
    @test -0.7263 < Ls[1] < -0.7261
    @test -0.3213 < Ls[2] < -0.3212

    @test τ_vals[2] == 18
    @test τ_vals[3] == 9

    @test length(ts_vals) == 3

    YY1 = DelayEmbeddings.hcat_lagged_values(Dataset(s), s, 21)
    YY2 = DelayEmbeddings.hcat_lagged_values(YY1, s, 13)
    @test length(YY1) == length(YY2)
    @test YY1 == YY2[:,1:2]

end

@testset "Multivariate example" begin
# Test of the proposed Pecora-Uzal-embedding-method (multivariate case)

    # For comparison reasons using Travis CI we carry out the integration on a UNIX
    # OS and load the resulting time series (see /test/timeseries/produce_timeseries.jl)

    include("../src/pecuzal_method.jl")
    include("../src/data_analysis_functions.jl")
    tr = readdlm("./test/timeseries/lorenz_pecora_multi.csv")
    tr = Dataset(tr) # input timeseries = x component of lorenz
    w1 = estimate_delay(tr[:,1], "mi_min")
    w2 = estimate_delay(tr[:,2], "mi_min")
    w3 = estimate_delay(tr[:,3], "mi_min")
    w = w1
    Tmax = 100
    samplesize = 1
    KNN = 3

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding_update(tr[1:5000,:];
                                         τs = 0:Tmax , w = w, samplesize = samplesize)
    L1 = sum(Ls)
    L2 = compute_delta_L(tr[1:5000,:], τ_vals, ts_vals, Tmax; KNN = KNN, w = w)
    @test L1 == L2
    @test length(ts_vals) == 5
    @test ts_vals[3] == ts_vals[4] == ts_vals[5] == 1
    @test ts_vals[1] == 3
    @test ts_vals[2] == 2
    @test τ_vals[2] == 0
    @test τ_vals[3] == 62
    @test τ_vals[4] == 48
    @test τ_vals[5] == 0
    @test -0.936 < Ls[1] < -0.935
    @test -0.356 < Ls[2] < -0.355
    @test -0.131 < Ls[3] < -0.130
    @test -0.015 < Ls[4] < -0.014

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding_update(tr[1:5000,:];
                        τs = 0:Tmax , w = w, L_threshold = 0.4, samplesize = samplesize)
    @test -0.936 < Ls[1] < -0.935
    @test length(Ls) == 1

    # Dummy input
    include("../src/pecuzal_method.jl")
    Random.seed!(1234)
    d1 = randn(1000)
    d2 = rand(1000)
    dummy_set = Dataset(hcat(d1,d2))
    w1 = estimate_delay(d1, "mi_min")
    w2 = estimate_delay(d2, "mi_min")
    w = w1

    @time Y, τ_vals, ts_vals, Ls , ε★ = pecuzal_embedding_update(dummy_set;
                                        τs = 0:Tmax , w = w, samplesize = samplesize)

    @test size(Y,2) == 1

end

end


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

__precompile__(true)

"""
Tools for distributed parallel processing.
"""
module Distributed

# imports for extension
import Base: getindex, wait, put!, take!, fetch, isready, push!, length,
             hash, ==, kill, close, isopen, showerror

# imports for use
using Base: Process, Semaphore, JLOptions, AnyDict, buffer_writes, wait_connected,
            VERSION_STRING, binding_module, notify_error, atexit, julia_exename,
            julia_cmd, AsyncGenerator, acquire, release, invokelatest,
            shell_escape_posixly, uv_error, coalesce, notnothing

using Serialization, Sockets
import Serialization: serialize, deserialize
import Sockets: connect

# NOTE: clusterserialize.jl imports additional symbols from Serialization for use

export
    @spawn,
    @spawnat,
    @fetch,
    @fetchfrom,
    @everywhere,
    @distributed,

    AbstractWorkerPool,
    addprocs,
    CachingPool,
    clear!,
    ClusterManager,
    default_worker_pool,
    init_worker,
    interrupt,
    launch,
    manage,
    myid,
    nprocs,
    nworkers,
    pmap,
    procs,
    remote,
    remotecall,
    remotecall_fetch,
    remotecall_wait,
    remote_do,
    rmprocs,
    workers,
    WorkerPool,
    RemoteChannel,
    Future,
    WorkerConfig,
    RemoteException,
    ProcessExitedException,

    process_messages,
    remoteref_id,
    channel_from_id,
    worker_id_from_socket,
    cluster_cookie,
    start_worker,

# Used only by shared arrays.
    check_same_host

function _require_callback(mod::Base.PkgId)
    if Base.toplevel_load[] && myid() == 1 && nprocs() > 1
        # broadcast top-level (e.g. from Main) import/using from node 1 (only)
        @sync for p in procs()
            p == 1 && continue
            @async remotecall_wait(p) do
                Base._require(mod)
                nothing
            end
        end
    end
end

include("clusterserialize.jl")
include("cluster.jl")   # cluster setup and management, addprocs
include("messages.jl")
include("process_messages.jl")  # process incoming messages
include("remotecall.jl")  # the remotecall* api
include("macros.jl")      # @spawn and friends
include("workerpool.jl")
include("pmap.jl")
include("managers.jl")    # LocalManager and SSHManager
include("precompile.jl")

# Deprecations

@eval @deprecate $(Symbol("@parallel")) $(Symbol("@distributed"))

# PR 26783
@deprecate pmap(p::AbstractWorkerPool, f, c; kwargs...) pmap(f, p, c; kwargs...)
@deprecate pmap(p::AbstractWorkerPool, f, c1, c...; kwargs...) pmap(f, p, c1, c...; kwargs...)

function __init__()
    push!(Base.package_callbacks, _require_callback)
    init_parallel()
end

end

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
# Code for Video 1.4: Controlling the flow

a = b = c = 4
n, m = 1, 2

area = (r = 3; π*r^2) # 28.274333882308138

x = begin
           a = 2
           b = 3
           a + b
    end

body_temp = 38.5 # 38.5
if body_temp > 38
  println("this person has fever")
elseif 35 < body_temp <= 37
  println("this person has normal temperature")
else
  println("this person has hypothermia")
end
# this person has fever

n = 10; m = 15
max = if n > m   n
      else       m
      end
# max == 15
max = n > m ? n : m  # 15

year = 2016
while year <= 2050
    println("increase year")
    year += 1
end

while year <= 2050
    println("increase year")
    year += 1
    if year == 2047
        println("reached 2047")
        break
    end
end

# while true
#     println("I do not stop!")
# end

blood_groups = ["A", "B", "AB", "O"]
# 4-element Array{ASCIIString,1}:
#  "A"
#  "B"
#  "AB"
#  "O"
for blg in blood_groups
    print("$blg / ")
end
# A / B / AB / O /
for (ix, blg) in enumerate(blood_groups)
    println("Blood-group no. $ix is $blg")
end
# Blood-group no. 1 is A
# Blood-group no. 2 is B
# Blood-group no. 3 is AB
# Blood-group no. 4 is O

for i in 1:5
    print("$i -")
end
# 1 - 2 - 3 - 4 - 5 -

# open("test.txt")
# ERROR: SystemError: opening file test.txt: No such file or directory
#  in open at iostream.jl:90
#  in open at iostream.jl:99
println("")

try
    open("test.txt")
catch ex
    showerror(STDOUT, ex)
    println("")
end
println("I recovered from the exception")

blood_groups = ["A", "B", "AB", "O"]
code = "X"
if code in blood_groups
   println("we have a supply of this blood")
else
#   throw(DomainError())
end
print("After error: not reached")
# ERROR: DomainError:
println("")

if code in blood_groups
   println("we have a supply of this blood")
else
#   error("blood group unknown!")
end
print("After error: not reached")
# ERROR: blood group unknown!
#  in error at error.jl:21
println("")

if code in blood_groups
   println("we have a supply of this blood")
else
   warn("blood group unknown!")
end
print("After error: reached!")

z = @async 26 * 78
println("In the meantime, I calculate other things!")
println(consume(z)) # 2028

max_stars_repo_path: Section 1/1_4.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_4.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_4.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
import SeisIO: demean, demean!, taper, taper!,detrend, detrend!
export detrend, detrend!, taper, taper!, demean, demean!, bandpass, bandpass!
export bandstop, bandstop!, lowpass, lowpass!, highpass, highpass!
export phase, phase!, hanningwindow
# Signal processing functions for arrays (rather than SeisData or SeisChannel)

"""
    detrend!(X::AbstractArray{<:Union{Float32,Float64},1})

Remove linear trend from array `X` using least-squares regression.
"""
function detrend!(X::AbstractArray{<:Union{Float32,Float64},1})
    T = eltype(X)
    N = length(X)
    A = ones(T,N,2)
    A[:,1] .= range(T(1/N),T(1),length=N)
    coeff = A \ X
    X[:] .= X .- A *coeff
    return nothing
end
detrend(A::AbstractArray{<:Union{Float32,Float64},1}) = (U = deepcopy(A);
        detrend!(U);return U)

"""
    detrend!(X::AbstractArray{<:Union{Float32,Float64},2})

Remove linear trend from columns of `X` using least-squares regression.
"""
function detrend!(X::AbstractArray{<:Union{Float32,Float64},2})
    M,N = size(X)
    T = eltype(X)
    A = ones(T,M,2)
    A[:,1] .= range(T(1),stop=T(M)) ./ T(M)

    # solve least-squares through qr decomposition
    Q,R = qr(A)
    rq = inv(factorize(R)) * Q'
    for ii = 1:N
        coeff = rq * X[:,ii]
        X[:,ii] .-=  A *coeff
    end
    return nothing
end
detrend(A::AbstractArray{<:Union{Float32,Float64},2}) = (U = deepcopy(A);
        detrend!(U);return U)
detrend!(R::RawData) = detrend!(R.x)
detrend(R::RawData) = (U = deepcopy(R); detrend!(U.x); return U)

"""
    demean!(A::AbstractArray{<:Union{Float32,Float64},1})

Remove mean from array `A`.
"""
function demean!(A::AbstractArray{<:Union{Float32,Float64},1})
      μ = mean(A)
      for ii = 1:length(A)
        A[ii] -= μ
      end
  return nothing
end
demean(A::AbstractArray{<:Union{Float32,Float64},1}) = (U = deepcopy(A);
       demean!(U);return U)

"""
   demean!(A::AbstractArray{<:Union{Float32,Float64},2})

Remove mean from columns of array `A`.
"""
function demean!(A::AbstractArray{<:Union{Float32,Float64},2})
      M,N = size(A)
      for ii = 1:N
        μ = mean(A[:,ii])
        for jj = 1:M
          A[jj,ii] -= μ
        end
      end
  return nothing
end
demean(A::AbstractArray{<:Union{Float32,Float64},2}) = (U = deepcopy(A);
       demean!(U);return U)
demean!(R::RawData) = demean!(R.x)
demean(R::RawData) = (U = deepcopy(R); demean!(U.x); return U)

"""
   taper!(A,fs; max_percentage=0.05, max_length=20.)

Taper a time series `A` with sampling_rate `fs`.
Defaults to 'hann' window. Uses smallest of `max_percentage` * `fs`
or `max_length`.

# Arguments
- `A::AbstractArray`: Time series.
- `fs::Float64`: Sampling rate of time series `A`.
- `max_percentage::float`: Decimal percentage of taper at one end (ranging
   from 0. to 0.5).
- `max_length::Float64`: Length of taper at one end in seconds.
"""
function taper!(A::AbstractArray{<:Union{Float32,Float64},1}, fs::Float64;
                max_percentage::Float64=0.05, max_length::Float64=20.)
   N = length(A)
   T = eltype(A)
   wlen = min(Int(floor(N * max_percentage)), Int(floor(max_length * fs)), Int(
         floor(N/2)))
   taper_sides = [-hanningwindow(T,2 * wlen -1, zerophase=true) .+ T(1);T(0)]
   A[1:wlen] .= A[1:wlen] .* taper_sides[1:wlen]
   A[end-wlen+1:end] .= A[end-wlen+1:end] .* taper_sides[wlen+1:end]
   return nothing
end
taper(A::AbstractArray{<:Union{Float32,Float64},1}, fs::Float64;
      max_percentage::Float64=0.05, max_length::Float64=20.) = (U = deepcopy(A);
      taper!(U,fs,max_percentage=max_percentage,max_length=max_length);return U)

function taper!(A::AbstractArray{<:Union{Float32,Float64},2}, fs::Float64;
                max_percentage::Float64=0.05, max_length::Float64=20.)
   M,N = size(A)
   T = eltype(A)
   wlen = min(Int(floor(M * max_percentage)), Int(floor(max_length * fs)), Int(
         floor(M/2)))
   taper_sides = [-hanningwindow(T,2 * wlen -1, zerophase=true) .+ T(1);T(0)]
   for ii = 1:N
       A[1:wlen,ii] .= A[1:wlen,ii] .* taper_sides[1:wlen]
       A[end-wlen+1:end,ii] .= A[end-wlen+1:end,ii] .* taper_sides[wlen+1:end]
   end
   return nothing
end
taper(A::AbstractArray{<:Union{Float32,Float64},2}, fs::Float64;
       max_percentage::Float64=0.05,max_length::Float64=20.) = (U = deepcopy(A);
       taper!(U,fs,max_percentage=max_percentage,max_length=max_length);return U)
taper!(R::RawData; max_percentage::Float64=0.05,
       max_length::Float64=20.) = taper!(R.x,R.fs,max_percentage=max_percentage,
       max_length=max_length)
taper(R::RawData; max_percentage::Float64=0.05,
       max_length::Float64=20.) = (U = deepcopy(R); taper!(U.x,U.fs); return U)

"""
    phase!(A::AbstractArray)

Extract instantaneous phase from signal A.

For time series `A`, its analytic representation ``S = A + H(A)``, where
``H(A)`` is the Hilbert transform of `A`. The instantaneous phase ``e^{iθ}``
of `A` is given by dividing ``S`` by its modulus: ``e^{iθ} = \\frac{S}{|S|}``
For more information on Phase Cross-Correlation, see:
[Ventosa et al., 2019](https://pubs.geoscienceworld.org/ssa/srl/article-standard/570273/towards-the-processing-of-large-data-volumes-with).
"""
function phase!(A::AbstractArray)
    A .= angle.(hilbert(A))
    return nothing
end
phase(A::AbstractArray) = (U = deepcopy(A);phase!(U);return U)
phase!(R::RawData) = phase!(R.x)
phase(R::RawData) = (U = deepcopy(R);phase!(U.x);return U)

"""
   bandpass!(A,freqmin,freqmax,fs,corners=4,zerophase=false)

Butterworth-Bandpass Filter.

Filter data `A` from `freqmin` to `freqmax` using `corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freqmin::Float64`: Pass band low corner frequency.
- `freqmax::Float64`: Pass band high corner frequency.
- `fs::Float64`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function bandpass!(A::AbstractArray{<:Union{Float32,Float64},1},
                   freqmin::Float64, freqmax::Float64, fs::Float64;
                   corners::Int=4, zerophase::Bool=false)
   fe = 0.5 * fs
   low = freqmin / fe
   high = freqmax / fe
   T = eltype(A)

   # warn if above Nyquist frequency
   if high - oneunit(high) > -1e-6
       @warn "Selected high corner frequency ($freqmax) of bandpass is at or
       above Nyquist ($fe). Applying a high-pass instead."
       highpass!(A,freqmin,fs,corners=corners,zerophase=zerophase)
       return nothing
   end

   # throw error if low above Nyquist frequency
   if low > 1
       ArgumentError("Selected low corner frequency is above Nyquist.")
   end

   # create filter
   responsetype = Bandpass(T(freqmin), T(freqmax); fs=fs)
   designmethod = Butterworth(T,corners)
   if zerophase
       A[:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:]))
   else
       A[:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:]))
   end

   return nothing
end
bandpass(A::AbstractArray{<:Union{Float32,Float64},1},freqmin::Float64,
         freqmax::Float64, fs::Float64; corners::Int=4,zerophase::Bool=false) =
         (U = deepcopy(A);bandpass!(U,freqmin,freqmax, fs, corners=corners,
         zerophase=zerophase);return U)

function bandpass!(A::AbstractArray{<:Union{Float32,Float64},2},
                   freqmin::Float64, freqmax::Float64, fs::Float64;
                   corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    low = freqmin / fe
    high = freqmax / fe
    M,N = size(A)
    T = eltype(A)

    # warn if above Nyquist frequency
    if high - oneunit(high) > -1e-6
        @warn "Selected high corner frequency ($freqmax) of bandpass is at or
        above Nyquist ($fe). Applying a high-pass instead."
        highpass!(A,freqmin,fs,corners=corners,zerophase=zerophase)
        return nothing
    end

    # throw error if low above Nyquist frequency
    if low > 1
        ArgumentError("Selected low corner frequency is above Nyquist.")
    end

    # create filter
    responsetype = Bandpass(T(freqmin), T(freqmax); fs=fs)
    designmethod = Butterworth(T,corners)
    filtx = digitalfilter(responsetype, designmethod)
    if zerophase
        for ii = 1:N
            A[:,ii] .= filtfilt(filtx, @view(A[:,ii]))
        end
    else
        for ii = 1:N
            A[:,ii] .= filt(filtx, @view(A[:,ii]))
        end
    end

    return nothing
end
bandpass(A::AbstractArray{<:Union{Float32,Float64},2},freqmin::Float64,
         freqmax::Float64, fs::Float64; corners::Int=4,zerophase::Bool=false) =
         (U = deepcopy(A);bandpass!(U,freqmin,freqmax,fs,corners=corners,
         zerophase=zerophase);return U)
bandpass!(R::RawData,freqmin::Float64,freqmax::Float64;
          corners::Int=4,zerophase::Bool=false) = (bandpass!(R.x,freqmin,freqmax,
          R.fs,corners=corners,zerophase=zerophase);setfield!(R,:freqmin,freqmin);
          setfield!(R,:freqmax,min(freqmax,R.fs/2));return nothing)
bandpass(R::RawData,freqmin::Float64,freqmax::Float64;
        corners::Int=4,zerophase::Bool=false) = (U = deepcopy(R);bandpass!(U.x,
        freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);
        setfield!(U,:freqmin,freqmin);
        setfield!(U,:freqmax,min(freqmax,R.fs/2));return U)

 """
     bandstop!(A,freqmin,freqmax,fs,corners=4,zerophase=false)

 Butterworth-Bandstop Filter.

 Filter data `A` removing data between frequencies `freqmin` to `freqmax` using
 `corners` corners.

 # Arguments
 - `A::AbstractArray`: Data to filter
 - `freqmin::Float64`: Stop band low corner frequency.
 - `freqmax::Float64`: Stop band high corner frequency.
 - `fs::Float64`: Sampling rate in Hz.
 - `fs::Int`: Filter corners / order.
 - `zerophase::Bool`: If True, apply filter once forwards and once backwards.
 This results in twice the filter order but zero phase shift in
 the resulting filtered trace.
 """
function bandstop!(A::AbstractArray{<:Union{Float32,Float64},1},
                    freqmin::Float64,freqmax::Float64,fs::Float64;
                    corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    low = freqmin / fe
    high = freqmax / fe
    T = eltype(A)

    # warn if above Nyquist frequency
    if high > 1
        @warn "Selected high corner frequency ($freqmax) is"
        "above Nyquist ($fe). Setting Nyquist as high corner."
        freqmax = fe
    end

    # throw error if low above Nyquist frequency
    if low > 1
        ArgumentError("Selected low corner frequency is above Nyquist.")
    end

    # create filter
    responsetype = Bandstop(T(freqmin), T(freqmax); fs=fs)
    designmethod = Butterworth(T,corners)
    if zerophase
        A[:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:]))
    else
        A[:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:]))
    end

    return nothing
end
bandstop(A::AbstractArray{<:Union{Float32,Float64},1},freqmin::Float64,
         freqmax::Float64, fs::Float64; corners::Int=4,zerophase::Bool=false) =
         (U = deepcopy(A);bandstop!(U,freqmin,freqmax,corners=corners,
         zerophase=zerophase);return U)

function bandstop!(A::AbstractArray{<:Union{Float32,Float64},2},
                   freqmin::Float64,freqmax::Float64,fs::Float64;
                   corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    low = freqmin / fe
    high = freqmax / fe
    M,N = size(A)
    T = eltype(A)

    # warn if above Nyquist frequency
    if high > 1
        @warn "Selected high corner frequency ($freqmax) is"
        "above Nyquist ($fe). Setting Nyquist as high corner."
        freqmax = fe
    end

    # throw error if low above Nyquist frequency
    if low > 1
        ArgumentError("Selected low corner frequency is above Nyquist.")
    end

    # create filter
    responsetype = Bandstop(T(freqmin), T(freqmax); fs=fs)
    designmethod = Butterworth(T,corners)
    filtx = digitalfilter(responsetype, designmethod)
    if zerophase
        for ii = 1:N
            A[:,ii] .= filtfilt(filtx, @view(A[:,ii]))
        end
    else
        for ii = 1:N
            A[:,ii] .= filt(filtx, @view(A[:,ii]))
        end
    end

    return nothing
end
bandstop(A::AbstractArray{<:Union{Float32,Float64},2},freqmin::Float64,
         freqmax::Float64, fs::Float64; corners::Int=4,zerophase::Bool=false) =
         (U = deepcopy(A);bandstop!(U,freqmin,freqmax,corners=corners,
         zerophase=zerophase);return U)
bandstop!(R::RawData,freqmin::Float64,freqmax::Float64;
    corners::Int=4,zerophase::Bool=false) = bandstop!(R.x,freqmin,freqmax,
    R.fs,corners=corners,zerophase=zerophase)
bandstop(R::RawData,freqmin::Float64,freqmax::Float64;
  corners::Int=4,zerophase::Bool=false) = (U = deepcopy(R);bandstop!(U.x,
  freqmin,freqmax,U.fs,corners=corners,zerophase=zerophase);return U)

"""
lowpass(A,freq,fs,corners=4,zerophase=false)

Butterworth-Lowpass Filter.

Filter data `A` over certain frequency `freq` using `corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freq::Float64`: Filter corner frequency.
- `fs::Float64`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function lowpass!(A::AbstractArray{<:Union{Float32,Float64},1},freq::Float64,
                  fs::Float64; corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    f = freq / fe
    T = eltype(A)

    # warn if above Nyquist frequency
    if f >= 1
        @warn """Selected corner frequency ($freq) is
        above Nyquist ($fe). Setting Nyquist as high corner."""
        freq = fe - 1. / fs
    end

    # create filter
    responsetype = Lowpass(T(freq); fs=fs)
    designmethod = Butterworth(T,corners)
    if zerophase
        A[:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:]))
    else
        A[:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:]))
    end
    return nothing
end
lowpass(A::AbstractArray{<:Union{Float32,Float64},1},freq::Float64, fs::Float64;
        corners::Int=4,zerophase::Bool=false) = (U = deepcopy(A);
        lowpass!(U,freq,fs,corners=corners,zerophase=zerophase);return U)

function lowpass!(A::AbstractArray{<:Union{Float32,Float64},2},freq::Float64,
                  fs::Float64; corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    f = freq / fe
    M,N = size(A)
    T = eltype(A)

    # warn if above Nyquist frequency
    if f >= 1
        @warn """Selected corner frequency ($freq) is
        above Nyquist ($fe). Setting Nyquist as high corner."""
        freq = fe - 1. / fs
    end

    # create filter
    responsetype = Lowpass(T(freq); fs=fs)
    designmethod = Butterworth(T,corners)
    filtx = digitalfilter(responsetype, designmethod)
    if zerophase
        for ii = 1:N
            A[:,ii] .= filtfilt(filtx, @view(A[:,ii]))
        end
    else
        for ii = 1:N
            A[:,ii] .= filt(filtx, @view(A[:,ii]))
        end
    end
    return nothing
end
lowpass(A::AbstractArray{<:Union{Float32,Float64},2},freq::Float64, fs::Float64;
        corners::Int=4,zerophase::Bool=false) = (U = deepcopy(A);
        lowpass!(U,freq,fs,corners=corners,zerophase=zerophase);return U)
lowpass!(R::RawData,freq::Float64; corners::Int=4,
         zerophase::Bool=false) = (lowpass!(R.x,freq,R.fs,corners=corners,
         zerophase=zerophase);setfield!(R,:freqmax,min(freqmax,R.fs/2));
         return nothing)
lowpass(R::RawData,freq::Float64; corners::Int=4,
         zerophase::Bool=false) = (U = deepcopy(R);lowpass!(U.x,freq,U.fs,
         corners=corners,zerophase=zerophase);
         setfield!(R,:freqmax,min(freqmax,R.fs/2));return U)

"""
    highpass(A,freq,fs,corners=4,zerophase=false)

Butterworth-Highpass Filter.

Filter data `A` removing data below certain frequency `freq` using `corners` corners.

# Arguments
- `A::AbstractArray`: Data to filter
- `freq::Float64`: Filter corner frequency.
- `fs::Float64`: Sampling rate in Hz.
- `fs::Int`: Filter corners / order.
- `zerophase::Bool`: If True, apply filter once forwards and once backwards.
This results in twice the filter order but zero phase shift in
the resulting filtered trace.
"""
function highpass!(A::AbstractArray{<:Union{Float32,Float64},1},freq::Float64,
                   fs::Float64; corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    f = freq / fe
    T = eltype(A)

    # warn if above Nyquist frequency
    if f > 1
        ArgumentError("Selected low corner frequency is above Nyquist.")
    end

    # create filter
    responsetype = Highpass(T(freq); fs=fs)
    designmethod = Butterworth(T,corners)
    if zerophase
        A[:] .= filtfilt(digitalfilter(responsetype, designmethod), @view(A[:]))
    else
        A[:] .= filt(digitalfilter(responsetype, designmethod), @view(A[:]))
    end
    return nothing
end
highpass(A::AbstractArray{<:Union{Float32,Float64},1},freq::Float64,fs::Float64;
         corners::Int=4,zerophase::Bool=false) = (U = deepcopy(A);
         highpass!(U,freq,fs,corners=corners,zerophase=zerophase);return U)

function highpass!(A::AbstractArray{<:Union{Float32,Float64},2},freq::Float64,
        fs::Float64; corners::Int=4, zerophase::Bool=false)
    fe = 0.5 * fs
    f = freq / fe
    M,N = size(A)
    T = eltype(A)

    # warn if above Nyquist frequency
    if f > 1
        ArgumentError("Selected low corner frequency is above Nyquist.")
    end

    # create filter
    responsetype = Highpass(T(freq); fs=fs)
    designmethod = Butterworth(T,corners)
    filtx = digitalfilter(responsetype, designmethod)
    if zerophase
        for ii = 1:N
            A[:,ii] .= filtfilt(filtx, @view(A[:,ii]))
        end
    else
        for ii = 1:N
            A[:,ii] .= filt(filtx, @view(A[:,ii]))
        end
    end
    return nothing
end
highpass(A::AbstractArray{<:Union{Float32,Float64},2},freq::Float64,fs::Float64;
corners::Int=4,zerophase::Bool=false) = (U = deepcopy(A);
highpass!(U,freq,fs,corners=corners,zerophase=zerophase);return U)
highpass!(R::RawData,freq::Float64; corners::Int=4,
         zerophase::Bool=false) = (highpass!(R.x,freq,R.fs,corners=corners,
         zerophase=zerophase);setfield!(R,:freqmin,freqmin);return nothing)
highpass(R::RawData,freq::Float64; corners::Int=4,
         zerophase::Bool=false) = (U = deepcopy(R);highpass!(U.x,freq,U.fs,
         corners=corners,zerophase=zerophase);setfield!(R,:freqmin,freqmin);
         return U)

function hanningwindow(::Type{T},n::Integer; padding::Integer=0, zerophase::Bool=false) where {T<:Real}
    if n < 0
        throw(ArgumentError("`n` must be nonnegative"))
    end
    if padding < 0
        throw(ArgumentError("`padding` must be nonnegative"))
    end
    win = zeros(T,n+padding)
    if n == 1
        win[1] = 0.5*(1+cos(2 * T(pi) *T(0.)))
    elseif zerophase
        # note that the endpoint of the window gets set in both lines. In the
        # unpadded case this will set the same index (which shouldn't make a
        # difference if the window is symmetric), but it's necessary for when
        # there's padding, which ends up in the center of the vector length
        # n÷2+1
        win[1:n÷2+1] .= 0.5 .* (1 .+ cos.(2 * T(pi) .* (range(T(0.0), stop=T(n÷2)/n, length=n÷2+1))))
        # length n÷2
        win[end-n÷2+1:end] .= 0.5 .* (1 .+ cos.(2 * T(pi) *(range(-T(n÷2)/n, stop=-1/T(n), length=n÷2))))
    else
        win[1:n] = 0.5 .*(1 .+ cos.(2 * T(pi) .* (range(-T(0.5), stop=T(0.5), length=n))))
    end
        win
end
hanningwindow(n;padding::Integer=0, zerophase::Bool=false) =
              hanningwindow(Float64,n,padding=padding,zerophase=zerophase)


text:
# Advent of Code 2020 - Day 25

test_list= split("5764801
17807724", '\n')

list = split("11239946
10464955", '\n')

list2 = split("11239946
10464955", '\n')

function loop(subject, value = 1)
    return rem(value*subject, 20201227)
end

function get_size(public_key)
    loop_size = 0
    value = 1
    while value != public_key
        loop_size += 1
        value = loop(7, value)
    end
    return loop_size
end

function transform(public_key, loop_size)
    value = 1
    for n = 1:loop_size
        value = loop(public_key, value)
    end
    return value
end

function get_encryptionkey(list)
    card_public_key, door_public_key = parse.(Int, list)
    card_loop_size = get_size(card_public_key)
    door_loop_size = get_size(door_public_key)
    min_loop = argmin([card_loop_size, door_loop_size])
    other_public = min_loop == 1 ? door_public_key : card_public_key
    loop_size = min_loop == 1 ? card_loop_size : door_loop_size
    key = transform(other_public, loop_size)
    return key
end

@show get_encryptionkey(test_list) == 14897079
@show get_encryptionkey(list) == 711945

nothing


text:
"""
Solving the clustering problem with Kmeans algorithm.
Introduce the set of (x,y) data and the nomber of K clusters in the function kmeansmy to get the groups of point that form part of every set Si.
The same function also give the centroids coordinates of each subgroup and the number of steps to achieve the agrupations.
The function graficaclusters give a plot of the data with the subgroups indicated and also display the centroids of each subgroup.
You can create a set of randn points with:
    X=vcat([(0.3*randn(),0.3*randn()) for i=1:100],[(1+.3*randn(),1+.3*randn()) for i=1:100],[(1+.3*randn(),-1+.3*randn()) for i=1:100]);
"""

function kmeansmy(data,k)
    M=[[(0.3*randn(),0.3*randn()) for i=1:k]]
    push!(M,[(0.3*randn(),0.3*randn()) for i=1:k])
    step=1;
    N=length(data)

    while M[end]!=M[end-1]
        Mnext=[]
        global S=[[] for i=1:k];

        distancias=[[] for i=1:N];

        for i=1:k,j=1:N
            push!(distancias[j],distanciacartesiana(data[j],M[end][i]))
        end

        for i=1:N,j=1:length(S)
            if minimum(distancias[i])==distancias[i][j]
                push!(S[j],data[i])
            end
        end

        for j=1:k
            push!(Mnext,(mean([S[j][i][1] for i=1:length(S[j])]),mean([S[j][i][2] for i=1:length(S[j])])))
        end

        push!(M,Mnext)
        step+=1
    end

    return S,M[end],step

end

function graficaclusters(data,k)
    S=kmeansmy(data,k)[1]
    C=kmeansmy(data,k)[2]

    scatter([S[1][j] for j=1:length(S[1])],markercolor=[:blue]);
    for i=2:k
        scatter!([S[i][j] for j=1:length(S[i])]);
    end
    scatter!([C[j] for j=1:length(C)],markercolor=[:yellow])
end

function distanciacartesiana(pto1,pto2)
    dis=sqrt((pto1[1]-pto2[1])^2+(pto1[2]-pto2[2])^2);
    return dis;
end


text:
include("eggify.jl")
using Metatheory.Library
using Metatheory.EGraphs.Schedulers

or_alg = @theory begin
    ((p ∨ q) ∨ r)       ==  (p ∨ (q ∨ r))
    (p ∨ q)             ==  (q ∨ p)
    (p ∨ p)             =>  p
    (p ∨ true)          =>  true
    (p ∨ false)         =>  p
end

and_alg = @theory begin
    ((p ∧ q) ∧ r)       ==  (p ∧ (q ∧ r))
    (p ∧ q)             ==  (q ∧ p)
    (p ∧ p)             =>  p
    (p ∧ true)          =>  p
    (p ∧ false)         =>  false
end

comb = @theory begin
    # DeMorgan
    ¬(p ∨ q)            ==  (¬p ∧ ¬q)
    ¬(p ∧ q)            ==  (¬p ∨ ¬q)
    # distrib
    (p ∧ (q ∨ r))       ==  ((p ∧ q) ∨ (p ∧ r))
    (p ∨ (q ∧ r))       ==  ((p ∨ q) ∧ (p ∨ r))
    # absorb
    (p ∧ (p ∨ q))       =>  p
    (p ∨ (p ∧ q))       =>  p
    # complement
    (p ∧ (¬p ∨ q))      =>  p ∧ q
    (p ∨ (¬p ∧ q))      =>  p ∨ q
end

negt = @theory begin
    (p ∧ ¬p)            =>  false
    (p ∨ ¬(p))          =>  true
    ¬(¬p)               ==  p
end

impl = @theory begin
    (p == ¬p)           =>  false
    (p == p)            =>  true
    (p == q)            =>  (¬p ∨ q) ∧ (¬q ∨ p)
    (p => q)            =>  (¬p ∨ q)
end

fold = @theory begin
    (true == false)     =>   false
    (false == true)     =>   false
    (true == true)      =>   true
    (false == false)    =>   true
    (true ∨ false)      =>   true
    (false ∨ true)      =>   true
    (true ∨ true)       =>   true
    (false ∨ false)     =>   false
    (true ∧ true)       =>   true
    (false ∧ true)      =>   false
    (true ∧ false)      =>   false
    (false ∧ false)     =>   false
    ¬(true)             =>   false
    ¬(false)            =>   true
end

theory = or_alg ∪ and_alg ∪ comb ∪ negt ∪ impl ∪ fold

query = :(¬(((¬p ∨ q) ∧ (¬r ∨ s)) ∧ (p ∨ r)) ∨ (q ∨ s))

###########################################

params = SaturationParams(timeout=22, eclasslimit=3051,
    scheduler=ScoredScheduler)#, schedulerparams=(1000,5, Schedulers.exprsize))

for i ∈ 1:2
    G = EGraph( query )
    report = saturate!(G, theory, params)
    ex = extract!(G, astsize)
    println( "Best found: $ex")
    println(report)
end

open("src/main.rs", "w") do f
    write(f, rust_code(theory, query, params))
end


text:
"""
    numericShift!(hh, qq, iq, qrows, qcols, neq, condn)

Compute the numeric shiftrights and store them in q.
"""
function numericShift!(hh::Array{Float64,2}, qq::Array{Float64,2}, iq::Int64, qRows::Int64, qCols::Int64, neq::Int64, condn::Float64)

    # total number of shifts
    nnumeric = 0

    # functions to seperate hh
    left = 1:qCols
    right = (qCols + 1):(qCols + neq)

    # preform QR factorization on right side of hh
    F = Compat.qr(hh[:, right], Val(true))

    # filter R only keeping rows that are zero
    Q, R = F
    zerorows = abs.(diag(R))
    zerorows = findall(x->(float(x) <= condn), zerorows)

    while (length(zerorows) != 0) && (iq <= qRows)
        # update hh with matrix multiplication of Q and hh
        hh = *(Q', hh)

        # need total number of zero rows
        nz = size(zerorows, 1)

        # update qq to keep track of rows that are zero
        qq[(iq + 1):(iq + nz), :] = hh[zerorows, left]

        # update hh by shifting right
        hh[zerorows, :] = shiftRight!(hh[zerorows, :], neq)

        # increment our variables by number of shifts
        iq = iq + nz
        nnumeric = nnumeric + nz

        # redo QR factorization and filter R as before
        F = Compat.qr(hh[:, right], Val(true))
	Q, R = F
        zerorows = abs.(diag(R))
        zerorows = findall(x->(float(x) <= condn), zerorows)

    end # while

    return(hh, qq, iq, nnumeric)

end # numericShift


text:
struct Stencil{W} end

width(::Stencil{W}) where {W} = W

# Marker and Cell
const mac = Stencil{1}()

stencil(::NTuple{N,Any}) where {N} = ntuple(i -> mac, Val(N))


text:
# Alias of Model{K,Kₑ,Kₘ, Kₑₜ,Kₘₜ, Kₑₗ,Kₘₗ, K², K₊₂,AK₊₂} for the TE waveguide mode
# mode equations for 1D cross section.
#
# The type parameter AK₊₂ specify device-specific arrays types (e.g., CuArray) and is
# user-defined in the constructor.
#
# Note that the TE mode has (Ey, Hx, Hz) as nonzero components; this is unlike TE Maxwell's
# equations that have two E-field components and one H-field component.

const ModelTE{AK₊₂} = Model{1,1,2, 1,1, 0,1, 1, 3,AK₊₂}

# Convenience constructor
function ModelTE(grid::Grid; Atype::Type=Array)
    ft_eq = EE

    cmpₛ = SInt(1)  # shapes in transverse dimension (x-axis)

    cmpₑₜ = SInt(2)  # transverse E-field (y-component)
    cmpₘₜ = SInt(1)  # transverse H-field (x-component)

    cmpₑₗ = SInt()  # longitudinal E-field (no component)
    cmpₘₗ = SInt(3)  # longitudinal H-field (z-component)

    iseₜ˔shp = true
    ishₜ˔shp = false

    return ModelTE{Atype{ComplexF,3}}(;ft_eq, grid, cmpₛ, cmpₑₜ, cmpₘₜ, cmpₑₗ, cmpₘₗ, iseₜ˔shp, ishₜ˔shp)
end

# Assign the material parameters on the grid and smooth them.
function calc_matparams!(mdl::ModelTE)
    # Take necessary fields of mdl.
    grid = mdl.grid

    isbloch = grid.isbloch
    N = grid.N
    l = grid.l
    σ = grid.σ

    lg = grid.ghosted.l
    τl = grid.ghosted.τl
    ∆τ = grid.ghosted.∆τ

    boundft = mdl.boundft
    gtₑ = ft2gt.(EE,boundft)
    gtₘ = ft2gt.(HH,boundft)

    εₜarr = mdl.εₜarr
    µₜarr = mdl.μₜarr
    µₗarr = mdl.μₗarr

    oind2shp = mdl.oind2shp
    oind2εind = mdl.oind2εind
    oind2μind = mdl.oind2μind
    εind2ε = mdl.εind2ε
    μind2μ = mdl.μind2μ

    cmpₑₜ = mdl.cmpₑₜ
    cmpₘₜ = mdl.cmpₘₜ
    cmpₘₗ = mdl.cmpₘₗ

    εind2εₜ = sub_pind2matprm(εind2ε, cmpₑₜ)
    μind2μₜ = sub_pind2matprm(μind2μ, cmpₘₜ)
    μind2μₗ = sub_pind2matprm(μind2μ, cmpₘₗ)

    iseₜ˔shp = mdl.iseₜ˔shp
    ishₜ˔shp = mdl.ishₜ˔shp
    ishₗ˔shp = true

    # Create temporary storages.
    # Because Ey and Hx are at the same locations, smoothing them can use the same object
    # index.
    pₜ_oind2d = create_oind_array(N)  # for εₜ (= εyy) and μₜ (= μxx)
    pₗ_oind2d = create_oind_array(N)  # for μₗ (= μzz)
    pₗ_oind2d′ = create_oind_array(N)  # for @assert below

    # Assign the material parameters.
    # Note that μₜarr uses gtₑ rather than gtₘ, because it is rank-0 tensor and its locations
    # are the E-field locations; see the definition MaxwellBase/assign_param!().
    # (See the figure in L10 - Eigenmode Analysis > 1.3 Waveguide mode analysis > Matrix
    # equation formulation; take the cross section through Ey, Hx, Hz.)
    assign_param!(εₜarr, tuple(pₗ_oind2d), oind2shp, oind2εind, εind2εₜ, gtₑ, τl, isbloch)  # εₜ tensors (rank-0, so scalars)
    assign_param!(µₜarr, tuple(pₗ_oind2d′), oind2shp, oind2μind, μind2μₜ, gtₑ, τl, isbloch)  # μₜ tensors (rank-0, so scalars)
    assign_param!(µₗarr, tuple(pₜ_oind2d), oind2shp, oind2μind, μind2μₗ, gtₘ, τl, isbloch)  # μₗ tensors (rank-0, so scalars)

    # Some temporary arrays should be identical.  Specifically, εyy_oind2d and μxx_oind2d
    # are evaluated at the locations surrounding the εyy and μxx locations, which are the
    # same locations.  (See the figure in L10 - Eigenmode Analysis > 1.3 Waveguide mode
    # analysis > Matrix equation formulation; take the cross section through Ey, Hx, Hz.)
    # Therefore, εyy_oind2d and μxx_oind2d are evaluated at the same locations, so they
    # should be the same object indices.
    @assert isequal(pₗ_oind2d′, pₗ_oind2d)

    # Smooth the material parameters.
    smooth_param!(εₜarr, tuple(pₜ_oind2d), oind2shp, oind2εind, εind2εₜ, gtₑ, l, lg, σ, ∆τ, iseₜ˔shp)  # εₜ tensors (rank-0, so scalars)
    smooth_param!(µₜarr, tuple(pₜ_oind2d), oind2shp, oind2μind, μind2μₜ, gtₑ, l, lg, σ, ∆τ, ishₜ˔shp)  # μₜ tensors (rank-0, so scalars)
    smooth_param!(µₗarr, tuple(pₗ_oind2d), oind2shp, oind2μind, μind2μₗ, gtₘ, l, lg, σ, ∆τ, ishₗ˔shp)  # μₗ tensors (rank-0, so scalars)

    return nothing
end

function complete_fields(β::Number,  # propagation constant
                         fₜ::AbsVecNumber,  # calculated transverse field Fₜ as column vector
                         ft::FieldType,  # type of input field Fₜ
                         ω::Number,
                         Ps::Tuple22{AbsMatNumber},
                         ∇̽s::Tuple22{AbsMatNumber},
                         πcmps::Tuple2{AbsMatNumber},
                         mdl::ModelTE)
    # Calculate the transverse components of the complementary field.
    βF′ₜgen = create_βF′ₜgen(ft, ω, Ps, ∇̽s, πcmps)
    f′ₜ = (βF′ₜgen * fₜ) ./ β

    eₜ, hₜ = ft==EE ? (fₜ,f′ₜ) : (f′ₜ,fₜ)

    # Calculate the longitudinal component of the H-field.
    iHₗgen = create_iHₗgen(ω, Ps, ∇̽s)
    hₗ = (iHₗgen * eₜ) ./ im

    # Take the Cartesian components of the fields and reshape them into arrays.
    sz_grid = size(mdl.grid)
    Ey = reshape(eₜ, sz_grid)
    Hx = reshape(hₜ, sz_grid)
    Hz = reshape(hₗ, sz_grid)

    return (E=Ey, H=(Hx, Hz))
end

# See the field configurations in L10 - Eigenmode Analysis > Matrix equation formulation.
function poynting(E::AbsVecNumber,  # Ey
                  H::Tuple2{AbsVecNumber},  # (Hx, Hz)
                  mdl::ModelTE)
    boundft = mdl.boundft
    grid = mdl.grid
    isbloch = grid.isbloch

    ∆l, ∆l⁻¹ = create_∆ls(grid, boundft)

    Ey = E
    Hx, Hz = H

    nX = 1

    Ftemp = ArrComplexF(undef, size(grid))  # vector; temporary storage

    # Calculate Sx.
    Sx = ArrComplexF(undef, size(grid))  # vector
    m̂xHz = Ftemp; interp_field!(m̂xHz, Hz, HH, nX, ∆l, ∆l⁻¹, boundft, isbloch)
    Sx .= Ey .* conj.(m̂xHz)

    # Calculate Sz.
    Sz = ArrComplexF(undef, size(grid))  # vector
    Sz₋ = Ftemp; Sz₋ .= Ey .* conj.(Hx)  # real array located at Ey
    interp_field!(Sz, Sz₋, EE, nX, ∆l, ∆l⁻¹, boundft, isbloch, α=-1.0)

    return Sx, Sz
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.subgraphview" begin

@testset "subgraphview" begin
    graph = pathgraph(5)
    subg = nodesubgraph(graph, Set([3, 4, 5]))
    @test_throws KeyError neighbors(subg, 1)[1]
    @test !hasedge(subg, 2, 3)
    @test hasedge(subg, 5, 4)

    # Node and edge key accessor should not be destructive
    nset = nodeset(subg)
    @test issetequal(nset, [3, 4, 5])
    pop!(nset)
    @test issetequal(nodeset(subg), [3, 4, 5])
    eset = edgeset(subg)
    @test issetequal(eset, [3, 4])
    pop!(eset)
    @test issetequal(edgeset(subg), [3, 4])
    @test nodecount(subg) == 3
    @test edgecount(subg) == 2

    graph = completegraph(5)
    nsubg = nodesubgraph(graph, Set([1, 2, 3, 4]))
    @test nodecount(nsubg) == 4
    @test edgecount(nsubg) == 6
    esubg = edgesubgraph(graph, Set([1, 5, 8, 10]))
    @test nodecount(esubg) == 5
    @test edgecount(esubg) == 4
    @test edgesubgraph(graph, [1, 5, 8, 10]) == esubg

    subgsubg = nodesubgraph(nsubg, Set([2, 3, 4]))
    @test issetequal(nodeset(subgsubg), [2, 3, 4])
    @test issetequal(edgeset(subgsubg), [5, 6, 8])
    @test degree(subgsubg, 4) == 2
    newg = plaingraph(subgsubg)
    @test issetequal(nodeset(newg), [1, 2, 3])
    @test issetequal(edgeset(newg), [1, 2, 3])
end

@testset "disubgraphview" begin
    graph = plaindigraph(5, [(1, 2), (2, 3), (3, 4), (4, 5)])
    subg = nodesubgraph(graph, Set([3, 4, 5]))
    @test_throws KeyError neighbors(subg, 1)[1]
    @test !hasedge(subg, 2, 3)
    @test hasedge(subg, 3, 4)
    @test !hasedge(subg, 5, 4)

    @test issetequal(nodeset(subg), [3, 4, 5])
    @test issetequal(edgeset(subg), [3, 4])
    @test nodecount(subg) == 3
    @test edgecount(subg) == 2

    graph = plaindigraph(5,
        [(1, 2), (2, 3), (3, 4), (4, 5), (5, 1),
         (1, 3), (3, 5), (5, 2), (2, 4), (4, 1)])
    nsubg = nodesubgraph(graph, Set([1, 2, 3, 4]))
    @test nodecount(nsubg) == 4
    @test edgecount(nsubg) == 6
    esubg = edgesubgraph(graph, Set([1, 2, 3, 4]))
    @test nodecount(esubg) == 5
    @test edgecount(esubg) == 4

    subgsubg = nodesubgraph(nsubg, Set([2, 3, 4]))
    @test issetequal(nodeset(subgsubg), [2, 3, 4])
    @test issetequal(edgeset(subgsubg), [2, 3, 9])
    @test degree(subgsubg, 4) == 2
    newg = plaindigraph(subgsubg)
    @test issetequal(nodeset(newg), [1, 2, 3])
    @test issetequal(edgeset(newg), [1, 2, 3])
end

end # graph.inducedsubgraph


text:
module utils

using LinearAlgebra
using Statistics: mean
export check_index,clear_output,_average_helper,convert2array,convert_1d,_check_curve,_calculate_for_curves, _calculate_for_curves_with_matrices, _trapz, _validate_distance_input

global CONVERT_ARRAY_TYPE = true
global ARRAY_TYPE = AbstractArray
global SUPRESS_WARNINGS = false

function check_index(x, none_accepted; class_name = nothing, ith_class = nothing, valid_modes = nothing, average = "macro", weights=nothing)
    if average == "sample-weights"; @assert weights != nothing && length(weights) == length(x) """If the average mode is weighted, weights that are the same size as the labels must be provided!
        @assert average in valid_modes "Unknown averaging mode. This function only supports the following types: " * join(valid_modes, ", ")
    If no precalculated weights can be provided but the class imbalance is to be taken into account try ' average = "weighted" ' or  ' average = "micro" ' """; end
    if !none_accepted; @assert class_name != nothing || ith_class != nothing "No class name or class indexing value provided"; end
    if none_accepted && class_name == nothing == ith_class
        return -1
    elseif class_name != nothing
        @assert class_name in x "There is no such class in the labels of the given confusion matrix"
        index = findfirst(x -> x == class_name, x)
        return index
    else
        @assert ith_class >= 0 && ith_class <= length(x) "ith_class value is not in range"
        return ith_class
    end
end

function clear_output(x, zero_division)
    if true in [isnan(i) || isinf(i) for i in x]
        if zero_division == "warn" || zero_division == "0"
            if zero_division == "warn"; @warn "Zero division, replacing NaN or Inf with 0"; end;
            if length(x) > 1
                return replace(x, NaN => 0)
            else
                return 0
            end
        else
            if length(x) > 1
                return replace(x, NaN => 1)
            else
                return 1
            end
        end
    else return x
    end
end

function _average_helper(numerator, denominator, weights, average, zero_division, normalize)
    if average == "macro"
        x = clear_output(denominator == nothing ? mean(numerator) : mean(numerator ./ denominator), zero_division)
    elseif average == "micro"
        x = clear_output(denominator == nothing ? sum(numerator) : sum(numerator) / sum(denominator) , zero_division)
    elseif average == "weighted" || average == "sample-weights"
        x = clear_output(denominator == nothing ? (sum(numerator .* weights) / sum(weights)) : (sum((numerator ./ denominator) .* weights) / sum(weights)), zero_division)
    else
        x = clear_output(denominator == nothing ? numerator : numerator ./ denominator, zero_division)
    end
    return normalize && length(x) > 1 ? LinearAlgebra.normalize(x) : x
end

function convert2array(x)
     if CONVERT_ARRAY_TYPE && !(typeof(x) <: ARRAY_TYPE)
         !SUPRESS_WARNINGS && @warn "Input has been converted to KnetMetrics.ARRAY_TYPE ($ARRAY_TYPE) . If this behaviour is not wanted,
         set KnetMetrics.CONVERT_ARRAY_TYPE to false. Warning can be supressed via setting Knet.SUPRESS_WARNINGS to
         true."
         return ARRAY_TYPE(x)
    else
        return x
    end
end

convert_1d(x) = reshape(convert2array(x), (length(x)))

function _check_curve(y_true, thresholding_func = nothing, thresholds = [], first_func = false_positive_rate, second_func = true_positive_rate; ith_class = nothing,
        class_name = nothing, average = "macro", normalize = false, weights = nothing)
    thresholding_func == nothing || thresholds == [] ? throw(DomainError("'thresholding_func' cannot be equal to nothing and thresholds must be provided.")) : nothing
    y_true, thresholds = convert_1d(y_true),  convert_1d(thresholds)
    @assert average in ["macro", "weighted", "micro", "sample-weights"] "Unknown averaging mode. This function only supports the following types: macro, weighted, sample-weights"
    if average == "sample-weights"; @assert weights != nothing && length(weights) == length(c.Labels) """If the average mode is weighted, weights that are the same size as the labels must be provided!
    If no precalculated weights can be provided but the class imbalance is to be taken into account try ' average = "weighted" ' or  ' average = "micro" ' """; end
end

function _calculate_for_curves(y_true, y_pred, thresholding_func = nothing, thresholds = [], first_func = false_positive_rate, second_func = true_positive_rate; ith_class = nothing,
        class_name = nothing, average = "macro", normalize = false, weights = nothing)
    vals = [thresholding_func(y_pred, i) for i in thresholds]
    matrices = [confusion_matrix(y_true, predictions) for predictions in vals]
    return _calculate_for_curves_with_matrices(matrices, thresholding_func, thresholds, first_func, second_func; ith_class = ith_class,
            class_name = class_name, average = average, normalize = normalize, weights = weights)

end

function _calculate_for_curves_with_matrices(c, thresholding_func = nothing, thresholds = [], first_func = false_positive_rate, second_func = true_positive_rate; ith_class = nothing,
        class_name = nothing, average = "macro", normalize = false, weights = nothing)
    options = (ith_class = ith_class, class_name = class_name, average = average, normalize = normalize, weights = weights)
    ratios = []
    ratios = [(first_func(c[i];options...), second_func(c[i];options...)) for i in 1:length(c)]
    return RocCurve(thresholds, ratios, options)
end

_trapz(y,x, dx = 1.0) = sum(diff(x) .* (y[1:end-1] .+ y[2:end] ) / 2.0 )

function _validate_distance_input(u,v,w; p = nothing, p_is_used = false, check_weight_length=false)
    if p_is_used; @assert p >= 1 "'p' value must be greater than or equal to one"; end
    check_weight_length && w != nothing && @assert length(w) == length(u) == length(v) "Given 'w' values must be of same size with distance vectors"
end

end


text:
"""
Along same lines as  `@test_throws` but to test if a macro throw an exception when it is
expanded.
"""
macro test_macro_throws(err_expr, expr)
    quote
        err = nothing
        try
            @macroexpand($(esc(expr)))
        catch load_err
            # all errors thrown at macro expansion time are LoadErrors, we need to unwrap
            @assert load_err isa LoadError
            err = load_err.error
        end
        # Reuse `@test_throws` logic
        if err!==nothing
            @test_throws $(esc(err_expr)) ($(Meta.quot(expr)); throw(err))
        else
            @test_throws $(esc(err_expr)) $(Meta.quot(expr))
        end
    end
end

@testset "rule_definition_tools.jl" begin
    @testset "@non_differentiable" begin
        @testset "two input one output function" begin
            nondiff_2_1(x, y) = fill(7.5, 100)[x + y]
            @non_differentiable nondiff_2_1(::Any, ::Any)
            @test frule((Zero(), 1.2, 2.3), nondiff_2_1, 3, 2) == (7.5, DoesNotExist())
            res, pullback = rrule(nondiff_2_1, 3, 2)
            @test res == 7.5
            @test pullback(4.5) == (NO_FIELDS, DoesNotExist(), DoesNotExist())
        end

        @testset "one input, 2-tuple output function" begin
            nondiff_1_2(x) = (5.0, 3.0)
            @non_differentiable nondiff_1_2(::Any)
            @test frule((Zero(), 1.2), nondiff_1_2, 3.1) == ((5.0, 3.0), DoesNotExist())
            res, pullback = rrule(nondiff_1_2, 3.1)
            @test res == (5.0, 3.0)
            @test isequal(
                pullback(Composite{Tuple{Float64, Float64}}(1.2, 3.2)),
                (NO_FIELDS, DoesNotExist()),
            )
        end

        @testset "constrained signature" begin
            nonembed_identity(x) = x
            @non_differentiable nonembed_identity(::Integer)

            @test frule((Zero(), 1.2), nonembed_identity, 2) == (2, DoesNotExist())
            @test frule((Zero(), 1.2), nonembed_identity, 2.0) == nothing

            res, pullback = rrule(nonembed_identity, 2)
            @test res == 2
            @test pullback(1.2) == (NO_FIELDS, DoesNotExist())

            @test rrule(nonembed_identity, 2.0) == nothing
        end

        @testset "Pointy UnionAll constraints" begin
            pointy_identity(x) = x
            @non_differentiable pointy_identity(::Vector{<:AbstractString})

            @test frule((Zero(), 1.2), pointy_identity, ["2"]) == (["2"], DoesNotExist())
            @test frule((Zero(), 1.2), pointy_identity, 2.0) == nothing

            res, pullback = rrule(pointy_identity, ["2"])
            @test res == ["2"]
            @test pullback(1.2) == (NO_FIELDS, DoesNotExist())

            @test rrule(pointy_identity, 2.0) == nothing
        end

        @testset "kwargs" begin
            kw_demo(x; kw=2.0) = x + kw
            @non_differentiable kw_demo(::Any)

            @testset "not setting kw" begin
                @assert kw_demo(1.5) == 3.5

                res, pullback = rrule(kw_demo, 1.5)
                @test res == 3.5
                @test pullback(4.1) == (NO_FIELDS, DoesNotExist())

                @test frule((Zero(), 11.1), kw_demo, 1.5) == (3.5, DoesNotExist())
            end

            @testset "setting kw" begin
                @assert kw_demo(1.5; kw=3.0) == 4.5

                res, pullback = rrule(kw_demo, 1.5; kw=3.0)
                @test res == 4.5
                @test pullback(1.1) == (NO_FIELDS, DoesNotExist())

                @test frule((Zero(), 11.1), kw_demo, 1.5; kw=3.0) == (4.5, DoesNotExist())
            end
        end

        @testset "Not supported (Yet)" begin
            # Varargs are not supported
            @test_macro_throws ErrorException @non_differentiable vararg1(xs...)
            @test_macro_throws ErrorException @non_differentiable vararg1(xs::Vararg)

            # Where clauses are not supported.
            @test_macro_throws(
                ErrorException,
                (@non_differentiable where_identity(::Vector{T}) where T<:AbstractString)
            )
        end
    end
end


text:
get_draw_predictions(b::Gtk.GtkBuilder)=get_gtk_property(b["dl_show_predictions"],:active,Bool)

function draw_predictions(han::Tracker_Handles)
    (preds,confidences) = predict_single_frame(han)
    _draw_predicted_whisker(preds[:,1] ./ 64 .* han.w,preds[:,2] ./ 64 .* han.h,confidences,han.c,han.nn.confidence_thres)
end

function draw_predicted_whisker(han::Tracker_Handles)
    d=han.displayed_frame
    x=han.nn.predicted[:,1,d]; y=han.nn.predicted[:,2,d]; conf=han.nn.predicted[:,3,d]
    _draw_predicted_whisker(x,y,conf,han.c,han.nn.confidence_thres)
end

function _draw_predicted_whisker(x,y,c,canvas,thres)

    circ_rad=5.0

    ctx=Gtk.getgc(canvas)
    num_points = length(x)

    for i=1:num_points
        if c[i] > thres
            Cairo.set_source_rgba(ctx,0,1,0,1-0.025*i)
            Cairo.arc(ctx, x[i],y[i], circ_rad, 0, 2*pi);
            Cairo.stroke(ctx);
        end
    end
    reveal(canvas)
end

function calculate_whisker_fit(pred_1::AbstractArray{T,2},img::AbstractArray{N,2}) where {T,N}

    calculate_whisker_fit(pred_1,size(img))
end

function calculate_whisker_fit(pred_1::AbstractArray{T,2},sz::Tuple,n_points_max=100) where T

    (rot_angle, loss) = quick_quad(pred_1)

    upsampled = StackedHourglass.upsample_pyramid(pred_1,sz)

    (x,y,conf) = get_points(upsampled)

    x = convert(Array{Float64,1},x)
    y = convert(Array{Float64,1},y)

    calculate_whisker_fit(x,y,conf,n_points_max,rot_angle)
end

function get_points(input::AbstractArray{T,2},conf_thres=0.5) where T
    input_s = input./maximum(input)
    points_1 = findall(input_s .> conf_thres)
    x = [points_1[i][1] for i=1:length(points_1)]
    y = [points_1[i][2] for i=1:length(points_1)]
    conf = input_s[points_1]
    (x,y,conf)
end

function quick_quad(pred::AbstractArray{T,2},angles=0.0:pi/12:pi) where T

    inds = findall(pred .> 0.1)

    x = zeros(Float64,length(inds))
    y = zeros(Float64,length(inds))

    for i=1:length(inds)
        x[i] = inds[i][1]
        y[i] = inds[i][2]
    end

    quick_quad(x,y,angles)
end

function quick_quad(x::AbstractArray{T,1},y::AbstractArray{T,1},angles=0.0:pi/12:pi) where T

    x_prime = zeros(Float64,length(x))
    y_prime = zeros(Float64,length(y))

    quad_mat = ones(Float64,length(x),3)

    losses = zeros(Float64,length(angles))

    for i = 1:length(angles)

        rotate_mat(x,y,x_prime,y_prime,angles[i])

        for j=1:length(x_prime)
            quad_mat[j,1] = x_prime[j]
            quad_mat[j,2] = x_prime[j] * x_prime[j]
        end
        out = quad_mat \ y_prime

        losses[i] = mean(((quad_mat * out) .- y_prime).^2)
    end

    (angles[findmin(losses)[2]],minimum(losses))
end

function calculate_whisker_fit(x::AbstractArray{T,1},y::AbstractArray{T,1},conf::AbstractArray{N,1},
    n_points_max=100,rot_angle=0.0) where {T,N}

    rotate_mat!(x,y,rot_angle)
    (new_x, new_y) = center_of_mass(x,y,conf,n_points_max)
    rotate_mat!(new_x,new_y,-1 * rot_angle)
    (new_x,new_y)
end

#=

=#
function center_of_mass(x::AbstractArray{T,1},y::AbstractArray{T,1},conf::AbstractArray{N,1},n_points_max=100) where {T,N}
    x_sort=sortperm(x)
    sorted_x=sort(x)
    sorted_y = y[x_sort]
    sorted_conf = conf[x_sort]

    out_length = round(Int64,sorted_x[end] - sorted_x[1])

    if out_length > n_points_max
        out_length = n_points_max
    end

    myedges = range(sorted_x[1],stop=sorted_x[end],length=out_length)

    out_x = zeros(Float64,length(myedges)-1)
    out_y = zeros(Float64,length(myedges)-1)
    out_weights = zeros(Float64,length(myedges)-1)

    start_ind = 1

    for i=1:(length(myedges)-1)
        a = myedges[i]
        b = myedges[i+1]
       for j=start_ind:length(x)
            this_x = sorted_x[j]
            if (a <= this_x < b)
                out_y[i] = out_y[i] + sorted_y[j] * sorted_conf[j]
                out_x[i] = out_x[i] + this_x * sorted_conf[j]
                out_weights[i] = out_weights[i] + sorted_conf[j]
            elseif (this_x >= b)
                start_ind = j
                break
            end
        end
    end

    out_y = out_y ./ out_weights
    out_x = out_x ./ out_weights

    nan_vals = isnan.(out_x) .| isnan.(out_y)
    deleteat!(out_x,nan_vals)
    deleteat!(out_y,nan_vals)

    (out_x, out_y)
end

function calculate_whisker_predictions(han::Tracker_Handles,hg)
    pred=StackedHourglass.predict_single_frame(hg,han.current_frame./255)
end

#=
Rotate matrices about specified angle
=#
function rotate_mat(x::AbstractArray{T,1},y::AbstractArray{T,1},x_prime::AbstractArray{T,1},y_prime::AbstractArray{T,1},theta) where T
    x_prime[:] = x .* cos(theta) .- y .* sin(theta)
    y_prime[:] = y .* cos(theta) .+ x .* sin(theta)
    nothing
end

function rotate_mat(x::AbstractArray{T,1},y::AbstractArray{T,1},theta::Real) where T
    x_prime = x .* cos(theta) .- y .* sin(theta)
    y_prime = y .* cos(theta) .+ x .* sin(theta)
    (x_prime,y_prime)
end

function rotate_mat!(x::AbstractArray{T,1},y::AbstractArray{T,1},theta::Real) where T
    for i=1:length(x)
        x_old = x[i]
        y_old = y[i]
        (x[i],y[i]) = rotate_mat(x_old,y_old,theta)
    end
    nothing
end

function rotate_mat(points::Array{CartesianIndex{2},1},theta::Real)
    x = zeros(Float64,length(points))
    y = zeros(Float64,length(points))
    for i=1:length(x)
        (x[i],y[i])=rotate_mat(points[i][1],points[i][2],theta)
    end
    (x,y)
end

function rotate_mat(x::Real,y::Real,theta::Real)
    x_out = x * cos(theta) - y * sin(theta)
    y_out = y * cos(theta) + x * sin(theta)
    (x_out,y_out)
end

function draw_prediction2(han::Tracker_Handles,hg,conf)

    colors=((1,0,0),(0,1,0),(0,1,1),(1,0,1))
    pred = calculate_whisker_predictions(han,hg)
    for i = 1:size(pred,3)
        (x,y) = calculate_whisker_fit(pred[:,:,i,1],han.current_frame)
        draw_points_2(han,y,x,colors[i])
    end

    reveal(han.c)
end

function draw_points_2(han::Tracker_Handles,x::Array{T,1},y::Array{T,1},cc) where T
    ctx=Gtk.getgc(han.c)

    set_source_rgb(ctx,cc...)

    set_line_width(ctx, 1.0);
    for i=1:length(x)
        arc(ctx, x[i],y[i], 1.0, 0, 2*pi);
        stroke(ctx);
    end
end

function get_whisker(upsampled::AbstractArray{T,2},rot_angle::Real,n_points=100) where T

    (xx,yy,confs) = WhiskerTracking.get_points(upsampled)

    xx = convert(Array{Float64,1},xx)
    yy = convert(Array{Float64,1},yy)

    (x_out,y_out) = WhiskerTracking.calculate_whisker_fit(xx,yy,confs,n_points,rot_angle)
end