text:
using PyCall
import SymPy: symbols, sympy, Sym
using GAlgebra
using Test

py"""
def vector(ga, components):
    basis = ga.mv()
    return sum([components[i] * e for i, e in enumerate(basis)])
"""
const vector = py"vector"

py"""
def signature(ga):
    basis = ga.mv()
    signs = [e * e for e in basis]
    p, q, r = 0, 0, 0
    for sign in signs:
        p += 1 if sign == 1 else 0
        q += 1 if sign == -1 else 0
        r += 1 if sign == 0 else 0

    return (p, q, r)
"""
const signature = py"signature"

macro ctest(assertion::Expr, sep::QuoteNode, condition::Expr)
    if sep == :(:if)
        esc(quote
            if $condition
                @test $assertion
            else
                @test_broken $assertion
            end
        end)
    elseif sep == :(:unless)
        esc(quote
            if $condition
                @test_broken $assertion
            else
                @test $assertion
            end
        end)
    else
        throw(DomainError(sep, "sep can only be one of :if, :unless"))
    end
end

@testset "GAlgebra.jl" begin
    # Basic
    Hyper   = G(1)       # Hyperbolic numbers.
    ℂ       = G(0,1)     # Complex numbers.
    Dual    = G(0,0,1)   # Dual numbers.
    ℍ       = G(0,2)     # Quaternions.

    # Clifford
    Cl2       = G(2)     # Clifford algebra for 2D vector space.
    Cl3       = G(3)     # Clifford algebra for 3D vector space.
    Spacetime = G(1,3)   # Clifford algebra for timespace vectors.

    # Geometric
    PGA2D = G(2,0,1)     # Projective Euclidean 2D plane. (dual)
    PGA3D = G(3,0,1)     # Projective Euclidean 3D space. (dual)
    CGA2D = G(3,1)       # Conformal 2D space.
    CGA3D = G(4,1)       # Conformal 3D space.

    function test_all(V)
        dimV = range(0, stop=V.n)
        I = V.I()

        α = V.mv("α", "scalar")
        β = V.mv("β", "scalar")
        γ = V.mv("γ", "scalar")
        λ = V.mv("λ", "scalar")

        u = V.mv("u", "vector")
        v = V.mv("v", "vector")
        w = V.mv("w", "vector")

        A = V.mv("A", "mv")
        B = V.mv("B", "mv")
        C = V.mv("C", "mv")
        D = V.mv("D", "mv")

        R = V.mv("R", "spinor")

        # Precalculte AB and BA
        AB = A * B
        BA = B * A

        # The following tests verified implementation correctness per definition

        @test u ⋅ v == u | v == (u < v) == (u > v) == u ⨼ v == u ⨽ v == u ⊙ v
        @test u ∧ v == u ⊠ v
        @test v ⨼ B == (v < B)
        @test v ⨽ B == (v > B)
        if V ∉ [PGA3D, CGA3D] # too slow
            @test A ⊙ B == A << B == (AB + BA) / 2
            @test A ⊠ B == A >> B == (AB - BA) / 2
            @test A ⊛ B == A % B
        end

        @test abs(v) == norm(v) == v.norm()
        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D]
            @test abs(R) == norm(R) == R.norm()
        end

        @test ~A == A[:~] == rev(A) == A.rev()

        if V ∉ [Dual, PGA2D, PGA3D, CGA2D, CGA3D]
            @test A' == dual(A) == A.dual() == adjoint(A) == A * I # Ga.dual_mode_value is default to "I+"
            @test (v)⁻¹ == v[:⁻¹] == v^-1 == inv(v) == v.inv()
            @test v^-2 == (v^2).inv()
        end

        @test (A)ˣ == A[:*] == involute(A) == (A)₊ - (A)₋ == A[:+] - A[:-] == A.even() - A.odd()
        @test (A)ǂ == A[:ǂ] == conj(A) == involute(A).rev()

        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D]
            @test R^-2 == (R^2).inv()
            @test (R)⁻¹ == R[:⁻¹] == R^-1 == inv(R) == R.inv()
            @test ((R)⁻¹)ˣ == ((R)ˣ)⁻¹
            @test ((R)⁻¹)ǂ == ((R)ǂ)⁻¹
        end

        if V ∈ [Cl2, Cl3]
            @test (v)⁻¹ == (~v) / norm(v)^2 == v / v^2
            @test (R)⁻¹ == (~R) / norm(R)^2 == R / R^2
        end

        @test v^0 == 1
        @test v^2 == v*v

        @test ((A)ˣ)ˣ == ~(~A) == A[:~][:~] == ((A)ǂ)ǂ == A
        @test ~((A)ˣ) == (~A)ˣ

        if V ∉ [Dual]
            @test proj(u, v) == v.project_in_blade(u)
            @test refl(u, v) == v.reflect_in_blade(u)
        end

        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
            @test rot(u ∧ v, A) == A.rotate_multivector(u ∧ v)
            @test exp(u ∧ v) == (u ∧ v).exp()
        end

        @test typeof(scalar(A)) == Sym
        @test typeof(A[0]) == Mv
        @test scalar(A) == A.scalar() == A[0].obj
        @test (A)₊ == A[:+] == even(A) == A.even()
        @test (A)₋ == A[:-] == odd(A) == A.odd()

        for r ∈ dimV
            A[r] == A.grade(r) == A.get_grade(r)
        end

        # The following tests verified many identities in Linear Algebra

        @test v + w == w + v
        @test (u + v) + w == u + (v + w)
        @test v + 0 == v
        @test 0 * v == 0
        @test 1 * v == v
        @test α * (β * v) == (α * β) * v
        @test α * (v + w) == α * v + α * w
        @test (α + β) * v == α * v + β * v
        @test v + (-1) * v == 0
        @test -v == -1 * v

        𝑶 = vector(V, fill(0, V.n))
        @test α * 𝑶 == 𝑶
        @test (-α) * v == α * (-v) == -α * v

        # The following tests verified many identities in https://arxiv.org/abs/1205.5935

        # Summary
        #
        # - The formulas after XXX are work-in-progress
        # - The following formulas are skip for now
        #   - A.4.27
        # - The following formulas are broken for now
        #   - A.4.{33, 35}

        @test v * v == (v * v).scalar()
        @test v * B == v ⋅ B + v ∧ B == v ⨼ B + v ∧ B

        @test u ∧ (v + λ * u) == u ∧ v

        @test v == v[1]
        if V.n >= 2
            G2 = V.mv("G2", "grade", 2)
            @test G2 == G2[2]
        end

        for r ∈ dimV

            @test (A + B)[r] == A[r] + B[r]
            @test (λ * A)[r] == (A * λ)[r] == λ * A[r]

            Ar = A[r]

            @test v ⨼ Ar == (v * Ar - (-1)^r * Ar * v) / 2
            @test Ar ⨽ v == (Ar * v - (-1)^r * v * Ar) / 2 == (-1)^(r-1) * (v ⨼ Ar)
            @test v ∧ Ar == (v * Ar + (-1)^r * Ar * v) / 2
            @test Ar ∧ v == (Ar * v + (-1)^r * v * Ar) / 2 == (-1)^r * (v ∧ Ar)

            @test v ⨼ Ar == (v * Ar)[r-1]
            @test v ∧ Ar == (v * Ar)[r+1]
            @test Ar ⨽ v == (Ar * v)[r-1]
            @test Ar ∧ v == (Ar * v)[r+1]
            @test v * Ar == v ⨼ Ar + v ∧ Ar
            @test Ar * v == Ar ⨽ v + Ar ∧ v

            Br = B[r]
            Ar ⨼ Br == Ar ⨽ Br == (Ar * Br).scalar()

            if r > 0
                a = V.mv("a", "vector")
                a₁₋ᵣ = [V.mv("a_$(i)", "vector") for i in 1:r]

                wedge_a₁₋ᵣ = reduce(∧, a₁₋ᵣ)
                prod_a₁₋ᵣ = reduce(*, a₁₋ᵣ)

                @test wedge_a₁₋ᵣ == prod_a₁₋ᵣ[r]                                                  # A.4.12

                A₁₋ᵣ = prod_a₁₋ᵣ

                for s ∈ dimV
                    Bs = B[s]
                    A_Bs_Aǂ = A₁₋ᵣ * Bs * (A₁₋ᵣ)ǂ

                    @test A_Bs_Aǂ == A_Bs_Aǂ[s]                                                                              # A.4.32
                end

                # TODO this is failing for now
                # @test (A₁₋ᵣ * B * (A₁₋ᵣ)ǂ) ∧ (A₁₋ᵣ * C * (A₁₋ᵣ)ǂ) == A₁₋ᵣ.norm()^2 * A₁₋ᵣ * (B ∧ C) * (A₁₋ᵣ)ǂ        # A.4.33

                if r > 1
                    function reduce_except(op, arr, j)
                        enumerate(arr) |> o -> Iterators.filter(e -> e[1] != j, o) |>
                            o -> Iterators.map(e-> e[2], o) |> o -> reduce(op, o)
                    end

                    @test a ⨼ wedge_a₁₋ᵣ == sum([(-1)^(j-1) * (a ⨼ a₁₋ᵣ[j]) * reduce_except(∧, a₁₋ᵣ, j) for j in 1:r])  # A.4.25
                    @test a₁₋ᵣ[1] ∧ reduce_except(∧, a₁₋ᵣ, 1) == wedge_a₁₋ᵣ                                             # A.4.26
                end
            end

            for s ∈ dimV
                @test A[r][s] == (if r == s; A[r] else 0 end)

                Bs = B[s]
                ArBs = Ar * Bs

                @test ArBs == sum([ArBs[abs(r - s) + 2j] for j=0:min(r, s)])                      # A.4.1
                @test Ar ⨼ Bs == (-1)^(r * (s - 1)) * Bs ⨽ Ar                                    # A.4.10
                @test Ar ∧ Bs == (-1)^(r * s) * Bs ∧ Ar                                          # A.4.11

                for j ∈ dimV
                    @test ArBs[r + s - 2j] == (-1)^(r * s - j) * (B[s] * A[r])[r + s - 2j]        # A.4.2
                end

                if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
                    @test v ⨼ ArBs == (v * ArBs - (-1)^(r+s) * ArBs * v)/2 ==
                        (v ⨼ Ar) * Bs + (-1)^r * Ar * (v ⨼ Bs) ==
                        (v ∧ Ar) * Bs - (-1)^r * Ar * (v ∧ Bs)
                    @test v ∧ ArBs == (v * ArBs + (-1)^(r+s) * ArBs * v)/2 ==
                        (v ∧ Ar) * Bs - (-1)^r * Ar * (v ⨼ Bs) ==
                        (v ⨼ Ar) * Bs + (-1)^r * Ar * (v ∧ Bs)
                end

                @test v ⨼ (Ar ∧ Bs) == (v ⨼ Ar) ∧ Bs + (-1)^r * Ar ∧ (v ⨼ Bs)
                @test v ∧ (Ar ⨽ Bs) == (v ∧ Ar) ⨽ Bs - (-1)^r * Ar ⨽ (v ⨼ Bs)
                @test v ∧ (Ar ⨼ Bs) == (v ⨼ Ar) ⨼ Bs + (-1)^r * Ar ⨼ (v ∧ Bs)

                if r > s
                    @test Ar ⨼ Bs == Bs ⨽ Ar == 0
                end

                if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
                    for t ∈ dimV
                        Ct = C[t]

                        Ar ∧ (Bs ∧ Ct) == (Ar * Bs * Ct)[r + s + t]
                    end
                end
            end
        end

        @test A == sum([A[r] for r ∈ dimV])
        @test A[-3] == 0

        @test v ⨼ A == (v * A - (A)ˣ * v)/2                                                      # A.4.13
        @test v ∧ A == (v * A + (A)ˣ * v)/2                                                      # A.4.14
        @test A ⨽ v == - v ⨼ (A)ˣ                                                                # A.4.15
        @test A ∧ v == v ∧ (A)ˣ                                                                  # A.4.16

        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
            @test v ⨼ (AB) == (v ⨼ A) * B + (A)ˣ * (v ⨼ B) == (v ∧ A) * B - (A)ˣ * (v ∧ B)       # A.4.18-19
            @test v ∧ (AB) == (v ∧ A) * B - (A)ˣ * (v ⨼ B) == (v ⨼ A) * B + (A)ˣ * (v ∧ B)       # A.4.20-21
        end

        @test v ⨼ (A ∧ B) == (v ⨼ A) ∧ B + (A)ˣ ∧ (v ⨼ B)                                       # A.4.22
        @test v ∧ (A ⨽ B) == (v ∧ A) ⨽ B - (A)ˣ ⨽ (v ⨼ B)                                       # A.4.23
        @test v ∧ (A ⨼ B) == (v ⨼ A) ⨼ B + (A)ˣ ⨼ (v ∧ B)                                       # A.4.24

        @test v ⨼ A[:+] == - (A[:+] ⨽ v)
        @test v ⨼ A[:-] == A[:-] ⨽ v
        @test v ∧ A[:+] == A[:+] ∧ v
        @test v ∧ A[:-] == - (A[:-] ∧ v)

        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
            @test (AB).scalar() == (BA).scalar() == (~A * ~B).scalar() ==
                ((A)ˣ * (B)ˣ).scalar() == ((A)ǂ * (B)ǂ).scalar()                                      # A.4.3-6
        end

        @test A ⨼ B == sum([sum([(A[r] * B[s])[s - r] for r ∈ dimV]) for s ∈ dimV])             # A.4.7
        @test A ⨽ B == sum([sum([(A[r] * B[s])[r - s] for r ∈ dimV]) for s ∈ dimV])             # A.4.8
        @test A ∧ B == sum([sum([(A[r] * B[s])[r + s] for r ∈ dimV]) for s ∈ dimV])             # A.4.9

        @test (A ∧ B) ∧ C == A ∧ (B ∧ C) == A ∧ B ∧ C                                       # A.4.28
        @test A ⨼ (B ⨽ C) == (A ⨼ B) ⨽ C                                                        # A.4.29
        @test A ⨼ (B ⨼ C) == (A ∧ B) ⨼ C                                                        # A.4.30
        @test A ⨽ (B ∧ C) == (A ⨽ B) ⨽ C                                                        # A.4.31
        @test (A ∧ B) ⨼ C == A ⨼ (B ⨼ C)                # commit 3ca67dea999f955c71a5dd2c18d819afbb971e85

        @test u ∧ A ∧ v == - v ∧ A ∧ u                                                           # A.4.17

        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
            @test AB == A ⊙ B + A ⊠ B
            @test A ⊙ B == B ⊙ A
            @test A ⊠ B == - B ⊠ A

            @test A ⊛ B == B ⊛ A

            @test A ⊛ B == ~A ⊛ ~B == A.rev() ⊛ B.rev()
            @test A ⊛ (B * C) == (~B * A) ⊛ C
            @test A ⊛ (B ⨽ C) == (~B ⨽ A) ⊛ C
            @test A ⊛ (B ⨼ C) == (~B ∧ A) ⊛ C
            @test A ⊛ (B ∧ C) == (~B ⨼ A) ⊛ C
        end

        if V ∉ [Spacetime, ℂ, ℍ, Dual, PGA2D, PGA3D, CGA2D, CGA3D]
            @test A ⊛ B == A' ⊛ B' == A.dual() ⊛ B.dual()
        end

        if V ∉ [Spacetime, PGA2D, PGA3D, CGA2D, CGA3D] # too slow
            @test AB ⋅ C ∧ D == ((AB) ⋅ C) ∧ D
        end

        if V ∉ [Dual, PGA2D, PGA3D, CGA2D, CGA3D]
            @test u.dual() == u * V.I()
            @test proj(u, v) == (v ⋅ u) / u == (v ⨼ u) ⨼ u.inv()
            @test proj(w, v) + proj(w, u) == proj(w, u + v)
        end

        if V == Cl3
            @test u × v == -I * (u ∧ v)
            @test_throws PyCall.PyError A × B

            Vr = u ∧ v
            @test proj(Vr, B) == B ⨼ Vr * (Vr)⁻¹ == (B ⨼ Vr) ⨼ (Vr)⁻¹                               # A.4.34
            # TODO this is failing for now
            @test_broken refl(Vr, B) == B ∧ Vr * (Vr)⁻¹ == (B ∧ Vr) ⨽ (Vr)⁻¹                               # A.4.35

            # The following tests verified interoperability with numeric and symbolic numbers
            (ex, ey, ez) = V.mv()

            uu = vector(V, [1, 2, 3])
            vv = vector(V, [4, 5, 6])
            ww = vector(V, [5, 6, 7])

            @test uu + vv == 5 * ex + 7 * ey + 9 * ez
            @test 7 * uu + 2 * ww == 17 * ex + 26 * ey + 35 * ez
            @test 7 * uu - 2 * ww == -3 * ex + 2 * ey + 7 * ez
            @test 3 * uu + 2 * vv + ww == 16 * ex + 22 * ey + 28 * ez
            @test (sympy.sqrt(2) * u + sympy.Rational(2, 3) * v) ⋅ ey ==
                sympy.sqrt(2) * (u ⋅ ey) + sympy.Rational(2, 3) * (v ⋅ ey)
        end
    end

    # [PGA2D, PGA3D, CGA2D, CGA3D] pass but take extremely long
    # even PGA2D takes 189.979728 seconds (40.41 k allocations: 1.395 MiB)
    # if not skipping the slow bit
    # CGA3D takes 152.269892 seconds (2.75 M allocations: 126.601 MiB)
    # even skipping the slow bit

    # [Hyper, Dual, ℂ] # for quick test any of p, q, r is 1

    setV = if "CI" in keys(ENV)
        if ENV["GITHUB_EVENT_NAME"] == "schedule"
            [Cl2, Cl3, ℂ, ℍ, Hyper, Dual, Spacetime, PGA2D, PGA3D, CGA2D] #, CGA3D]
        else
            [Cl2, Cl3, ℂ, ℍ, Hyper, Dual, Spacetime]
        end
    else
        [Cl2, Cl3, ℂ, ℍ, Hyper, Dual]
    end

    for V ∈ setV
        sigV = signature(V)
        @time @testset "G$sigV" begin
            # @test_broken 1==1 # for triggering error
            test_all(V)
        end

        # The following takes forever

        # High-Dimensional GA
        # DCGA3D = G(6,2)      # Double Conformal 3D Space.
        # TCGA3D = G(9,3)      # Tripple Conformal 3D Space.
        # DCGSTA = G(4,8)      # Double Conformal Geometric Space Time Algebra.
        # QCGA   = G(9,6)      # Quadric Conformal Geometric Algebra.
    end
end


text:
module AlignUtils

using Pkg

haskey(Pkg.installed(),"DISSEQT") || warn("Module DISSEQT not installed. Please refer to AlignUtils installation instructions at https://github.com/rasmushenningsson/AlignUtils.jl")
haskey(Pkg.installed(),"SynapseClient") || warn("Module SynapseClient not installed. Please refer to AlignUtils installation instructions at https://github.com/rasmushenningsson/AlignUtils.jl")
haskey(Pkg.installed(),"SynapseTools") || warn("Module SynapseTools not installed. Please refer to AlignUtils installation instructions at https://github.com/rasmushenningsson/AlignUtils.jl")
haskey(Pkg.installed(),"BamReader") || warn("Module BamReader not installed. Please refer to AlignUtils installation instructions at https://github.com/rasmushenningsson/AlignUtils.jl")

using Distributed
using DataFrames
using Gadfly
using Colors
using JLD
using Levenshtein
haskey(Pkg.installed(),"SynapseClient") && using SynapseClient
haskey(Pkg.installed(),"SynapseTools") && using SynapseTools
using BamReader
using DISSEQT

export
	Sample,
	find_samples,
	find_aligned,
	assign_reference!,
	align_sample!,
	align_samples!,
	askforconfirmation,
	makecleanfolder,
	getreferenceinfo,
	reference_sanity_check,
	computecodonfrequencies,
	computenucleotidefrequencies,
	uploadaligned,
	uploadswarms,
	coverageplots

include("log.jl")
include("sample.jl")
include("synapseutils.jl")
include("align.jl")
include("misc.jl")
include("consensus.jl")
include("codonfrequency.jl")
include("nucleotidefrequency.jl")
include("readcoverage.jl")

end

ext: jl


text:
@testset "990.satisfiability-of-equality-equations.jl" begin
    @test equations_possible(["a==b", "b!=a"]) == false
    @test equations_possible(["b==a", "a==b"]) == true
    @test equations_possible(["a==b", "b==c", "a==c"]) == true
    @test equations_possible(["a==b", "b!=c", "c==a"]) == false
    @test equations_possible(["c==c", "b==d", "x!=z"]) == true
end


text:
@testset "testing constructor" begin

    @test Δ_max(test_options) ≡ Δₘ

    @test δ_tol(test_options) ≡ δ

    @test samples(test_options) ≡ 2w

    @test max_iterations(test_options) ≡ N

    @test weave_level(test_options) ≡ WEAVE_LEVEL

    @test log_level(test_options) ≡ LOG_LEVEL

    restricted_type_test(test_options)
end

@testset "testing behavior" begin

    options = DriverOptions()

    # Mutating option fields only if positive
    @test !(Δ_max!(options, 0.0) ≡ 0.0)

    @test !(δ_tol!(options, 0.0) ≡ 0.0)

    @test !(ϵ_tol!(options, 0.0) ≡ 0.0)

    @test !(samples!(options, 0) ≡ 0)

    @test !(samples!(options, 1) ≡ 1)

    @test Δ_max!(options, 1.0) ≡ 1.0

    @test δ_tol!(options, 1.0) ≡ 1.0

    @test ϵ_tol!(options, 1.0) ≡ 1.0

    @test samples!(options, 2) ≡ 2

    @test weave_level!(options, NONE) ≡ NONE

    @test log_level!(options, ERROR) ≡ ERROR
end


text:
using BinaryProvider

const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
const libname = Sys.iswindows() ? "gmsh-4.4" : "libgmsh"

products = Product[
    LibraryProduct(prefix, libname, :libgmsh),
]

bin_prefix = "http://gmsh.info/bin"
version = "4.4.1"

download_info = Dict(
    # since v4.2.3, Linux SDK will cause segment fault if being `dlopen`
    Linux(:x86_64, :glibc) => ("$bin_prefix/Linux/gmsh-4.2.2-Linux64-sdk.tgz", "ea6a6d36da41b9e777111e055c416ffe994d57c7e3debf174b98e4c09b3b33d7"),
    Windows(:x86_64) => ("$bin_prefix/Windows/gmsh-$version-Windows64-sdk.zip", "094207b56e23e462f2e11ffc2d7006f88c641b62fa9d01522f731dcf00e321a9"),
    MacOS(:x86_64) => ("$bin_prefix/MacOSX/gmsh-$version-MacOSX-sdk.tgz", "40c13c22f0bff840fc827e5f4530668b2818c1472593370ebf302555df498f9e"),
)

if haskey(ENV, "GMSH_LIB_PATH")
    products = Product[LibraryProduct(ENV["GMSH_LIB_PATH"], libname, :libgmsh),]
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
else
    if any(!satisfied(p; verbose=verbose) for p in products)
        try
            # download and install binaries
            url, tarball_hash = choose_download(download_info)
            try
                install(url, tarball_hash; prefix=prefix, force=true, verbose=true)
            catch e
                # cannot list content of .zip, manually unzip
                tarball_path = joinpath(prefix, "downloads", basename(url))
                run(`unzip $(tarball_path) -d $(prefix.path)`)
            end

            # strip the top directory
            content_path = joinpath(prefix, splitext(basename(url))[1])
            foreach(
                (x) -> mv(joinpath(content_path, x), joinpath(prefix, x); force=true),
                readdir(content_path)
                )
            rm(content_path; force=true, recursive=true)

            # BinaryProvider will search $prefix/bin instead of $prefix/lib
            if Sys.iswindows()
                dir_path = libdir(products[1].prefix, platform_key_abi())
                lib_path = joinpath(dirname(dir_path), "lib")
                run(`cp -L $(lib_path)/* $(dir_path)`)
            end
        catch e
            if typeof(e) <: ArgumentError
                error("Your platform $(Sys.MACHINE) is not supported by this package!")
            else
                rethrow(e)
            end
        end
    end
    write_deps_file(joinpath(@__DIR__, "deps.jl"), products)
end


text:
using Formatting
using MKL
using Base.LinAlg.BLAS
using Gadfly
using DataFrames
using Colors

# ------------------------------------------------------------------------------
#                              INDIVIDUAL BENCH
# ------------------------------------------------------------------------------

# custom colors
Tblue = RGB(104/255,135/255,255/255)
Tred = RGB(255/255,110/255,103/255)
Tlightgray = RGB(0.7,0.7,0.7)
Tdarkgray = RGB(0.5,0.5,0.5)

"""
`bstat`retrieves basic statistics for a collection of numbers.
'α' is the risk error (defaults to 0.05).
`v_conf` is the *confidence interval* regarding `α` : v = μ +/- δμ with probability 1-α
"""
function bstat{T<:Number}(v::Vector{T}, α=0.05)
    n = length(v)
    μ = mean(v) # mean
    v_min , v_max = extrema(v) # min and max
    σ =  stdm(v, μ) # standard deviation
    δμ = sqrt(2)*erfinv(1-α)*(σ/sqrt(n)) # confidence interval : v = μ +/- δμ
    return μ, v_min, v_max, σ, δμ
end

"""
`ìbench` performs a benchmark on a function `f` with arguments `arg`.
`f` is evaluated sevral times so the benchmark last approximatively `timing`.
`JIT` is triggered and results are stored in the return arrays at item `1`.
`gcdisable` allows to disable the GC during evaluations (defaults to `true`).
"""
function ibench(f, args, timing=0.0001, gcdisable=true ; α=0.05, print_header = true, print_eval=false, print_stats=false, print_detail=false)

    # TRIGGER JIT : first compilation will be ignored
    res = @timed f(args...)

    # eval nrun to match timing
    t = @elapsed f(args...)
    nrun = max(2,convert(Int64,div(timing,t))) # at least one run

    # allocate results arrays
    res_cpu = Array(Float64,nrun+1)
    res_gc = Array(Float64,nrun+1)
    res_alloc = Array(Int64,nrun+1)

    # write JIT results
    res_cpu[1] = res[2]
    res_gc[1] = res[4]
    res_alloc[1] = res[3]

    # write PROFILING results : call f n times
    gcdisable ? gc() : ""
    for i=1:nrun
        # gcdisable ? gc_enable(false) : ""
        res = @timed f(args...)
        # gcdisable ? gc_enable(true) : ""
        res_cpu[i+1] = res[2]
        res_gc[i+1] = res[4]
        res_alloc[i+1] = res[3]
    end

    # if print_header==true || print_stats==true || print_detail==true
    #     idisplay(f, args, timing, res_cpu, res_gc, res_alloc ; α=α, print_eval=print_eval, print_stats=print_stats, print_detail=print_detail)
    # end
    return res_cpu, res_gc, res_alloc
end

"""
`ìdisplay` prints the results of an individual bench `ìbench` performed on a
function `f` with the given input arguments `arg`.
`detail` printing defaults to `false`.
`res` prints a sample evaluation (for checking) is `true` ; defaults to `false`
"""
function idisplay(f, args, timing, res_cpu::Vector{Float64}, res_gc::Vector{Float64}, res_alloc::Vector{Int64} ; α=0.05, print_eval=false, print_stats=true, print_detail=false)

    nrun = length(res_cpu)-1

    # FORMATER
    fe_sep1 = FormatExpr("{:=>44}")
    fe_sep2 = FormatExpr("{:->44}")
    fe_head = FormatExpr("{:<10s}{:>12s}{:>8s}{:>14s}")
    fe_line = FormatExpr("{:<10s}{:>12.2e}{:>8.2f}{:>14d}")

    # GET STATS
    stat_cpu = bstat(res_cpu[2:end],α)
    stat_gc = bstat(res_gc[2:end],α)
    stat_alloc = bstat(res_alloc[2:end],α)

    # SUMMARY
    println("")
    printfmtln(fe_sep1,"")
    println("Running $f for about $timing s")
    println("Evaluations : $nrun")
    println("CPU time : ",sum(res_cpu),"s")
    # PRINT EVAL
    if print_eval
        println(string(f)," = ")
        println(f(args...))
    end

    # HEADER
    printfmtln(fe_sep1,"")
    printfmtln(fe_head,"NAME","CPU (s)", "GC (%)", "ALLOC (byt)")
    printfmtln(fe_sep2,"")

    # RESULTS : MEAN
    conf = convert(Int32,100*(1-α))
    printfmtln(fe_line, "MEAN", stat_cpu[1], stat_gc[1], stat_alloc[1])
    if print_stats
        printfmtln(fe_line, "MIN", stat_cpu[2], stat_gc[2], stat_alloc[2])
        printfmtln(fe_line, "MAX", stat_cpu[3], stat_gc[3], stat_alloc[3])
        printfmtln(fe_line, "STD", stat_cpu[4], stat_gc[4], stat_alloc[4])
        printfmtln(fe_line, "CONF($conf%)", stat_cpu[5], stat_gc[5], stat_alloc[5])
        printfmtln(fe_sep2,"")
    end

    # RESULTS : JIT
    printfmtln(fe_line, "JIT", res_cpu[1], res_gc[1], res_alloc[1])
    printfmtln(fe_sep2,"")

    # RESULTS : details
    if print_detail == true
        for i=1:nrun
            printfmtln(fe_line, i, res_cpu[i+1], res_gc[i+1], res_alloc[i+1])
        end
    end
    return
end

function idisplay2(f, args, n, timing, neval, res_cpu::Vector{Float64}, res_gc::Vector{Float64}, res_alloc::Vector{Int64}; print_eval=false, print_scale=true)

    # FORMATER
    fe_sep1 = FormatExpr("{:=>44}")
    fe_sep2 = FormatExpr("{:->44}")
    fe_head = FormatExpr("{:<10s}{:>12s}{:>8s}{:>14s}")

    # SUMMARY
    println("")
    printfmtln(fe_sep1,"")
    println("Running '$f' for about $timing s")
    println("Evaluations : $neval")
    # PRINT EVAL
    if print_eval
        println(string(f)," = ")
        println(f(args...))
    end

    # HEADER

    if print_scale
        fe_line = FormatExpr("{:<10s}{:>12.1f}{:>8.1f}{:>14d}")
        printfmtln(fe_sep1,"")
        printfmtln(fe_head,"n=$n","CPU (ns/el)", "GC (%)", "ALLOC (by/el)")
        printfmtln(fe_sep2,"")
        printfmtln(fe_line, "MEAN", res_cpu[2]*(1e9/n), res_gc[2]*(1e9/n), res_alloc[2]÷n)
        printfmtln(fe_line, "JIT", res_cpu[1]*(1e9/n), res_gc[1]*(1e9/n), res_alloc[1]÷n)
        printfmtln(fe_sep2,"")
    else
        fe_line = FormatExpr("{:<10s}{:>12.1e}{:>8.1f}{:>14d}")
        printfmtln(fe_sep1,"")
        printfmtln(fe_head,"n=$n","CPU (s)", "GC (%)", "ALLOC (by)")
        printfmtln(fe_sep2,"")
        printfmtln(fe_line, "MEAN", res_cpu[2], res_gc[2], res_alloc[2])
        printfmtln(fe_line, "JIT", res_cpu[1], res_gc[1], res_alloc[1])
        printfmtln(fe_sep2,"")
    end

    return
end

function ibench2(f, args, n, nrep=500, gcdisable=true)

    # TRIGGER JIT : first compilation will be ignored
    res = @timed f(args...)

    # eval nrun to match timing
    t = @elapsed f(args...)

    # allocate results arrays
    res_cpu = Array(Float64,2)
    res_gc = Array(Float64,2)
    res_alloc = Array(Int64,2)

    # write JIT results
    res_cpu[1] = res[2]
    res_gc[1] = res[4]
    res_alloc[1] = res[3]

    tcpu = 0.0
    tgc = 0.0
    alloc = 0

    gc()

    if n <= 10_000
        nchunk = 10_000 ÷ n + 1
        @inbounds for j in 1:nrep
            res = @timed for k in 1:nchunk f(args...) end
            tcpu += res[2]
            tgc += res[4]
            alloc += res[3]
        end
    else
        nchunk = 1.0
        @inbounds for j in 1:nrep
            res = @timed f(args...)
            tcpu += res[2]
            tgc += res[4]
            alloc += res[3]
        end
    end

    res_cpu[2] = tcpu / (nrep * nchunk)
    res_gc[2] = tgc / (nrep * nchunk)
    res_alloc[2] = alloc ÷ (nrep * nchunk)

    return res_cpu, res_gc, res_alloc
end

# ------------------------------------------------------------------------------
#                              FUNCTIONS
# ------------------------------------------------------------------------------

# FOR NON EXISTING BLAS FUNCTIONS
nanf() = ""

# DEFINE VECTORIZED JULIA FONCTIONS
vj_dot(a,b) = sum(a.*b)
vj_axpy(λ,a,b) = λ*a+b
vj_invsqrt{T<:Number}(a::Vector{T}) = 1 ./ sqrt(a)

# sqrt(X)
function lj_sqrt{T<:Number}(dest,x::T)
    for i in eachindex(x) dest[i] = sqrt(x[i]) end
    dest
end
function ilj_sqrt{T<:Number}(dest,x::T)
    @inbounds for i in eachindex(x) dest[i] = sqrt(x[i]) end
    dest
end


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

mutable struct RouteNextHopType2 <: SwaggerModel

    function RouteNextHopType2(;)
        o = new()
        o
    end
end # type RouteNextHopType2

const _property_map_RouteNextHopType2 = Dict{Symbol,Symbol}()
const _property_types_RouteNextHopType2 = Dict{Symbol,String}()
Base.propertynames(::Type{ RouteNextHopType2 }) = collect(keys(_property_map_RouteNextHopType2))
Swagger.property_type(::Type{ RouteNextHopType2 }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_RouteNextHopType2[name]))}
Swagger.field_name(::Type{ RouteNextHopType2 }, property_name::Symbol) =  _property_map_RouteNextHopType2[property_name]

function check_required(o::RouteNextHopType2)
    true
end

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


text:
function bob(stimulus)
    stripulus = strip(stimulus)
    question = "Sure."
    qyells = "Calm down, I know what I'm doing!"
    yells = "Whoa, chill out!"
    silences = "Fine. Be that way!"
    if last(stripulus, 1) == "?"
        if stripulus == uppercase(stripulus) && occursin(r"[a-zA-Z]", stripulus)
        #This can't be the easiest way to evaluate all caps. Please advise..
            return qyells
            end
        return question
    elseif stripulus == ""
        return silences
    elseif stripulus == uppercase(stripulus) && occursin(r"[a-zA-Z]", stripulus)
    #This can't be the easiest way to evaluate all caps. Please advise.
        return yells
    else
        return "Whatever."
    end
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
function colors!(context::NamedTuple,labels::DataFrame,groups::DataFrame,colors::NamedTuple)
    request = context.request
	try
		response = HTTP.Response(200)
        label = JSON.parse(String(request.body))
        t = @elapsed colors!(label,labels,groups,colors)

		@info """$(request.method) $(replace(request.target, r"\?.+" => "")) | $t seconds"""
		return response

	catch exception
		printstyled("ERROR $(replace(request.target, r"\?.+" => "")) | $exception\n",color=:red)
		return HTTP.Response(500)
	end
end

function colors!(label::Dict,labels::DataFrame,groups::DataFrame,colors::NamedTuple)
	colors.labels.rows[1:3, label["name"] ∈ Base.names(labels) ? labels[!,label["name"]] : groups[!,label["name"]]] .= channelview([parse(RGB,label["color"])])
	return colors
end

function toIndex(x::AbstractVector{<:Union{Number,Missing}};nlevels::Int=10,p::Real=0.1)

	min,max = quantile(skipmissing(x),(p,1-p))
	scaled = ( x .- min ) / ( max - min )

	@. scaled[(scaled<0)&(~ismissing(scaled))] = 0
	@. scaled[(1<scaled)&(~ismissing(scaled))] = 1

	scaled = @. trunc(Union{Int,Missing},(nlevels-1)*scaled) + 1
	return replace!(scaled,missing=>nlevels+1)
end

max_stars_repo_name: gszep/flow-explorer

max_issues_repo_name: gszep/flow-explorer

max_forks_repo_name: gszep/flow-explorer


text:
# # Phillips model of Baroclinic Instability
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/multilayerqg_2layer.ipynb).
#
# A simulation of the growth of barolinic instability in the Phillips 2-layer model
# when we impose a vertical mean flow shear as a difference ``\Delta U`` in the
# imposed, domain-averaged, zonal flow at each layer.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Plots, Printf"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using NPZ

using GeophysicalFlows, Plots, Printf

using Random: seed!

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical parameters and time-stepping parameters

n = 128                  # 2D resolution = n²
stepper = "FilteredRK4"  # timestepper
     dt = 2.5e-3           # timestep
 nsteps = 200000          # total number of time-steps
 nsubs  = 100             # number of time-steps for plotting (nsteps must be multiple of nsubs)
nothing # hide

# ## Physical parameters
L = 2π                   # domain size
μ = 5e-2                 # bottom drag
β = 5                    # the y-gradient of planetary PV
ν = 0                 # hyperviscosity

nlayers = 2              # number of layers
f₀, g = 1, 1             # Coriolis parameter and gravitational constant
 H = [0.2, 0.8]          # the rest depths of each layer
 ρ = [4.0, 5.0]          # the density of each layer

 U = zeros(nlayers) # the imposed mean zonal flow in each layer
 U[1] = 1.0
 U[2] = 0.0
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. In this example we don't
# do any dealiasing to our solution by providing `aliased_fraction=0`.
prob = MultiLayerQG.Problem(nlayers, dev;
                            nx=n, Lx=L, f₀=f₀, g=g, H=H, ρ=ρ, U=U, μ=μ, β=β, ν=ν,
			    dt=dt, stepper=stepper, aliased_fraction=0)
nothing # hide

# and define some shortcuts.
sol, clock, params, vars, grid = prob.sol, prob.clock, prob.params, prob.vars, prob.grid
x, y = grid.x, grid.y
nothing # hide

# ## Setting initial conditions

# Our initial condition is some small-amplitude random noise. We smooth our initial
# condidtion using the `timestepper`'s high-wavenumber `filter`.
#
# `ArrayType()` function returns the array type appropriate for the device, i.e., `Array` for
# `dev = CPU()` and `CuArray` for `dev = GPU()`.

seed!(1234) # reset of the random number generator for reproducibility
q₀  = 1e-2 * ArrayType(dev)(randn((grid.nx, grid.ny, nlayers)))
q₀h = prob.timestepper.filter .* rfft(q₀, (1, 2)) # apply rfft  only in dims=1, 2
q₀  = irfft(q₀h, grid.nx, (1, 2))                 # apply irfft only in dims=1, 2

MultiLayerQG.set_q!(prob, q₀)
nothing # hide

# ## Diagnostics

# Create Diagnostics -- `energies` function is imported at the top.
E = Diagnostic(MultiLayerQG.energies, prob; nsteps=nsteps)
diags = [E] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_2layer"
plotname = "snapshots"
filename = joinpath(filepath, "2layer.jld2")
nothing # hide

# Do some basic file management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# And then create Output
get_sol(prob) = sol # extracts the Fourier-transformed solution
function get_u(prob)
  sol, params, vars, grid = prob.sol, prob.params, prob.vars, prob.grid

  @. vars.qh = sol
  streamfunctionfrompv!(vars.ψh, vars.qh, params, grid)
  @. vars.uh = -im * grid.l * vars.ψh
  invtransform!(vars.u, vars.uh, params)

  return vars.u
end

out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
nothing # hide

startwalltime = time()

u_data = zeros(round(Int, nsteps / nsubs)+1, n, n, nlayers)
v_data = zeros(round(Int, nsteps / nsubs)+1, n, n, nlayers)
q_data = zeros(round(Int, nsteps / nsubs)+1, n, n, nlayers)
ψ_data = zeros(round(Int, nsteps / nsubs)+1, n, n, nlayers)

for j = 0:round(Int, nsteps / nsubs)

  u_data[j+1, :, :, :] = prob.vars.u
  v_data[j+1, :, :, :] = prob.vars.v
  q_data[j+1, :, :, :] = prob.vars.q
  ψ_data[j+1, :, :, :] = prob.vars.ψ

  if j % (1000 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %.1f, cfl: %.2f, KE₁: %.3e, KE₂: %.3e, PE: %.3e, walltime: %.2f min", clock.step, clock.t, cfl, E.data[E.i][1][1], E.data[E.i][1][2], E.data[E.i][2][1], (time()-startwalltime)/60)

    println(log)
  end

  stepforward!(prob, diags, nsubs)
  MultiLayerQG.updatevars!(prob)
end

npzwrite("u_data.npy", u_data)
npzwrite("v_data.npy", v_data)
npzwrite("q_data.npy", q_data)
npzwrite("psi_data.npy", ψ_data)

# ## Save
# Finally, we can save, e.g., the last snapshot via
# ```julia
# savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), clock.step)
# savefig(savename)
# ```


text:
export @var, @new

# Unfortunately there's no good way to deprecate these macros in a way that
# actually yields warnings since they're now parsed during the `crawl` process
# and are never a actually invoked during normal usage.
macro var(x)
    Base.depwarn("Importing @var from JSExpr.jl is deprecated.", :jsexpr_var_deprecation)
    esc(x)
end
macro new(x)
    Base.depwarn("Importing @new from JSExpr.jl is deprecated.", :jsexpr_new_deprecation)
    esc(x)
end


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

module TestDiagonal

using Test, LinearAlgebra, Random
using LinearAlgebra: BlasFloat, BlasComplex

n=12 #Size of matrix problem to test
Random.seed!(1)

@testset for relty in (Float32, Float64, BigFloat), elty in (relty, Complex{relty})
    dd=convert(Vector{elty}, randn(n))
    vv=convert(Vector{elty}, randn(n))
    UU=convert(Matrix{elty}, randn(n,n))
    if elty <: Complex
        dd+=im*convert(Vector{elty}, randn(n))
        vv+=im*convert(Vector{elty}, randn(n))
        UU+=im*convert(Matrix{elty}, randn(n,n))
    end
    D = Diagonal(dd)
    DM = Matrix(Diagonal(dd))

    @testset "constructor" begin
        for x in (dd, GenericArray(dd))
            @test Diagonal(x)::Diagonal{elty,typeof(x)} == DM
            @test Diagonal(x).diag === x
            @test Diagonal{elty}(x)::Diagonal{elty,typeof(x)} == DM
            @test Diagonal{elty}(x).diag === x
            @test Diagonal{elty}(D) === D
        end
        @test eltype(Diagonal{elty}([1,2,3,4])) == elty
        @test isa(Diagonal{elty,Vector{elty}}(GenericArray([1,2,3,4])), Diagonal{elty,Vector{elty}})
        DI = Diagonal([1,2,3,4])
        @test Diagonal(DI) === DI
        @test isa(Diagonal{elty}(DI), Diagonal{elty})
        # issue #26178
        @test_throws MethodError convert(Diagonal, [1, 2, 3, 4])
    end

    @testset "Basic properties" begin
        @test_throws ArgumentError size(D,0)
        @test typeof(convert(Diagonal{ComplexF32},D)) <: Diagonal{ComplexF32}
        @test typeof(convert(AbstractMatrix{ComplexF32},D)) <: Diagonal{ComplexF32}

        @test Array(real(D)) == real(DM)
        @test Array(abs.(D)) == abs.(DM)
        @test Array(imag(D)) == imag(DM)

        @test parent(D) == dd
        @test D[1,1] == dd[1]
        @test D[1,2] == 0

        @test issymmetric(D)
        @test isdiag(D)
        @test isdiag(Diagonal([[1 0; 0 1], [1 0; 0 1]]))
        @test !isdiag(Diagonal([[1 0; 0 1], [1 0; 1 1]]))
        @test istriu(D)
        @test istriu(D, -1)
        @test !istriu(D, 1)
        @test istriu(Diagonal(zero(diag(D))), 1)
        @test istril(D)
        @test !istril(D, -1)
        @test istril(D, 1)
        @test istril(Diagonal(zero(diag(D))), -1)
        if elty <: Real
            @test ishermitian(D)
        end
    end

    @testset "diag" begin
        @test_throws ArgumentError diag(D,  n+1)
        @test_throws ArgumentError diag(D, -n-1)
        @test (@inferred diag(D))::typeof(dd) == dd
        @test (@inferred diag(D, 0))::typeof(dd) == dd
        @test (@inferred diag(D, 1))::typeof(dd) == zeros(elty, n-1)
        DG = Diagonal(GenericArray(dd))
        @test (@inferred diag(DG))::typeof(GenericArray(dd)) == GenericArray(dd)
        @test (@inferred diag(DG, 1))::typeof(GenericArray(dd)) == GenericArray(zeros(elty, n-1))
    end

    @testset "Simple unary functions" begin
        for op in (-,)
            @test op(D)==op(DM)
        end

        for func in (det, tr)
            @test func(D) ≈ func(DM) atol=n^2*eps(relty)*(1+(elty<:Complex))
        end
        if relty <: BlasFloat
            for func in (exp, cis, sinh, cosh, tanh, sech, csch, coth)
                @test func(D) ≈ func(DM) atol=n^3*eps(relty)
            end
            @test log(Diagonal(abs.(D.diag))) ≈ log(abs.(DM)) atol=n^3*eps(relty)
        end
        if elty <: BlasComplex
            for func in (logdet, sqrt, sin, cos, tan, sec, csc, cot,
                         asin, acos, atan, asec, acsc, acot,
                         asinh, acosh, atanh, asech, acsch, acoth)
                @test func(D) ≈ func(DM) atol=n^2*eps(relty)*2
            end
        end
    end

    @testset "Two-dimensional Euler formula for Diagonal" begin
        @test cis(Diagonal([π, π])) ≈ -I
    end

    @testset "Linear solve" begin
        for (v, U) in ((vv, UU), (view(vv, 1:n), view(UU, 1:n, 1:2)))
            @test D*v ≈ DM*v atol=n*eps(relty)*(1+(elty<:Complex))
            @test D*U ≈ DM*U atol=n^2*eps(relty)*(1+(elty<:Complex))

            @test transpose(U)*D ≈ transpose(U)*Array(D)
            @test U'*D ≈ U'*Array(D)

            if relty != BigFloat
                atol_two = 2n^2 * eps(relty) * (1 + (elty <: Complex))
                atol_three = 2n^3 * eps(relty) * (1 + (elty <: Complex))
                @test D\v ≈ DM\v atol=atol_two
                @test D\U ≈ DM\U atol=atol_three
                @test ldiv!(D, copy(v)) ≈ DM\v atol=atol_two
                @test ldiv!(transpose(D), copy(v)) ≈ DM\v atol=atol_two
                @test ldiv!(adjoint(conj(D)), copy(v)) ≈ DM\v atol=atol_two
                @test ldiv!(D, copy(U)) ≈ DM\U atol=atol_three
                @test ldiv!(transpose(D), copy(U)) ≈ DM\U atol=atol_three
                @test ldiv!(adjoint(conj(D)), copy(U)) ≈ DM\U atol=atol_three
                # this method tests AbstractMatrix/AbstractVec for second arg
                Usym_bad = Symmetric(ones(elty, n+1, n+1))
                @test_throws DimensionMismatch ldiv!(D, copy(Usym_bad))

                @test ldiv!(zero(v), D, copy(v)) ≈ DM\v atol=atol_two
                @test ldiv!(zero(v), transpose(D), copy(v)) ≈ DM\v atol=atol_two
                @test ldiv!(zero(v), adjoint(conj(D)), copy(v)) ≈ DM\v atol=atol_two
                @test ldiv!(zero(U), D, copy(U)) ≈ DM\U atol=atol_three
                @test ldiv!(zero(U), transpose(D), copy(U)) ≈ DM\U atol=atol_three
                @test ldiv!(zero(U), adjoint(conj(D)), copy(U)) ≈ DM\U atol=atol_three

                Uc = copy(U')
                target = rmul!(Uc, Diagonal(inv.(D.diag)))
                @test rdiv!(Uc, D) ≈ target atol=atol_three
                @test_throws DimensionMismatch rdiv!(Matrix{elty}(I, n-1, n-1), D)
                @test_throws SingularException rdiv!(Uc, Diagonal(fill!(similar(D.diag), 0)))
                @test rdiv!(Uc, transpose(D)) ≈ target atol=atol_three
                @test rdiv!(Uc, adjoint(conj(D))) ≈ target atol=atol_three
                @test ldiv!(D, Matrix{eltype(D)}(I, size(D))) ≈ D \ Matrix{eltype(D)}(I, size(D)) atol=atol_three
                @test_throws DimensionMismatch ldiv!(D, fill(elty(1), n + 1))
                @test_throws SingularException ldiv!(Diagonal(zeros(relty, n)), copy(v))
                b = rand(elty, n, n)
                @test ldiv!(D, copy(b)) ≈ Array(D)\Array(b)
                @test_throws SingularException ldiv!(Diagonal(zeros(elty, n)), copy(b))
                b = view(rand(elty, n), Vector(1:n))
                b2 = copy(b)
                c = ldiv!(D, b)
                d = Array(D)\b2
                @test c ≈ d
                @test_throws SingularException ldiv!(Diagonal(zeros(elty, n)), b)
                b = rand(elty, n+1, n+1)
                @test_throws DimensionMismatch ldiv!(D, copy(b))
                b = view(rand(elty, n+1), Vector(1:n+1))
                @test_throws DimensionMismatch ldiv!(D, b)
            end
        end
    end
    d = convert(Vector{elty}, randn(n))
    D2 = Diagonal(d)
    DM2= Matrix(Diagonal(d))
    @testset "Binary operations" begin
        for op in (+, -, *)
            @test Array(op(D, D2)) ≈ op(DM, DM2)
        end
        @testset "with plain numbers" begin
            a = rand()
            @test Array(a*D) ≈ a*DM
            @test Array(D*a) ≈ DM*a
            @test Array(D/a) ≈ DM/a
            if elty <: Real
                @test Array(abs.(D)^a) ≈ abs.(DM)^a
            else
                @test Array(D^a) ≈ DM^a
            end
            @test Diagonal(1:100)^2 == Diagonal((1:100).^2)
            p = 3
            @test Diagonal(1:100)^p == Diagonal((1:100).^p)
            @test Diagonal(1:100)^(-1) == Diagonal(inv.(1:100))
            @test Diagonal(1:100)^2.0 == Diagonal((1:100).^2.0)
            @test Diagonal(1:100)^(2.0+0im) == Diagonal((1:100).^(2.0+0im))
        end

        if relty <: BlasFloat
            for b in (rand(elty,n,n), rand(elty,n))
                @test lmul!(copy(D), copy(b)) ≈ Array(D)*Array(b)
                @test lmul!(transpose(copy(D)), copy(b)) ≈ transpose(Array(D))*Array(b)
                @test lmul!(adjoint(copy(D)), copy(b)) ≈ Array(D)'*Array(b)
            end
        end

        #a few missing mults
        bd = Bidiagonal(D2)
        @test D*transpose(D2) ≈ Array(D)*transpose(Array(D2))
        @test D2*transpose(D) ≈ Array(D2)*transpose(Array(D))
        @test D2*D' ≈ Array(D2)*Array(D)'

        #division of two Diagonals
        @test D/D2 ≈ Diagonal(D.diag./D2.diag)
        @test D\D2 ≈ Diagonal(D2.diag./D.diag)

        # QR \ Diagonal
        A = rand(elty, n, n)
        qrA = qr(A)
        @test qrA \ D ≈ A \ D

        # HermOrSym
        A     = rand(elty, n, n)
        Asym  = Symmetric(A + transpose(A), :U)
        Aherm = Hermitian(A + adjoint(A), :U)
        for op in (+, -)
            @test op(Asym, D) isa Symmetric
            @test Array(op(Asym, D)) ≈ Array(Symmetric(op(Array(Asym), Array(D))))
            @test op(D, Asym) isa Symmetric
            @test Array(op(D, Asym)) ≈ Array(Symmetric(op(Array(D), Array(Asym))))
            if !(elty <: Real)
                Dr = real(D)
                @test op(Aherm, Dr) isa Hermitian
                @test Array(op(Aherm, Dr)) ≈ Array(Hermitian(op(Array(Aherm), Array(Dr))))
                @test op(Dr, Aherm) isa Hermitian
                @test Array(op(Dr, Aherm)) ≈ Array(Hermitian(op(Array(Dr), Array(Aherm))))
            end
        end
        @test Array(D*transpose(Asym)) ≈ Array(D) * Array(transpose(Asym))
        @test Array(D*adjoint(Asym)) ≈ Array(D) * Array(adjoint(Asym))
        @test Array(D*transpose(Aherm)) ≈ Array(D) * Array(transpose(Aherm))
        @test Array(D*adjoint(Aherm)) ≈ Array(D) * Array(adjoint(Aherm))
        @test Array(transpose(Asym)*transpose(D)) ≈ Array(transpose(Asym)) * Array(transpose(D))
        @test Array(transpose(D)*transpose(Asym)) ≈ Array(transpose(D)) * Array(transpose(Asym))
        @test Array(adjoint(Aherm)*adjoint(D)) ≈ Array(adjoint(Aherm)) * Array(adjoint(D))
        @test Array(adjoint(D)*adjoint(Aherm)) ≈ Array(adjoint(D)) * Array(adjoint(Aherm))

        # Performance specialisations for A*_mul_B!
        vvv = similar(vv)
        @test (r = Matrix(D) * vv   ; mul!(vvv, D, vv)  ≈ r ≈ vvv)
        @test (r = Matrix(D)' * vv  ; mul!(vvv, adjoint(D), vv) ≈ r ≈ vvv)
        @test (r = transpose(Matrix(D)) * vv ; mul!(vvv, transpose(D), vv) ≈ r ≈ vvv)

        UUU = similar(UU)
        for transformA in (identity, adjoint, transpose)
            for transformD in (identity, adjoint, transpose)
                @test mul!(UUU, transformA(UU), transformD(D)) ≈  transformA(UU) * Matrix(transformD(D))
                @test mul!(UUU, transformD(D), transformA(UU)) ≈  Matrix(transformD(D)) * transformA(UU)
            end
        end

        alpha = elty(randn())  # randn(elty) does not work with BigFloat
        beta = elty(randn())
        @test begin
            vvv = similar(vv)
            vvv .= randn(size(vvv))  # randn!(vvv) does not work with BigFloat
            r = alpha * Matrix(D) * vv + beta * vvv
            mul!(vvv, D, vv, alpha, beta)  ≈ r ≈ vvv
        end
        @test begin
            vvv = similar(vv)
            vvv .= randn(size(vvv))  # randn!(vvv) does not work with BigFloat
            r = alpha * Matrix(D)' * vv + beta * vvv
            mul!(vvv, adjoint(D), vv, alpha, beta) ≈ r ≈ vvv
        end
        @test begin
            vvv = similar(vv)
            vvv .= randn(size(vvv))  # randn!(vvv) does not work with BigFloat
            r = alpha * transpose(Matrix(D)) * vv + beta * vvv
            mul!(vvv, transpose(D), vv, alpha, beta) ≈ r ≈ vvv
        end

        @test begin
            UUU = similar(UU)
            UUU .= randn(size(UUU))  # randn!(UUU) does not work with BigFloat
            r = alpha * Matrix(D) * UU + beta * UUU
            mul!(UUU, D, UU, alpha, beta) ≈ r ≈ UUU
        end
        @test begin
            UUU = similar(UU)
            UUU .= randn(size(UUU))  # randn!(UUU) does not work with BigFloat
            r = alpha * Matrix(D)' * UU + beta * UUU
            mul!(UUU, adjoint(D), UU, alpha, beta) ≈ r ≈ UUU
        end
        @test begin
            UUU = similar(UU)
            UUU .= randn(size(UUU))  # randn!(UUU) does not work with BigFloat
            r = alpha * transpose(Matrix(D)) * UU + beta * UUU
            mul!(UUU, transpose(D), UU, alpha, beta) ≈ r ≈ UUU
        end

        # make sure that mul!(A, {Adj|Trans}(B)) works with B as a Diagonal
        VV = Array(D)
        DD = copy(D)
        r  = VV * Matrix(D)
        @test Array(rmul!(VV, DD)) ≈ r ≈ Array(D)*Array(D)
        DD = copy(D)
        r  = VV * transpose(Array(D))
        @test Array(rmul!(VV, transpose(DD))) ≈ r
        DD = copy(D)
        r  = VV * Array(D)'
        @test Array(rmul!(VV, adjoint(DD))) ≈ r

        # kron
        D3 = Diagonal(convert(Vector{elty}, rand(n÷2)))
        DM3= Matrix(D3)
        @test Matrix(kron(D, D3)) ≈ kron(DM, DM3)
        M4 = rand(elty, n÷2, n÷2)
        @test kron(D3, M4) ≈ kron(DM3, M4)
        @test kron(M4, D3) ≈ kron(M4, DM3)
        X = [ones(1,1) for i in 1:2, j in 1:2]
        @test kron(I(2), X)[1,3] == zeros(1,1)
        X = [ones(2,2) for i in 1:2, j in 1:2]
        @test kron(I(2), X)[1,3] == zeros(2,2)
    end
    @testset "iszero, isone, triu, tril" begin
        Dzero = Diagonal(zeros(elty, 10))
        Done = Diagonal(ones(elty, 10))
        Dmix = Diagonal(zeros(elty, 10))
        Dmix[end,end] = one(elty)
        @test iszero(Dzero)
        @test !isone(Dzero)
        @test !iszero(Done)
        @test isone(Done)
        @test !iszero(Dmix)
        @test !isone(Dmix)
        @test istriu(D)
        @test istril(D)
        @test iszero(triu(D,1))
        @test triu(D,0)  == D
        @test triu(D,-1) == D
        @test tril(D,1)  == D
        @test iszero(tril(D,-1))
        @test tril(D,0)  == D
        @test_throws ArgumentError tril(D, -n - 2)
        @test_throws ArgumentError tril(D, n)
        @test_throws ArgumentError triu(D, -n)
        @test_throws ArgumentError triu(D, n + 2)
    end

    # factorize
    @test factorize(D) == D

    @testset "Eigensystem" begin
        eigD = eigen(D)
        @test Diagonal(eigD.values) ≈ D
        @test eigD.vectors == Matrix(I, size(D))
    end

    @testset "ldiv" begin
        v = rand(n + 1)
        @test_throws DimensionMismatch D\v
        v = rand(n)
        @test D\v ≈ DM\v
        V = rand(n + 1, n)
        @test_throws DimensionMismatch D\V
        V = rand(n, n)
        @test D\V ≈ DM\V
    end

    @testset "conj and transpose" begin
        @test transpose(D) == D
        if elty <: BlasComplex
            @test Array(conj(D)) ≈ conj(DM)
            @test adjoint(D) == conj(D)
        end
        # Translates to Ac/t_mul_B, which is specialized after issue 21286
        @test(D' * vv == conj(D) * vv)
        @test(transpose(D) * vv == D * vv)
    end

    # logdet and logabsdet
    if relty <: Real
        lD = Diagonal(convert(Vector{relty}, rand(n)))
        lM = Matrix(lD)
        @test logdet(lD) ≈ logdet(lM)
        d1, s1 = @inferred logabsdet(lD)
        d2, s2 = logabsdet(lM)
        @test d1 ≈ d2
        @test s1 == s2
        @test logdet(Diagonal(relty[-1,-2])) ≈ log(2)
        @test_throws DomainError logdet(Diagonal(relty[-1,-2,-3]))
    end

    @testset "similar" begin
        @test isa(similar(D), Diagonal{elty})
        @test isa(similar(D, Int), Diagonal{Int})
        @test isa(similar(D, (3,2)), Matrix{elty})
        @test isa(similar(D, Int, (3,2)), Matrix{Int})
    end

    # Issue number 10036
    # make sure issymmetric/ishermitian work for
    # non-real diagonal matrices
    @testset "issymmetric/hermitian for complex Diagonal" begin
        @test issymmetric(D2)
        @test ishermitian(D2)
        if elty <: Complex
            dc = d .+ elty(1im)
            D3 = Diagonal(dc)
            @test issymmetric(D3)
            @test !ishermitian(D3)
        end
    end

    @testset "svd (#11120/#11247)" begin
        U, s, V = svd(D)
        @test (U*Diagonal(s))*V' ≈ D
        @test svdvals(D) == s
        @test svd(D).V == V
    end

end

@testset "rdiv! (#40887)" begin
    @test rdiv!(Matrix(Diagonal([2.0, 3.0])), Diagonal(2:3)) == Diagonal([1.0, 1.0])
    @test rdiv!(fill(3.0, 3, 3), 3.0I(3)) == ones(3,3)
end

@testset "kron (issue #40595)" begin
    # custom array type to test that kron on Diagonal matrices preserves types of the parents if possible
    struct KronTestArray{T, N, AT} <: AbstractArray{T, N}
        data::AT
    end
    KronTestArray(data::AbstractArray) = KronTestArray{eltype(data), ndims(data), typeof(data)}(data)
    Base.size(A::KronTestArray) = size(A.data)
    LinearAlgebra.kron(A::KronTestArray, B::KronTestArray) = KronTestArray(kron(A.data, B.data))
    Base.getindex(K::KronTestArray{<:Any,N}, i::Vararg{Int,N}) where {N} = K.data[i...]

    A = KronTestArray([1, 2, 3]);
    @test kron(A, A) isa KronTestArray
    Ad = Diagonal(A);
    @test kron(Ad, Ad).diag isa KronTestArray
    @test kron(Ad, Ad).diag == kron([1, 2, 3], [1, 2, 3])
end

@testset "svdvals and eigvals (#11120/#11247)" begin
    D = Diagonal(Matrix{Float64}[randn(3,3), randn(2,2)])
    @test sort([svdvals(D)...;], rev = true) ≈ svdvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]])
    @test sort([eigvals(D)...;], by=LinearAlgebra.eigsortby) ≈ eigvals([D.diag[1] zeros(3,2); zeros(2,3) D.diag[2]])
end

@testset "eigmin (#27847)" begin
    for _ in 1:100
        d = randn(rand(1:10))
        D = Diagonal(d)
        @test eigmin(D) == minimum(d)
    end
end

@testset "isposdef" begin
    @test isposdef(Diagonal(1.0 .+ rand(n)))
    @test !isposdef(Diagonal(-1.0 * rand(n)))
    @test isposdef(Diagonal(complex(1.0, 0.0) .+ rand(n)))
    @test !isposdef(Diagonal(complex(1.0, 1.0) .+ rand(n)))
    @test isposdef(Diagonal([[1 0; 0 1], [1 0; 0 1]]))
    @test !isposdef(Diagonal([[1 0; 0 1], [1 0; 1 1]]))
end

@testset "getindex" begin
    d = randn(n)
    D = Diagonal(d)
    # getindex bounds checking
    @test_throws BoundsError D[0, 0]
    @test_throws BoundsError D[-1, -2]
    @test_throws BoundsError D[n, n + 1]
    @test_throws BoundsError D[n + 1, n]
    @test_throws BoundsError D[n + 1, n + 1]
    # getindex on and off the diagonal
    for i in 1:n, j in 1:n
        @test D[i, j] == (i == j ? d[i] : 0)
    end
end

@testset "setindex!" begin
    d = randn(n)
    D = Diagonal(d)
    # setindex! bounds checking
    @test_throws BoundsError D[0, 0] = 0
    @test_throws BoundsError D[-1 , -2] = 0
    @test_throws BoundsError D[n, n + 1] = 0
    @test_throws BoundsError D[n + 1, n] = 0
    @test_throws BoundsError D[n + 1, n + 1] = 0
    for i in 1:n, j in 1:n
        if i == j
            # setindex on! the diagonal
            @test ((D[i, j] = i) == i; D[i, j] == i)
        else
            # setindex! off the diagonal
            @test ((D[i, j] = 0) == 0; iszero(D[i, j]))
            @test_throws ArgumentError D[i, j] = 1
        end
    end
end

@testset "inverse" begin
    for d in (randn(n), [1, 2, 3], [1im, 2im, 3im])
        D = Diagonal(d)
        @test inv(D) ≈ inv(Array(D))
    end
    @test_throws SingularException inv(Diagonal(zeros(n)))
    @test_throws SingularException inv(Diagonal([0, 1, 2]))
    @test_throws SingularException inv(Diagonal([0im, 1im, 2im]))
end

# allow construct from range
@test all(Diagonal(range(1, stop=3, length=3)) .== Diagonal([1.0,2.0,3.0]))

# Issue 12803
for t in (Float32, Float64, Int, ComplexF64, Rational{Int})
    @test Diagonal(Matrix{t}[fill(t(1), 2, 2), fill(t(1), 3, 3)])[2,1] == zeros(t, 3, 2)
end

# Issue 15401
@test Matrix(1.0I, 5, 5) \ Diagonal(fill(1.,5)) == Matrix(I, 5, 5)

@testset "Triangular and Diagonal" begin
    function _test_matrix(type)
        if type == Int
            return rand(1:9, 5, 5)
        else
            return randn(type, 5, 5)
        end
    end
    types = (Float64, Int, ComplexF64)
    for ta in types
        D = Diagonal(_test_matrix(ta))
        for tb in types
            B = _test_matrix(tb)
            Tmats = (LowerTriangular(B), UnitLowerTriangular(B), UpperTriangular(B), UnitUpperTriangular(B))
            restypes = (LowerTriangular, LowerTriangular, UpperTriangular, UpperTriangular)
            for (T, rtype) in zip(Tmats, restypes)
                adjtype = (rtype == LowerTriangular) ? UpperTriangular : LowerTriangular

                # Triangular * Diagonal
                R = T * D
                @test R ≈ Array(T) * Array(D)
                @test isa(R, rtype)

                # Diagonal * Triangular
                R = D * T
                @test R ≈ Array(D) * Array(T)
                @test isa(R, rtype)

                # Adjoint of Triangular * Diagonal
                R = T' * D
                @test R ≈ Array(T)' * Array(D)
                @test isa(R, adjtype)

                # Diagonal * Adjoint of Triangular
                R = D * T'
                @test R ≈ Array(D) * Array(T)'
                @test isa(R, adjtype)

                # Transpose of Triangular * Diagonal
                R = transpose(T) * D
                @test R ≈ transpose(Array(T)) * Array(D)
                @test isa(R, adjtype)

                # Diagonal * Transpose of Triangular
                R = D * transpose(T)
                @test R ≈ Array(D) * transpose(Array(T))
                @test isa(R, adjtype)
            end
        end
    end
end

let D1 = Diagonal(rand(5)), D2 = Diagonal(rand(5))
    @test LinearAlgebra.rmul!(copy(D1),D2) == D1*D2
    @test LinearAlgebra.lmul!(D1,copy(D2)) == D1*D2
    @test LinearAlgebra.rmul!(copy(D1),transpose(D2)) == D1*transpose(D2)
    @test LinearAlgebra.lmul!(transpose(D1),copy(D2)) == transpose(D1)*D2
    @test LinearAlgebra.rmul!(copy(D1),adjoint(D2)) == D1*adjoint(D2)
    @test LinearAlgebra.lmul!(adjoint(D1),copy(D2)) == adjoint(D1)*D2
end

@testset "multiplication of a Diagonal with a Matrix" begin
    A = collect(reshape(1:8, 4, 2));
    B = BigFloat.(A);
    DL = Diagonal(collect(axes(A, 1)));
    DR = Diagonal(Float16.(collect(axes(A, 2))));

    @test DL * A == collect(DL) * A
    @test A * DR == A * collect(DR)
    @test DL * B == collect(DL) * B
    @test B * DR == B * collect(DR)

    A = reshape([ones(2,2), ones(2,2)*2, ones(2,2)*3, ones(2,2)*4], 2, 2)
    Ac = collect(A)
    D = Diagonal([collect(reshape(1:4, 2, 2)), collect(reshape(5:8, 2, 2))])
    Dc = collect(D)
    @test A * D == Ac * Dc
    @test D * A == Dc * Ac
    @test D * D == Dc * Dc

    AS = similar(A)
    mul!(AS, A, D, true, false)
    @test AS == A * D

    D2 = similar(D)
    mul!(D2, D, D)
    @test D2 == D * D

    copyto!(D2, D)
    lmul!(D, D2)
    @test D2 == D * D
    copyto!(D2, D)
    rmul!(D2, D)
    @test D2 == D * D
end

@testset "multiplication of QR Q-factor and Diagonal (#16615 spot test)" begin
    D = Diagonal(randn(5))
    Q = qr(randn(5, 5)).Q
    @test D * Q' == Array(D) * Q'
    Q = qr(randn(5, 5), ColumnNorm()).Q
    @test_throws ArgumentError lmul!(Q, D)
end

@testset "block diagonal matrices" begin
    D = Diagonal([[1 2; 3 4], [1 2; 3 4]])
    Dherm = Diagonal([[1 1+im; 1-im 1], [1 1+im; 1-im 1]])
    Dsym = Diagonal([[1 1+im; 1+im 1], [1 1+im; 1+im 1]])
    @test adjoint(D) == Diagonal([[1 3; 2 4], [1 3; 2 4]])
    @test transpose(D) == Diagonal([[1 3; 2 4], [1 3; 2 4]])
    @test adjoint(Dherm) == Dherm
    @test transpose(Dherm) == Diagonal([[1 1-im; 1+im 1], [1 1-im; 1+im 1]])
    @test adjoint(Dsym) == Diagonal([[1 1-im; 1-im 1], [1 1-im; 1-im 1]])
    @test transpose(Dsym) == Dsym

    v = [[1, 2], [3, 4]]
    @test Dherm' * v == Dherm * v
    @test transpose(D) * v == [[7, 10], [15, 22]]

    @test issymmetric(D) == false
    @test issymmetric(Dherm) == false
    @test issymmetric(Dsym) == true

    @test ishermitian(D) == false
    @test ishermitian(Dherm) == true
    @test ishermitian(Dsym) == false

    @test exp(D) == Diagonal([exp([1 2; 3 4]), exp([1 2; 3 4])])
    @test cis(D) == Diagonal([cis([1 2; 3 4]), cis([1 2; 3 4])])
    @test log(D) == Diagonal([log([1 2; 3 4]), log([1 2; 3 4])])
    @test sqrt(D) == Diagonal([sqrt([1 2; 3 4]), sqrt([1 2; 3 4])])

    @test tr(D) == 10
    @test det(D) == 4
end

@testset "linear solve for block diagonal matrices" begin
    D = Diagonal([rand(2,2) for _ in 1:5])
    b = [rand(2,2) for _ in 1:5]
    B = [rand(2,2) for _ in 1:5, _ in 1:5]
    @test ldiv!(D, copy(b)) ≈ Diagonal(inv.(D.diag)) * b
    @test ldiv!(D, copy(B)) ≈ Diagonal(inv.(D.diag)) * B
    @test rdiv!(copy(B), D) ≈ B * Diagonal(inv.(D.diag))
end

@testset "multiplication with Symmetric/Hermitian" begin
    for T in (Float64, ComplexF64)
        D = Diagonal(randn(T, n))
        A = randn(T, n, n); A = A'A
        S = Symmetric(A)
        H = Hermitian(A)
        for (transform1, transform2) in ((identity,  identity),
                (identity,  adjoint  ), (adjoint,   identity ), (adjoint,   adjoint  ),
                (identity,  transpose), (transpose, identity ), (transpose, transpose) )
            @test *(transform1(D), transform2(S)) ≈ *(transform1(Matrix(D)), transform2(Matrix(S)))
            @test *(transform1(D), transform2(H)) ≈ *(transform1(Matrix(D)), transform2(Matrix(H)))
            @test *(transform1(S), transform2(D)) ≈ *(transform1(Matrix(S)), transform2(Matrix(D)))
            @test *(transform1(S), transform2(H)) ≈ *(transform1(Matrix(S)), transform2(Matrix(H)))
        end
    end
end

@testset "multiplication of transposes of Diagonal (#22428)" begin
    for T in (Float64, ComplexF64)
        D = Diagonal(randn(T, 5, 5))
        B = Diagonal(randn(T, 5, 5))
        DD = Diagonal([randn(T, 2, 2), rand(T, 2, 2)])
        BB = Diagonal([randn(T, 2, 2), rand(T, 2, 2)])
        fullDD = copyto!(Matrix{Matrix{T}}(undef, 2, 2), DD)
        fullBB = copyto!(Matrix{Matrix{T}}(undef, 2, 2), BB)
        for (transform1, transform2) in ((identity,  identity),
                (identity,  adjoint  ), (adjoint,   identity ), (adjoint,   adjoint  ),
                (identity,  transpose), (transpose, identity ), (transpose, transpose))
            @test *(transform1(D), transform2(B))::typeof(D) ≈ *(transform1(Matrix(D)), transform2(Matrix(B))) atol=2 * eps()
            @test *(transform1(DD), transform2(BB))::typeof(DD) == *(transform1(fullDD), transform2(fullBB))
        end
        M = randn(T, 5, 5)
        MM = [randn(T, 2, 2) for _ in 1:2, _ in 1:2]
        for transform in (identity, adjoint, transpose)
            @test lmul!(transform(D), copy(M)) ≈ *(transform(Matrix(D)), M)
            @test rmul!(copy(M), transform(D)) ≈ *(M, transform(Matrix(D)))
            @test lmul!(transform(DD), copy(MM)) ≈ *(transform(fullDD), MM)
            @test rmul!(copy(MM), transform(DD)) ≈ *(MM, transform(fullDD))
        end
    end
end

@testset "Diagonal of adjoint/transpose vectors (#23649)" begin
    @test Diagonal(adjoint([1, 2, 3])) == Diagonal([1 2 3])
    @test Diagonal(transpose([1, 2, 3])) == Diagonal([1 2 3])
end

@testset "Multiplication with adjoint and transpose vectors (#26863)" begin
    x = collect(1:2)
    xt = transpose(x)
    A = reshape([[1 2; 3 4], zeros(Int,2,2), zeros(Int, 2, 2), [5 6; 7 8]], 2, 2)
    D = Diagonal(A)
    @test x'*D == x'*A == collect(x')*D == collect(x')*A
    @test xt*D == xt*A == collect(xt)*D == collect(xt)*A
    outadjxD = similar(x'*D); outtrxD = similar(xt*D);
    mul!(outadjxD, x', D)
    @test outadjxD == x'*D
    mul!(outtrxD, xt, D)
    @test outtrxD == xt*D

    D1 = Diagonal([[1 2; 3 4]])
    @test D1 * x' == D1 * collect(x') == collect(D1) * collect(x')
    @test D1 * xt == D1 * collect(xt) == collect(D1) * collect(xt)
    outD1adjx = similar(D1 * x'); outD1trx = similar(D1 * xt);
    mul!(outadjxD, D1, x')
    @test outadjxD == D1*x'
    mul!(outtrxD, D1, xt)
    @test outtrxD == D1*xt

    y = [x, x]
    yt = transpose(y)
    @test y'*D*y == (y'*D)*y == (y'*A)*y
    @test yt*D*y == (yt*D)*y == (yt*A)*y
    outadjyD = similar(y'*D); outtryD = similar(yt*D);
    outadjyD2 = similar(collect(y'*D)); outtryD2 = similar(collect(yt*D));
    mul!(outadjyD, y', D)
    mul!(outadjyD2, y', D)
    @test outadjyD == outadjyD2 == y'*D
    mul!(outtryD, yt, D)
    mul!(outtryD2, yt, D)
    @test outtryD == outtryD2 == yt*D
end

@testset "Multiplication of single element Diagonal (#36746, #40726)" begin
    @test_throws DimensionMismatch Diagonal(randn(1)) * randn(5)
    @test_throws DimensionMismatch Diagonal(randn(1)) * Diagonal(randn(3, 3))
    A = [1 0; 0 2]
    v = [3, 4]
    @test Diagonal(A) * v == A * v
    @test Diagonal(A) * Diagonal(A) == A * A
    @test_throws DimensionMismatch [1 0;0 1] * Diagonal([2 3])   # Issue #40726
    @test_throws DimensionMismatch lmul!(Diagonal([1]), [1,2,3]) # nearby
end

@testset "Triangular division by Diagonal #27989" begin
    K = 5
    for elty in (Float32, Float64, ComplexF32, ComplexF64)
        U = UpperTriangular(randn(elty, K, K))
        L = LowerTriangular(randn(elty, K, K))
        D = Diagonal(randn(elty, K))
        @test (U / D)::UpperTriangular{elty} == UpperTriangular(Matrix(U) / Matrix(D))
        @test (L / D)::LowerTriangular{elty} == LowerTriangular(Matrix(L) / Matrix(D))
        @test (D \ U)::UpperTriangular{elty} == UpperTriangular(Matrix(D) \ Matrix(U))
        @test (D \ L)::LowerTriangular{elty} == LowerTriangular(Matrix(D) \ Matrix(L))
    end
end

@testset "(Sym)Tridiagonal division by Diagonal" begin
    for K in (5, 1), elty in (Float64, ComplexF32), overlength in (1, 0)
        S = SymTridiagonal(randn(elty, K), randn(elty, K-overlength))
        T = Tridiagonal(randn(elty, K-1), randn(elty, K), randn(elty, K-1))
        D = Diagonal(randn(elty, K))
        D0 = Diagonal(zeros(elty, K))
        @test (D \ S)::Tridiagonal{elty} == Tridiagonal(Matrix(D) \ Matrix(S))
        @test (D \ T)::Tridiagonal{elty} == Tridiagonal(Matrix(D) \ Matrix(T))
        @test (S / D)::Tridiagonal{elty} == Tridiagonal(Matrix(S) / Matrix(D))
        @test (T / D)::Tridiagonal{elty} == Tridiagonal(Matrix(T) / Matrix(D))
        @test_throws SingularException D0 \ S
        @test_throws SingularException D0 \ T
        @test_throws SingularException S / D0
        @test_throws SingularException T / D0
    end
    # 0-length case
    S = SymTridiagonal(Float64[], Float64[])
    T = Tridiagonal(Float64[], Float64[], Float64[])
    D = Diagonal(Float64[])
    @test (D \ S)::Tridiagonal{Float64} == T
    @test (D \ T)::Tridiagonal{Float64} == T
    @test (S / D)::Tridiagonal{Float64} == T
    @test (T / D)::Tridiagonal{Float64} == T
    # matrix eltype case
    K = 5
    for elty in (Float64, ComplexF32), overlength in (1, 0)
        S = SymTridiagonal([rand(elty, 2, 2) for _ in 1:K], [rand(elty, 2, 2) for _ in 1:K-overlength])
        T = Tridiagonal([rand(elty, 2, 2) for _ in 1:K-1], [rand(elty, 2, 2) for _ in 1:K], [rand(elty, 2, 2) for _ in 1:K-1])
        D = Diagonal(randn(elty, K))
        SM = fill(zeros(elty, 2, 2), K, K)
        TM = copy(SM)
        SM[1,1] = S[1,1]; TM[1,1] = T[1,1]
        for j in 2:K
            SM[j,j-1] = S[j,j-1]; SM[j,j] = S[j,j]; SM[j-1,j] = S[j-1,j]
            TM[j,j-1] = T[j,j-1]; TM[j,j] = T[j,j]; TM[j-1,j] = T[j-1,j]
        end
        for (M, Mm) in ((S, SM), (T, TM))
            DS = D \ M
            @test DS isa Tridiagonal
            DM = D \ Mm
            for i in -1:1; @test diag(DS, i) ≈ diag(DM, i) end
        end
    end
    # eltype promotion case
    S = SymTridiagonal(rand(-20:20, K), rand(-20:20, K-1))
    T = Tridiagonal(rand(-20:20, K-1), rand(-20:20, K), rand(-20:20, K-1))
    D = Diagonal(rand(1:20, K))
    @test (D \ S)::Tridiagonal{Float64} == Tridiagonal(Matrix(D) \ Matrix(S))
    @test (D \ T)::Tridiagonal{Float64} == Tridiagonal(Matrix(D) \ Matrix(T))
    @test (S / D)::Tridiagonal{Float64} == Tridiagonal(Matrix(S) / Matrix(D))
    @test (T / D)::Tridiagonal{Float64} == Tridiagonal(Matrix(T) / Matrix(D))
end

@testset "eigenvalue sorting" begin
    D = Diagonal([0.4, 0.2, -1.3])
    @test eigvals(D) == eigen(D).values == [0.4, 0.2, -1.3] # not sorted by default
    @test eigvals(Matrix(D)) == eigen(Matrix(D)).values == [-1.3, 0.2, 0.4] # sorted even if diagonal special case is detected
    E = eigen(D, sortby=abs) # sortby keyword supported for eigen(::Diagonal)
    @test E.values == [0.2, 0.4, -1.3]
    @test E.vectors == [0 1 0; 1 0 0; 0 0 1]
end

@testset "sum, mapreduce" begin
    D = Diagonal([1,2,3])
    Ddense = Matrix(D)
    @test sum(D) == 6
    @test_throws ArgumentError sum(D, dims=0)
    @test sum(D, dims=1) == sum(Ddense, dims=1)
    @test sum(D, dims=2) == sum(Ddense, dims=2)
    @test sum(D, dims=3) == sum(Ddense, dims=3)
    @test typeof(sum(D, dims=1)) == typeof(sum(Ddense, dims=1))
    @test mapreduce(one, min, D, dims=1) == mapreduce(one, min, Ddense, dims=1)
    @test mapreduce(one, min, D, dims=2) == mapreduce(one, min, Ddense, dims=2)
    @test mapreduce(one, min, D, dims=3) == mapreduce(one, min, Ddense, dims=3)
    @test typeof(mapreduce(one, min, D, dims=1)) == typeof(mapreduce(one, min, Ddense, dims=1))
    @test mapreduce(zero, max, D, dims=1) == mapreduce(zero, max, Ddense, dims=1)
    @test mapreduce(zero, max, D, dims=2) == mapreduce(zero, max, Ddense, dims=2)
    @test mapreduce(zero, max, D, dims=3) == mapreduce(zero, max, Ddense, dims=3)
    @test typeof(mapreduce(zero, max, D, dims=1)) == typeof(mapreduce(zero, max, Ddense, dims=1))

    D = Diagonal(Int[])
    Ddense = Matrix(D)
    @test sum(D) == 0
    @test_throws ArgumentError sum(D, dims=0)
    @test sum(D, dims=1) == sum(Ddense, dims=1)
    @test sum(D, dims=2) == sum(Ddense, dims=2)
    @test sum(D, dims=3) == sum(Ddense, dims=3)
    @test typeof(sum(D, dims=1)) == typeof(sum(Ddense, dims=1))

    D = Diagonal(Int[2])
    Ddense = Matrix(D)
    @test sum(D) == 2
    @test_throws ArgumentError sum(D, dims=0)
    @test sum(D, dims=1) == sum(Ddense, dims=1)
    @test sum(D, dims=2) == sum(Ddense, dims=2)
    @test sum(D, dims=3) == sum(Ddense, dims=3)
    @test typeof(sum(D, dims=1)) == typeof(sum(Ddense, dims=1))
end

@testset "logabsdet for generic eltype" begin
    d = Any[1, -2.0, -3.0]
    D = Diagonal(d)
    d1, s1 = logabsdet(D)
    @test d1 ≈ sum(log ∘ abs, d)
    @test s1 == prod(sign, d)
end

@testset "Empty (#35424)" begin
    @test zeros(0)'*Diagonal(zeros(0))*zeros(0) === 0.0
    @test transpose(zeros(0))*Diagonal(zeros(Complex{Int}, 0))*zeros(0) === 0.0 + 0.0im
    @test dot(zeros(Int32, 0), Diagonal(zeros(Int, 0)), zeros(Int16, 0)) === 0
end

@testset "Diagonal(undef)" begin
    d = Diagonal{Float32}(undef, 2)
    @test length(d.diag) == 2
end

@testset "permutedims (#39447)" begin
    for D in (Diagonal(zeros(5)), Diagonal(zeros(5) .+ 1im), Diagonal([[1,2],[3,4]]))
        @test permutedims(D) === permutedims(D,(1,2)) === permutedims(D,(2,1)) === D
        @test_throws ArgumentError permutedims(D,(1,3))
    end
end

@testset "Inner product" begin
    A = Diagonal(rand(10) .+ im)
    B = Diagonal(rand(10) .+ im)
    @test dot(A, B) ≈ dot(Matrix(A), B)
    @test dot(A, B) ≈ dot(A, Matrix(B))
    @test dot(A, B) ≈ dot(Matrix(A), Matrix(B))
    @test dot(A, B) ≈ conj(dot(B, A))
end

@testset "eltype relaxation(#41015)" begin
    A = rand(3,3)
    for trans in (identity, adjoint, transpose)
        @test ldiv!(trans(I(3)), A) == A
        @test rdiv!(A, trans(I(3))) == A
    end
end

const BASE_TEST_PATH = joinpath(Sys.BINDIR, "..", "share", "julia", "test")
isdefined(Main, :ImmutableArrays) || @eval Main include(joinpath($(BASE_TEST_PATH), "testhelpers", "ImmutableArrays.jl"))
using .Main.ImmutableArrays

@testset "Conversion to AbstractArray" begin
    # tests corresponding to #34995
    d = ImmutableArray([1, 2, 3, 4])
    D = Diagonal(d)

    @test convert(AbstractArray{Float64}, D)::Diagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == D
    @test convert(AbstractMatrix{Float64}, D)::Diagonal{Float64,ImmutableArray{Float64,1,Array{Float64,1}}} == D
end

@testset "divisions functionality" for elty in (Int, Float64, ComplexF64)
    B = Diagonal(rand(elty,5,5))
    x = rand(elty)
    @test \(x, B) == /(B, x)
end

@testset "promotion" begin
    for (v1, v2) in (([true], [1]), ([zeros(2,2)], [zeros(Int, 2,2)]))
        T = promote_type(eltype(v1), eltype(v2))
        V = promote_type(typeof(v1), typeof(v2))
        d1 = Diagonal(v1)
        d2 = Diagonal(v2)
        v = [d1, d2]
        @test (@inferred eltype(v)) == Diagonal{T, V}
    end
    # test for a type for which promote_type doesn't lead to a concrete eltype
    struct MyArrayWrapper{T,N,A<:AbstractArray{T,N}} <: AbstractArray{T,N}
       a :: A
    end
    Base.size(M::MyArrayWrapper) = size(M.a)
    Base.axes(M::MyArrayWrapper) = axes(M.a)
    Base.length(M::MyArrayWrapper) = length(M.a)
    Base.getindex(M::MyArrayWrapper, i::Int...) = M.a[i...]
    Base.setindex!(M::MyArrayWrapper, v, i::Int...) = M.a[i...] = v
    d1 = Diagonal(MyArrayWrapper(1:3))
    d2 = Diagonal(MyArrayWrapper(1.0:3.0))
    c = [d1, d2]
    @test c[1] == d1
    @test c[2] == d2
end

@testset "zero and one" begin
    D1 = Diagonal(rand(3))
    @test D1 + zero(D1) == D1
    @test D1 * one(D1) == D1
    @test D1 * oneunit(D1) == D1
    @test oneunit(D1) isa typeof(D1)
    D2 = Diagonal([collect(reshape(1:4, 2, 2)), collect(reshape(5:8, 2, 2))])
    @test D2 + zero(D2) == D2
    @test D2 * one(D2) == D2
    @test D2 * oneunit(D2) == D2
    @test oneunit(D2) isa typeof(D2)
    D3 = Diagonal([D2, D2]);
    @test D3 + zero(D3) == D3
    @test D3 * one(D3) == D3
    @test D3 * oneunit(D3) == D3
    @test oneunit(D3) isa typeof(D3)
end

@testset "AbstractTriangular" for (Tri, UTri) in ((UpperTriangular, UnitUpperTriangular), (LowerTriangular, UnitLowerTriangular))
    A = randn(4, 4)
    TriA = Tri(A)
    UTriA = UTri(A)
    D = Diagonal(1.0:4.0)
    DM = Matrix(D)
    DMF = factorize(DM)
    outTri = similar(TriA)
    out = similar(A)
    # 2 args
    for fun in (*, rmul!, rdiv!, /)
        @test fun(copy(TriA), D)::Tri == fun(Matrix(TriA), D)
        @test fun(copy(UTriA), D)::Tri == fun(Matrix(UTriA), D)
    end
    for fun in (*, lmul!, ldiv!, \)
        @test fun(D, copy(TriA))::Tri == fun(D, Matrix(TriA))
        @test fun(D, copy(UTriA))::Tri == fun(D, Matrix(UTriA))
    end
    # 3 args
    @test outTri === ldiv!(outTri, D, TriA)::Tri == ldiv!(out, D, Matrix(TriA))
    @test outTri === ldiv!(outTri, D, UTriA)::Tri == ldiv!(out, D, Matrix(UTriA))
    @test outTri === mul!(outTri, D, TriA)::Tri == mul!(out, D, Matrix(TriA))
    @test outTri === mul!(outTri, D, UTriA)::Tri == mul!(out, D, Matrix(UTriA))
    @test outTri === mul!(outTri, TriA, D)::Tri == mul!(out, Matrix(TriA), D)
    @test outTri === mul!(outTri, UTriA, D)::Tri == mul!(out, Matrix(UTriA), D)
    # 5 args
    @test outTri === mul!(outTri, D, TriA, 2, 1)::Tri == mul!(out, D, Matrix(TriA), 2, 1)
    @test outTri === mul!(outTri, D, UTriA, 2, 1)::Tri == mul!(out, D, Matrix(UTriA), 2, 1)
    @test outTri === mul!(outTri, TriA, D, 2, 1)::Tri == mul!(out, Matrix(TriA), D, 2, 1)
    @test outTri === mul!(outTri, UTriA, D, 2, 1)::Tri == mul!(out, Matrix(UTriA), D, 2, 1)
end

end # module TestDiagonal

max_stars_repo_name: wentasah/julia

max_issues_repo_name: wentasah/julia

max_forks_repo_name: wentasah/julia


text:
# this is the default stack page
# ultimately this should be a property of either the CPU or Chipset
const STACK_PAGE = 0x0100

# TODO WTF is with the break (B) flag?
mutable struct FlagsRegister
    C::Bool
    Z::Bool
    I::Bool
    D::Bool
    B::Bool
    Ω::Bool
    V::Bool
    N::Bool

    FlagsRegister() = new(false, false, false, false, false, false, false, false)
end

mutable struct CPU
    A::UInt8  # accumulator register
    X::UInt8  # index registers
    Y::UInt8
    SP::UInt8  # stack pointer register, I'm tempted to make this 16-bit
    PC::UInt16  # program counter register
    flags::FlagsRegister

    function CPU(A::UInt8=0x00, X::UInt8=0x00, Y::UInt8=0x00,
                 SP::UInt8=0xff, PC::UInt16=0x0600, flags::FlagsRegister=FlagsRegister())
        new(A, X, Y, SP, PC, flags)
    end
end

@inline Base.fetch(c::CPU, reg::Symbol) = getfield(c, reg)
@inline store!(c::CPU, reg::Symbol, val::Integer) = setfield!(c, reg, val)

counter!(c::CPU, ℓ::Unsigned) = (c.PC += ℓ)

function describe(c::CPU)
    st = status_string(c)
    A = hexstring(c.A); X = hexstring(c.X); Y = hexstring(c.Y)
    SP = hexstring(c.SP); PC = hexstring(c.PC)
    """
    MOS6502 CPU
        Registers:
            A=$A  X=$X  Y=$Y
            SP=$SP  PC=$PC
        Status: NV-BDIZC
                $st
    """
end

Base.show(io::IO, c::CPU) = print(describe(c))

#==============================================================================================
    <status register access>
==============================================================================================#
status(fr::FlagsRegister, flag::Symbol)::Bool = getfield(fr, flag)
status(c::CPU, flag::Symbol)::Bool = getfield(fr, flag)

# these are toggles
status!(fr::FlagsRegister, flag::Symbol) = setfield!(fr, ~getfield(fr, flag))
status!(c::CPU, flag::Symbol) = status!(c.flags, flag)

status!(fr::FlagsRegister, flag::Symbol, val::Bool) = setfield!(fr, flag, val)
status!(c::CPU, flag::Symbol, val::Bool) = status!(c.flags, flag, val)

function status_string(fr::FlagsRegister)
    string((Int(getfield(fr, f)) for f ∈ reverse(fieldnames(typeof(fr))))...)
end
status_string(c::CPU) = status_string(c.flags)

#==============================================================================================
    </status register access>
==============================================================================================#

#==============================================================================================
    <special pointers>
==============================================================================================#
stackpointer(c::CPU) = Π(STACK_PAGE + c.SP)
#==============================================================================================
    </special pointers>
==============================================================================================#


text:
immutable State1D
    s::Float64 # position
    v::Float64 # speed [m/s]
end
Base.write(io::IO, ::MIME"text/plain", s::State1D) = @printf(io, "%.16e %.16e", s.s, s.v)
function Base.read(io::IO, ::MIME"text/plain", ::Type{State1D})
    i = 0
    tokens = split(strip(readline(io)), ' ')
    s = parse(Float64, tokens[i+=1])
    v = parse(Float64, tokens[i+=1])
    return State1D(s,v)
end


text:
"""
  set_simulate_empirical_moments!(sMMProblem::MSMProblem, f::Function)

Function to set the field simulate_empirical_moments for a MSMProblem.
The function simulate_empirical_moments takes parameter values and return
the corresponding simulate moments values.
"""
function set_simulate_empirical_moments!(sMMProblem::MSMProblem, f::Function)

  # set the function that returns an ordered dictionary
  sMMProblem.simulate_empirical_moments = f

  # set the function that returns an array (respecting the order of the ordered dict)
  # this function is used to calculate the jacobian
  function simulate_empirical_moments_array(x)

    momentsODict = sMMProblem.simulate_empirical_moments(x)
    momentsArray = Array{Float64}(undef, length(momentsODict))

    for (i, k) in enumerate(keys(momentsODict))
        momentsArray[i] = momentsODict[k]
    end

    return momentsArray

  end

  sMMProblem.simulate_empirical_moments_array = simulate_empirical_moments_array

end

"""
  construct_objective_function!(sMMProblem::MSMProblem)

Function that construct an objective function, using the function
MSMProblem.simulate_empirical_moments.
"""
function construct_objective_function!(sMMProblem::MSMProblem)

    function objective_function_MSM(x)

          # Initialization
          #----------------
          distanceEmpSimMoments = sMMProblem.options.penaltyValue

          #---------------------------------------------------------------------
          # A. Generate simulated moments
          #---------------------------------------------------------------------
          simulatedMoments, convergence = try

          sMMProblem.simulate_empirical_moments(x), 1

      catch errorSimulation

            info("An error occured with parameter values = $(x)")
            info("$(errorSimulation)")

            OrderedDict{String,Array{Float64,1}}(), 0

      end

      #------------------------------------------------------------------------
      # B. If generating moment was successful, calculate distance between empirical
      # and simulated moments
      # (If no convergence, returns penalty value : sMMProblem.options.penaltyValue)
      #------------------------------------------------------------------------
      if convergence == 1

        # to store the distance between empirical and simulated moments
        arrayDistance = zeros(length(keys(sMMProblem.empiricalMoments)))

        for (indexMoment, k) in enumerate(keys(sMMProblem.empiricalMoments))

          # * sMMProblem.empiricalMoments[k][1] is the empirical moments
          #---------------------------------------------------------------------
          arrayDistance[indexMoment] = (sMMProblem.empiricalMoments[k][1] - simulatedMoments[k])

        end

        # formula is (m - m*)'*W*(m - m*)'
        distanceEmpSimMoments = transpose(arrayDistance)*sMMProblem.W*arrayDistance

        if sMMProblem.options.showDistance == true
          println("distance = $(distanceEmpSimMoments)")
        end

      end

      return distanceEmpSimMoments

    end

    # Attach the objective function
    sMMProblem.objective_function = objective_function_MSM

end

"""
  set_priors!(sMMProblem::MSMProblem, priors::OrderedDict{String,Array{Float64,1}})

Function to change the field sMMProblem.priors
"""
function set_priors!(sMMProblem::MSMProblem, priors::OrderedDict{String,Array{Float64,1}})

  sMMProblem.priors = priors

end

"""
   set_weight_matrix!(sMMProblem::MSMProblem, W::Matrix{Float64})

Function to change the field sMMProblem.empiricalMoments
"""
function set_weight_matrix!(sMMProblem::MSMProblem, W::Matrix{Float64})

  sMMProblem.W = W

end

"""
   set_empirical_moments!(sMMProblem::MSMProblem, empiricalMoments::OrderedDict{String,Array{Float64,1}})

Function to change the field sMMProblem.empiricalMoments
"""
function set_empirical_moments!(sMMProblem::MSMProblem, empiricalMoments::OrderedDict{String,Array{Float64,1}})

  sMMProblem.empiricalMoments = empiricalMoments

end

"""
   set_Sigma0!(sMMProblem::MSMProblem, Sigma0::Array{Float64,2})

Function to change the field sMMProblem.Sigma0, where Sigma0 is the distance matrix,
in the terminology of Duffie and Singleton (1993)
"""
function  set_Sigma0!(sMMProblem::MSMProblem, Sigma0::Array{Float64,2})

  sMMProblem.Sigma0 = Sigma0

end

"""
  set_global_optimizer!(sMMProblem::MSMProblem)

Function to set the fields corresponding to the global
optimizer problem.
"""
function set_global_optimizer!(sMMProblem::MSMProblem)

  if is_bb_optimizer(sMMProblem.options.globalOptimizer) == true

    set_bbSetup!(sMMProblem)

  else

    Base.error("sMMProblem.options.globalOptimizer = $(sMMProblem.options.globalOptimizer) is not supported.")

  end

end

"""
  set_bbSetup!(sMMProblem::MSMProblem)

Function to set the field bbSetup for a MSMProblem.
"""
function set_bbSetup!(sMMProblem::MSMProblem)

  # A. using sMMProblem.priors, generate searchRange:
  #-------------------------------------------------
  mySearchRange = generate_bbSearchRange(sMMProblem)

  info("$(nworkers()) worker(s) detected")
   # Debug:
   #-------
   @sync for (idx, pid) in enumerate(workers())
     @async @spawnat(pid, println("hello"))
   end

  if nworkers() == 1
    info("Starting optimization in serial")
    sMMProblem.bbSetup = bbsetup(sMMProblem.objective_function;
                              Method = sMMProblem.options.globalOptimizer,
                              SearchRange = mySearchRange,
                              MaxFuncEvals = sMMProblem.options.maxFuncEvals,
                              TraceMode = :verbose,
                              PopulationSize = sMMProblem.options.populationSize,
                              NumDimensions = length(keys(sMMProblem.priors)))
  else
    info("Starting optimization in parallel")
    sMMProblem.bbSetup = bbsetup(sMMProblem.objective_function;
                                Method = sMMProblem.options.globalOptimizer,
                                SearchRange = mySearchRange,
                                MaxFuncEvals = sMMProblem.options.maxFuncEvals,
                                Workers = workers(),
                                PopulationSize = sMMProblem.options.populationSize,
                                TraceMode = :verbose,
                                NumDimensions = length(keys(sMMProblem.priors)))
  end

end

"""
  generate_bbSearchRange(sMMProblem::MSMProblem)

Function to generate a search range that matches the convention used by
BlackBoxOptim.
"""
function generate_bbSearchRange(sMMProblem::MSMProblem)

  # sMMProblem.priors["key"][1] contains the initial guess
  # sMMProblem.priors["key"][2] contains the lower bound
  # sMMProblem.priors["key"][3] contains the upper bound
  #-----------------------------------------------------
  [(sMMProblem.priors[k][2], sMMProblem.priors[k][3]) for k in keys(sMMProblem.priors)]
end

"""
  create_lower_bound(sMMProblem::MSMProblem)

Function to generate a lower bound used by Optim when minimizing with Fminbox.
The lower bound is of type Array{Float64,1}.
"""
function create_lower_bound(sMMProblem::MSMProblem)

  # sMMProblem.priors["key"][1] contains the initial guess
  # sMMProblem.priors["key"][2] contains the lower bound
  # sMMProblem.priors["key"][3] contains the upper bound
  #-----------------------------------------------------
  [sMMProblem.priors[k][2] for k in keys(sMMProblem.priors)]
end

"""
  create_upper_bound(sMMProblem::MSMProblem)

Function to generate a lower bound used by Optim when minimizing with Fminbox.
The upper bound is of type Array{Float64,1}.
"""
function create_upper_bound(sMMProblem::MSMProblem)

  # sMMProblem.priors["key"][1] contains the initial guess
  # sMMProblem.priors["key"][2] contains the lower bound
  # sMMProblem.priors["key"][3] contains the upper bound
  #-----------------------------------------------------
  [sMMProblem.priors[k][3] for k in keys(sMMProblem.priors)]
end

# Function coded by Robert Feldt. All credits to him. I only changed two lines:
# * cubedim = Vector{T}(n) instead of cubedim = Vector{T}(undef, n)
# * I return return transpose(result) instead of result
# source: https://github.com/robertfeldt/BlackBoxOptim.jl/blob/master/src/utilities/latin_hypercube_sampling.jl
"""
    latin_hypercube_sampling(mins, maxs, n)

Randomly sample `n` vectors from the parallelogram defined
by `mins` and `maxs` using the Latin hypercube algorithm.
Returns `dims`×`n` matrix.
"""
function latin_hypercube_sampling(mins::AbstractVector{T},
                                  maxs::AbstractVector{T},
                                  n::Integer) where T<:Number
    length(mins) == length(maxs) ||
        throw(DimensionMismatch("mins and maxs should have the same length"))
    all(xy -> xy[1] <= xy[2], zip(mins, maxs)) ||
        throw(ArgumentError("mins[i] should not exceed maxs[i]"))
    dims = length(mins)
    result = zeros(T, dims, n)
    # Julia 0.7
    #----------
    # cubedim = Vector{T}(undef, n)
    # Julia 0.6.4
    #------------
    cubedim = Vector{T}(undef, n)
    @inbounds for i in 1:dims
        imin = mins[i]
        dimstep = (maxs[i] - imin) / n
        for j in 1:n
            cubedim[j] = imin + dimstep * (j - 1 + rand(T))
        end
        result[i, :] .= shuffle!(cubedim)
    end
    return transpose(result)
end

"""
    latin_hypercube_sampling(mySearchRange::Vector{Tuple{Float64, Float64}}, n::Integer; gens::Integer = 100)

Function to create an optimised Latin Hypercube Sampling Plan
## Input
* mySearchRange: a vector of tuple of the form (xj_min, xj_max)
* n: number of points to draw
* gens: (optional) optimization is run for gens generations. See LatinHypercubeSampling.jl

## Output
* Matrix{Float64}: row = observation; column = dimension
"""
function latin_hypercube_sampling(mySearchRange::Vector{Tuple{Float64, Float64}}, n::Integer; gens::Integer = 100)

        d = length(mySearchRange) #dimension
        # See: https://mrurq.github.io/LatinHypercubeSampling.jl/stable/man/lhcoptim/
        plan, _ = LHCoptim(n,d,gens)
        # Rescale plan:
        scaled_plan = scaleLHC(plan, mySearchRange)

        return scaled_plan
end

function sobol_sampling(lb::Vector{Float64},
                        ub::Vector{Float64},
                        n::Integer)
    s = SobolSeq(lb,ub)
    # "If you know in advance the number n of points that you plan to generate,
    # some authors suggest that better uniformity can be attained by first skipping
    # the initial portion of the LDS.". See: https://github.com/stevengj/Sobol.jl
    skip(s,n)
    d = length(lb) #dimension
    if d == 1
        result = transpose(reduce(hcat, next!(s)[1] for i = 1:n))
    else
        result = transpose(reduce(hcat, next!(s) for i = 1:n))
    end
    return result
end

"""
  get_now()

Returns date and time in a manner that does not clash with Windows, Linux and OSX
"""
function get_now()

  "$(Dates.today())--$(Dates.hour(Dates.now()))h-$(Dates.minute(Dates.now()))m-$(Dates.second(Dates.now()))s"

end

"""
    info(text)

To display information to user
"""
function info(text)
    @info text
end

"""
  linspace(z_n::Int64, z_start::Real, z_end::Real)

Similar behavior of Base.linespace on julia v. < 0.6
"""
function linspace(z_start::Real, z_end::Real, z_n::Int64)
    return collect(range(z_start,stop=z_end,length=z_n))
end