text:
function test_interface(
    rng::AbstractRNG, lik, k::Kernel, x::AbstractVector; functor_args=()
)
    gp = GP(k)
    lgp = LatentGP(gp, lik, 1e-5)
    lfgp = lgp(x)

    # Check if likelihood produces a distribution
    @test lik(rand(rng, lfgp.fx)) isa Distribution

    N = length(x)
    y = rand(rng, lfgp.fx)

    if x isa MOInput
        # TODO: replace with mo_inverse_transform
        N = length(x.x)
        y = [y[[i + j * N for j in 0:(x.out_dim - 1)]] for i in 1:N]
    end

    # Check if the likelihood samples are of correct length
    @test length(rand(rng, lik(y))) == N

    # Check if functor works properly
    xs, re = Functors.functor(lik)
    @test lik == re(xs)
    if isempty(functor_args)
        @test xs === ()
    else
        @test keys(xs) == functor_args
    end

    return nothing
end

"""
    test_interface(lik, k::Kernel, x::AbstractVector; functor_args=())

This function provides unified method to check the interface of the various likelihoods
defined. It checks if the likelihood produces a distribution, length of likelihood
samples is correct and if the functor works as intended.
...
# Arguments
- `lik`: the likelihood to test the interface of
- `k::Kernel`: the kernel to use for the GP
- `x::AbstractVector`: intputs to compute the likelihood on
- `functor_args=()`: a collection of symbols of arguments to match functor parameters with.
...
"""
function test_interface(lik, k::KernelFunctions.Kernel, x::AbstractVector; kwargs...)
    return test_interface(Random.GLOBAL_RNG, lik, k, x; kwargs...)
end


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

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

# --- Imports

# Standard library
import InteractiveUtils.subtypes
import LinearAlgebra
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- File inclusions

# Test utilities
include("test_utils.jl")

# --- Tests

# ------ M::Multivector

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

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

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

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

# --- B::Blade

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

    test_dim = 20

    C = Blade(rand(test_dim, 7))

    # --- Tests

    # ------ dim(B) == dim(C), grade(B) < grade(C)

    for grade_B in 2:5
        B_vectors = basis(C)[:, 1:grade_B] * rand(grade_B, grade_B)
        B = Blade(B_vectors)

        B_dual_C = dual(B, C)

        # --- Check dim, grade, and norm

        @test dim(B_dual_C) == test_dim
        @test grade(B_dual_C) == grade(C) - grade_B
        @test norm(B_dual_C) == norm(B)

        # --- Check that B and dual(B) are orthogonal complements

        @test LinearAlgebra.norm(transpose(basis(B_dual_C)) * basis(B)) <
            10 * eps(Float64)

        # --- Check sign(dual(B, C))

        # Compute sign of I_C formed from basis(B) and basis(dual(B, C))
        sign_Q = sign(LinearAlgebra.det(
            transpose(hcat(basis(B), basis(B_dual_C))) * basis(C)))

        # Compute expected_sign
        expected_sign = sign(B) * sign_Q

        # Account for sign of I_C^{-1} relative to I_C
        if mod(grade(C), 4) >= 2
            expected_sign = -expected_sign
        end

        # Account for reversals required to eliminate B
        if mod(grade(B), 4) >= 2
            expected_sign = -expected_sign
        end

        @test sign(B_dual_C) == expected_sign

        # Check dual(dual(B, C), C) = (-1)^(grade(C) * (grade(C) - 1) / 2) B
        if mod(grade(C), 4) < 2
            @test dual(B_dual_C, C) ≈ B
        else
            @test dual(B_dual_C, C) ≈ -B
        end
    end

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

    coefficients = rand(grade(C), grade(C))
    B_vectors = basis(C) * coefficients
    B = Blade(B_vectors)

    relative_orientation =
        sign(LinearAlgebra.det(transpose(basis(B)) * basis(C)))
    expected_result = mod(grade(B), 4) < 2 ?
        Scalar(relative_orientation * volume(B)) :
        Scalar(-relative_orientation * volume(B))

    @test dual(B, C) ≈ expected_result

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

    B = Blade(rand(test_dim + 1, grade(C) - 1))
    @test_throws DimensionMismatch dual(B, C)

    # ------ B not contained in C

    for grade_B in 2:5
        # grade(B) > grade(C)
        B = Blade(rand(test_dim, grade(C) + 1))
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)

        # grade(B) <= grade(C), B not contained in C
        B = Blade(rand(test_dim, grade_B))
        while LinearAlgebra.norm(reject(basis(B), C)) < sqrt(eps(Float64))
            B = Blade(rand(test_dim, grade_B))
        end
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

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

    test_dim = 10

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

    # --- Tests

    # dim(B) == dim(C)
    for test_grade in 2:5
        B = Blade(rand(test_dim, test_grade))
        C = Pseudoscalar(test_dim, test_value)

        B_dual_C = dual(B, C)
        @test B_dual_C == dual(B)

        # Check dual(dual(B, C), C) = (-1)^(grade(C) * (grade(C) - 1) / 2) B
        if mod(grade(C), 4) < 2
            @test dual(B_dual_C, C) ≈ B
        else
            @test dual(B_dual_C, C) ≈ -B
        end
    end

    # dim(B) != dim(C)
    B = Blade(rand(test_dim, 3))
    C = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch dual(B, C)
end

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

    for test_grade in 5:8
        B = Blade(randn(test_dim, test_grade))
        C = Scalar(test_value)
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

@testset "dual(B::Blade, C::One)" begin
    test_dim = 10
    for test_grade in 5:8
        B = Blade(randn(test_dim, test_grade))
        C = One()
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

@testset "dual(B::Blade, C::Real)" begin
    test_dim = 9
    test_value = get_random_value()

    for test_grade in 5:8
        B = Blade(randn(test_dim, test_grade))
        C = test_value
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

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

    test_dim = 20
    grade_C = 1

    C = rand(test_dim, grade_C)[:, grade_C]

    # --- Tests

    # ------ dim(B) == dim(C), C is proportional to B

    B_vectors = C * rand()
    B = Blade(B_vectors)

    relative_orientation =
        sign(LinearAlgebra.det(transpose(basis(B)) * reshape(C, length(C), 1)))
    expected_result = Scalar(relative_orientation * volume(B))

    B_dual_C = dual(B, C)
    @test B_dual_C ≈ expected_result

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

    B = Blade(rand(test_dim + 1, grade_C))
    @test_throws DimensionMismatch dual(B, C)

    # ------ B not contained in C

    # grade(B) > grade(C)
    B = Blade(rand(test_dim, grade_C + 1))
    @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)

    # grade(B) == grade(C), B not contained in C
    B_basis = copy(C)
    B_basis[1] += 1
    B = Blade(B_basis)
    @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
end

# --- B::Pseudoscalar

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

    test_dim = 10

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

    # --- Tests

    # dim(B) == dim(C)
    for test_grade in 2:5
        B = Pseudoscalar(test_dim, test_value)
        C = Blade(rand(test_dim, test_grade))
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end

    # dim(B) != dim(C)
    B = Pseudoscalar(test_dim + 1, test_value)
    C = Blade(rand(test_dim, 3))
    @test_throws DimensionMismatch dual(B, C)
end

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

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

    # --- Tests

    # dim(B) == dim(C), test_value_1 != 1
    test_value_1 = get_random_value(2)  # add 2 to keep value away from 0 and 1

    for test_dim in 10:13
        B = Pseudoscalar(test_dim, test_value_1)
        C = Pseudoscalar(test_dim, test_value_2)
        B_dual_C = dual(B, C)
        @test B_dual_C isa Scalar
        @test B_dual_C == test_value_1

        # Check dual(dual(B, C), C) = (-1)^(grade(C) * (grade(C) - 1) / 2) B
        if mod(grade(C), 4) < 2
            @test dual(B_dual_C, C) == B
        else
            @test dual(B_dual_C, C) == -B
        end
    end

    # dim(B) == dim(C), test_value_1 == 1
    test_value_1 = 1

    for test_dim in 10:13
        B = Pseudoscalar(test_dim, test_value_1)
        C = Pseudoscalar(test_dim, test_value_2)
        B_dual_C = dual(B, C)
        @test B_dual_C isa One

        # Check dual(dual(B, C), C) = (-1)^(grade(C) * (grade(C) - 1) / 2) B
        if mod(grade(C), 4) < 2
            @test dual(B_dual_C, C) == B
        else
            @test dual(B_dual_C, C) == -B
        end
    end

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

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

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

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

    # --- Tests

    for test_dim in 10:13
        B = Pseudoscalar(test_dim, test_value_1)
        C = Scalar(test_value_2)
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

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

    for test_dim in 10:13
        B = Pseudoscalar(test_dim, test_value)
        C = One()
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

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

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

    # --- Tests

    for test_dim in 10:13
        B = Pseudoscalar(test_dim, test_value_1)
        C = test_value_2
        @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)
    end
end

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

    test_dim = 10
    test_value = get_random_value()

    # --- Tests

    # dim(B) == dim(C) > 1
    B = Pseudoscalar(test_dim, test_value)
    C = rand(test_dim, 1)[:, 1]
    @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)

    # dim(B) == dim(C) == 1
    B = Pseudoscalar(1, test_value)
    C = [rand()]

    B_dual_C = dual(B, C)
    @test B_dual_C isa Scalar
    @test dual(B, C) == Scalar(test_value)

    # dim(B) != dim(C)
    B = Pseudoscalar(test_dim + 1, test_value)
    C = rand(test_dim, 1)[:, 1]
    @test_throws DimensionMismatch dual(B, C)
end

# --- B::Scalar

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

    test_dim = 15

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

    # --- Tests

    for test_grade in 5:8
        B = Scalar(test_value)
        C = Blade(randn(test_dim, test_grade))

        B_dual_C = dual(B, C)
        expected_result = mod(test_grade, 4) < 2 ?
            Blade(C, volume=test_value) :
            Blade(C, volume=-test_value)

        @test B_dual_C isa Blade
        @test B_dual_C == expected_result
    end
end

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

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

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

    # --- Tests

    for test_dim in 5:8
        B = Scalar(test_value_1)
        C = Pseudoscalar(test_dim, test_value_2)

        expected_result = mod(test_dim, 4) < 2 ?
            Pseudoscalar(test_dim, test_value_1) :
            Pseudoscalar(test_dim, -test_value_1)

        @test dual(B, C) == expected_result
    end
end

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

    B = Scalar(test_value_1)
    C = Scalar(test_value_2)

    B_dual_C = dual(B, C)
    @test B_dual_C isa Scalar
    @test B_dual_C == test_value_1
end

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

    C = One()

    B_dual_C = dual(B, C)
    @test B_dual_C === B
end

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

    B = Scalar(test_value_1)
    C = test_value_2

    B_dual_C = dual(B, C)
    @test B_dual_C isa Scalar
    @test B_dual_C == test_value_1
end

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

    test_dim = 15
    test_value = get_random_value()

    # --- Tests

    B = Scalar(test_value)
    C = randn(test_dim, 1)[:, 1]

    B_dual_C = dual(B, C)
    expected_result = Blade(C, volume=test_value)

    @test B_dual_C isa Blade
    @test B_dual_C == expected_result
end

# --- B::One

@testset "dual(B::One, C::Blade)" begin
    # Preparations
    test_dim = 12

    for test_grade in 5:8
        B = One()
        C = Blade(randn(test_dim, test_grade))

        B_dual_C = dual(B, C)
        expected_result = mod(test_grade, 4) < 2 ?
            Blade(C, volume=1, copy_basis=false) :
            Blade(C, volume=-1, copy_basis=false)

        @test B_dual_C isa Blade
        @test B_dual_C == expected_result
    end
end

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

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

    # --- Tests

    for test_dim in 15:18
        B = One()
        C = Pseudoscalar(test_dim, test_value)

        expected_result = mod(test_dim, 4) < 2 ?
            Pseudoscalar(test_dim, 1) :
            Pseudoscalar(test_dim, -1)

        @test dual(B, C) == expected_result
    end
end

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

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

    B_dual_C = dual(B, C)
    @test B_dual_C === B
end

@testset "dual(B::One, C::One)" begin
    B = One()
    C = One()
    B_dual_C = dual(B, C)
    @test B_dual_C === B
end

@testset "dual(B::One, C::Vector)" begin
    # Preparations
    test_dim = 12

    B = One()
    C = randn(test_dim, 1)[:, 1]

    B_dual_C = dual(B, C)
    expected_result = Blade(C, volume=1)

    @test B_dual_C isa Blade
    @test B_dual_C == expected_result
end

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

    B_dual_C = dual(B, C)
    @test B_dual_C === B
end

# --- B::Zero or C::Zero

@testset "dual(B::Zero, C)" begin
    # Note: the case C::Zero is excluded from this set of tests

    # --- Preparations

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

    B = Zero()
    expected_error = "The dual of Zero is not well-defined"

    # --- Tests

    # C::Real
    C = test_value
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::Vector
    C = rand(4, 1)[:, 1]
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::One
    C = One()
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::Scalar
    C = Scalar(test_value)
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::Blade
    C = Blade(randn(4, 3))
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::Pseudoscalar
    C = Pseudoscalar(5, test_value)
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::Multivector
    # C = Multivector(5, test_value)
    @test_skip 1
end

@testset "dual(B, C::Zero)" begin
    # --- Preparations

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

    C = Zero()
    expected_error = "The dual of anything relative to Zero is not well-defined"

    # --- Tests

    # B::Real
    B = test_value
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::Vector is not zero
    B = rand(4, 1)[:, 1]
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::Vector is zero
    B = [0; 0; 0]
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::Zero
    B = Zero()
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::One
    B = One()
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::Scalar
    B = Scalar(test_value)
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::Blade
    B = Blade(randn(4, 3))
    @test_throws ErrorException(expected_error) dual(B, C)

    # B::Pseudoscalar
    B = Pseudoscalar(5, test_value)
    @test_throws ErrorException(expected_error) dual(B, C)

    # C::Multivector
    # C = Multivector(5, test_value)
    @test_skip 1
end

# --- B::Real

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

    test_dim = 15
    test_value = get_random_value()

    # --- Tests

    for test_grade in 5:8
        B = test_value
        C = Blade(randn(test_dim, test_grade))

        B_dual_C = dual(B, C)
        expected_result = mod(test_grade, 4) < 2 ?
            Blade(C, volume=test_value) :
            Blade(C, volume=-test_value)

        @test B_dual_C isa Blade
        @test B_dual_C == expected_result
    end
end

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

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

    # --- Tests

    for test_dim in 5:8
        B = test_value_1
        C = Pseudoscalar(test_dim, test_value_2)

        expected_result = mod(test_dim, 4) < 2 ?
            Pseudoscalar(test_dim, test_value_1) :
            Pseudoscalar(test_dim, -test_value_1)

        @test dual(B, C) == expected_result
    end
end

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

    B = test_value_1
    C = Scalar(test_value_2)

    B_dual_C = dual(B, C)
    @test B_dual_C isa Scalar
    @test B_dual_C == test_value_1
end

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

    B = test_value
    C = One()

    B_dual_C = dual(B, C)
    @test B_dual_C isa Scalar
    @test B_dual_C == test_value
end

# --- B::Vector

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

    test_dim = 20

    C = Blade(rand(test_dim, 7))

    # --- Tests

    # ------ dim(B) == dim(C), grade(C) > 1, B is contained in C

    grade_B = 1
    B = basis(C)[:, grade_B] * rand()

    B_dual_C = dual(B, C)

    # --- Check dim, grade, and norm

    @test dim(B_dual_C) == test_dim
    @test grade(B_dual_C) == grade(C) - grade_B
    @test norm(B_dual_C) == norm(B)

    # --- Check that B and dual(B) are orthogonal complements

    @test LinearAlgebra.norm(transpose(basis(B_dual_C)) * reshape(B, length(B), 1)) <
        10 * eps(Float64)

    # --- Check sign(dual(B, C))

    # Compute sign of I_C formed from basis(B) and basis(dual(B, C))
    expected_sign = sign(LinearAlgebra.det(
        transpose(hcat(B, basis(B_dual_C))) * basis(C)))

    # Account for sign of I_C^{-1} relative to I_C
    if mod(grade(C), 4) >= 2
        expected_sign = -expected_sign
    end

    @test sign(B_dual_C) == expected_sign

    # Check dual(dual(B, C), C) = (-1)^(grade(C) * (grade(C) - 1) / 2) B
    if mod(grade(C), 4) < 2
        @test dual(B_dual_C, C) ≈ B
    else
        @test dual(B_dual_C, C) ≈ -B
    end

    # ------ dim(B) == dim(C), grade(C) == 1, B is contained in C

    C_grade_1 = Blade(rand(test_dim, grade_B))
    coefficients = rand()
    B = basis(C_grade_1)[:, grade_B] * coefficients

    relative_orientation =
        sign(LinearAlgebra.det(reshape(B, 1, length(B)) * basis(C_grade_1)))
    expected_result = Scalar(relative_orientation * LinearAlgebra.norm(B))

    @test dual(B, C_grade_1) ≈ expected_result

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

    B = rand(test_dim + 1, grade_B)[:, grade_B]
    @test_throws DimensionMismatch dual(B, C)

    # ------ B not contained in C

    # grade(B) <= grade(C), B not contained in C
    B = rand(test_dim, grade_B)[:, grade_B]
    while LinearAlgebra.norm(reject(reshape(B, length(B), 1), C)) < sqrt(eps(Float64))
        B = rand(test_dim, grade_B)[:, grade_B]
    end
    @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)

    # ------ B is a vector of 0s

    B = fill(0, test_dim)
    @test_throws ErrorException("The dual of Zero is not well-defined") dual(B, C)
end

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

    test_dim = 10
    test_value = get_random_value()

    # --- Tests

    # dim(B) == dim(C)
    B = rand(test_dim, 1)[:, 1]
    C = Pseudoscalar(test_dim, test_value)

    B_dual_C = dual(B, C)
    @test B_dual_C isa Blade
    @test B_dual_C == dual(Blade(B))

    # Check dual(dual(B, C), C) = (-1)^(grade(C) * (grade(C) - 1) / 2) B
    if mod(grade(C), 4) < 2
        @test dual(B_dual_C, C) ≈ B
    else
        @test dual(B_dual_C, C) ≈ -B
    end

    # dim(B) != dim(C)
    B = rand(test_dim, 1)[:, 1]
    C = Pseudoscalar(test_dim + 1, test_value)
    @test_throws DimensionMismatch dual(B, C)
end

@testset "dual(B::Vector, C::Scalar)" begin
    C = Scalar(get_random_value())

    # B is not zero
    B = randn(9, 1)[:, 1]
    @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)

    # B is zero
    B = [0; 0; 0]
    @test_throws ErrorException("The dual of Zero is not well-defined") dual(B, C)
end

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

    # B is not zero
    B = randn(10, 1)[:, 1]
    @test_throws ArgumentError("`B` not contained in `C`") dual(B, C)

    # B is zero
    B = [0; 0; 0]
    @test_throws ErrorException("The dual of Zero is not well-defined") dual(B, C)
end

ext: jl


text:
# Example how to use the residual block
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using Test, LinearAlgebra, InvertibleNetworks, NNlib

# Input
nx1 = 32
nx2 = 32
nx_in = 8
n_hidden = 16   # same for x and y
batchsize = 2

ny1 = 64
ny2 = 44
ny_in = 4

# Input image
X = glorot_uniform(nx1, nx2, nx_in, batchsize)
D = glorot_uniform(ny1, ny2, ny_in, batchsize)

# Residual blocks
RB = ConditionalResidualBlock(nx1, nx2, nx_in, ny1, ny2, ny_in, n_hidden, batchsize)

# Observed data
Y, D_ = RB.forward(X, D)

# Set data residual to zero
ΔY = Y.*0f0; ΔD = D.*0f0

# Backward pass
ΔX, ΔD = RB.backward(ΔY, ΔD, X, D)


text:
module NoveltyColors

export ColorDict
import JSON: parsefile

directory = dirname(@__FILE__) * "/../data/"

ColorDict = Dict{AbstractString,Any}()
for file in readdir(directory)
    if file[end-4:end] == ".json"
        merge!(ColorDict, Dict{AbstractString,Any}(file[1:end-5] => (parsefile(directory * file))))
    end
end

end # module


text:
using Pkg
Pkg.activate(pwd())

# # Scope of variables
# ## Local scope

function f()
    z = 42
    return
end

#+

f()
z

#+

function f()
    global z = 42
    return
end

#+

f()
z

#+

function f()
    z = 42
    return z
end

#+

z = f()
z

# ## Global scope

module A
    a = 1 # a global in A's scope
    b = 2 # b global in A's scope
end

a # errors as Main's global scope is separate from A's

#+

using .A: b # make variable b from module A available

A.a
b

#+

b = 4
c = 10

#+

foo(x) = x + c

#+

foo(1)

#+

x = rand(10);
y = rand(10);

f_global() = x .+ y
f_local(x, y) = x .+ y

#+

hcat(f_global(), f_local(x, y))

#+

@time f_global();
@time f_local(x, y);

#+

a, b = 1:10, 11:20;

@time f_local(a, b);
@time f_local(a, b);


text:
using Unitful, ImagineInterface, Imagine, Statistics
using Test

#NOTE: the piezo should be turned on, connected, and in closed loop mode in order for this test to pass
srate = 100000*Unitful.s^-1
pmin = 0.0*Unitful.μm
pmax = 200.0*Unitful.μm
stack_time = 3.0*Unitful.s
reset_time = 2.0*Unitful.s
exp_time = 0.011*Unitful.s
flash_frac = 0.1 #fraction of time to keep laser on during exposure
z_spacing = 3.1*Unitful.μm
z_pad = 5.0*Unitful.μm
nstacks = 10
d = gen_unidirectional_stack(pmin, pmax, z_spacing, stack_time, stack_time, exp_time, srate, flash_frac; z_pad = z_pad)

ocpi2 = rigtemplate("ocpi-2"; sample_rate = srate)
pos = getpositioners(ocpi2)[1]
pos_mon = getname(ocpi2, monitor_name(pos))
append!(pos, "uni_stack_pos", d["positioner"])
replicate!(pos, nstacks-1)
#write_commands("test.json", ocpi2, nstacks, nframes, exp_time; isbidi = false)

aos = [pos;]
ais = [pos_mon;]

sigs = run_imagine("", vcat(aos,ais); ai_trig_dest = "PFI2", ao_trig_dest = "PFI1", trigger_source = "Port2/Line0")

c = Statistics.cor(ustrip.(get_samples(sigs[1]; sampmap=:volts)), ustrip.(get_samples(pos, sampmap=:volts)))
@test c >= 0.99

#ao2 = getoutputs(getanalog(ocpi2))[2]
#ai = getinputs(getanalog(ocpi2))[3] #This is a generic input.  pos_mon currently has a problem with negative voltages
#ai2 = getinputs(getanalog(ocpi2))[4]

#nsamps_on = 10001
#nsamps_off = 10001
#nsamps_back_on = 10001
#nsamps = nsamps_on+nsamps_off+nsamps_back_on
#append!(pos, "on", fill(3.0*Unitful.V, nsamps_on))
#append!(pos, "off", fill(0.0*Unitful.V, nsamps_off))
#append!(pos, "back_on", fill(3.0*Unitful.V, nsamps_back_on))
#append!(ao2, "on2", fill(3.0*Unitful.V, nsamps_on))
#append!(ao2, "off2", fill(0.0*Unitful.V, nsamps_off))
#append!(ao2, "back_on2", fill(3.0*Unitful.V, nsamps_back_on))
#base_name = "" #"test_imaginectrl5"

#using Plots
#plot(ustrip(get_samples(sigs[1])))

#for testing
#tsk = prepare_ao(aos, 1000, "PFI1")
#tsk2 = prepare_ai(ais, 1000, 3003, "disable")
#After much pain, realized that can only synchronize start of AO and AI acq on usb 6xxx devices by using two different pfi channels
#NIDAQ.catch_error(NIDAQ.ExportSignal(tsk.th, NIDAQ.Val_StartTrigger, b"/Dev1/ai/StartTrigger"))
#rr2 = record_signals(base_name, ais, nsamps; trigger_terminal = "PFI1", run_locally = false)
#rr = output_signals(aos; trigger_terminal = "PFI0", run_locally = false)
#sleep(5.0)
#print("delivering pulses\n")
#ttl_pulse() #P0.0 is wired to PFI0
#(proc_id2, ai_sigs) = fetch(rr2)
#proc_id1 = fetch(rr)
#print("Freeing workers...\n")
#free_workers(proc_id1)
#free_workers(proc_id2)
##print("Loading ai file...\n")
##aidata = ImagineInterface.parse_ai(base_name*".ai", map(daq_channel, ais), "dummy-6002", srate)
#using Plots
#plot(ustrip(get_samples(sigs[1])))


text:
function raindrops(number)
    plingplangplong= ""

    if number % 3 == 0
        plingplangplong = "Pling"
    end
    if number % 5 == 0
        plingplangplong *= "Plang"
    end
    if number % 7 == 0
        plingplangplong *= "Plong"
    end

    plingplangplong == "" ? string(number) : plingplangplong

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:
export AbstractSymmetryOperationEmbedding
export AbstractSpaceSymmetryOperationEmbedding

export SitePermutation

export embed
export isidentity

abstract type AbstractSymmetryOperationEmbedding <: AbstractSymmetryOperation end
abstract type AbstractSpaceSymmetryOperationEmbedding <: AbstractSymmetryOperationEmbedding end

"""
    SitePermutation

Represents a permutation of sites as a symmetry operation of a lattice.

# Fields
* `permutation::Permutation`
"""
struct SitePermutation <: AbstractSpaceSymmetryOperationEmbedding
    permutation::Permutation
    SitePermutation(p::Permutation) = new(p)
    SitePermutation(p::AbstractVector{<:Integer}) = new(Permutation(p))
end

Base.:(==)(lhs::SitePermutation, rhs::SitePermutation) = lhs.permutation == rhs.permutation
Base.:(*)(lhs::SitePermutation, rhs::SitePermutation) = SitePermutation(lhs.permutation * rhs.permutation)
Base.:(^)(lhs::SitePermutation, rhs::Integer) = SitePermutation(lhs.permutation^rhs)
Base.inv(sp::SitePermutation) = SitePermutation(inv(sp.permutation))

# import Base.isless
# isless(lhs::SitePermutation, rhs::SitePermutation) = isless(lhs.permutation, rhs.permutation)

Base.hash(arg::SitePermutation, h::UInt) = hash(SitePermutation, hash(arg.permutation, h))

(p::SitePermutation)(i::Integer) = p.permutation(i)

Base.size(p::SitePermutation, args...; kwargs...) = size(p.permutation, args...; kwargs...)

"""
    embed(lattice, translation_operation)

Embed the simplest version of integer translation (no mapping between sites etc.)
"""
function embed(lattice::Lattice, top::TranslationOperation{<:Integer})
    p = zeros(Int, numsite(lattice.supercell))
    for (site_index1, ((site_name1, uc_coord1), _)) in enumerate(lattice.supercell.sites)
        _, uc_coord2 = lattice.hypercube.wrap( top(uc_coord1) )
        site_index1 = getsiteindex(lattice.supercell, (site_name1, uc_coord1))
        site_index2 = getsiteindex(lattice.supercell, (site_name1, uc_coord2))
        p[site_index1] = site_index2
    end
    return SitePermutation(p)
end

"""
    embed(lattice, point_operation)

Embed the simplest version of point operation. (no local unitary operation)
"""
function embed(lattice::Lattice, pop::PointOperation{<:Integer})
    site_map = findsitemap(lattice.unitcell, pop)
    if isnothing(site_map)
        throw(ArgumentError("lattice not compatible with $pop"))
    end
    p = zeros(Int, numsite(lattice.supercell))
    for (i, (j, dR)) in enumerate(site_map)
        namei = getsitename(lattice.unitcell, i)
        namej = getsitename(lattice.unitcell, j)
        for Ri in lattice.bravais_coordinates
            _, Rj = lattice.hypercube.wrap(pop.matrix * Ri + dR)
            i_super = lattice.supercell.siteindices[(namei, Ri)]
            j_super = lattice.supercell.siteindices[(namej, Rj)]
            p[i_super] = j_super
        end
    end
    return SitePermutation(p)
end

"""
    embed(lattice, sop::SpaceOperation{<:Integer, <:Integer})
"""
function embed(lattice::Lattice, sop::SpaceOperation{<:Integer, <:Integer})
    return embed(lattice, PointOperation(sop.matrix)) * embed(lattice, TranslationOperation(sop.displacement))
end

"""
    isidentity(perm::SitePermutation)

Test whether `perm` is an identity.
"""
isidentity(perm::SitePermutation) = isidentity(perm.permutation)


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 07ca21ee-80d8-11eb-26f5-3927b8d5af66
begin
    using PlutoUI
    using Zygote
    using Plots
    using VisualOptim
    plotly()
end

# ╔═╡ 4419ca0c-80ea-11eb-38e0-9dc5d25fa2c2
md"Add Library for Dark mode"

# ╔═╡ dc72ce94-8758-11eb-0846-3b66ba7f962c
html"<button onclick=present()> present </button>"

# ╔═╡ b7d9dd1e-80d7-11eb-0fc1-2d52414a9e70
begin
    using DarkMode: DarkMode
    DarkMode.enable()
    # OR DarkMode.Toolbox(theme="default")
end

# ╔═╡ 0008da96-80eb-11eb-1b44-15ac852eae50
md"
These sliders control the learning rate $\eta$ and the initial value of $x$. Please play around with them to get a sense of how the optimizer works
"

# ╔═╡ b74929b4-80d8-11eb-08d6-f52ce81e4da4
@bind η Slider(0.1:0.05:1; show_value=true, default=0.2)

# ╔═╡ 62e0efe8-8981-11eb-03db-9da240116edb

# ╔═╡ Cell order:
# ╟─4419ca0c-80ea-11eb-38e0-9dc5d25fa2c2
# ╠═dc72ce94-8758-11eb-0846-3b66ba7f962c
# ╠═b7d9dd1e-80d7-11eb-0fc1-2d52414a9e70
# ╠═07ca21ee-80d8-11eb-26f5-3927b8d5af66
# ╟─0008da96-80eb-11eb-1b44-15ac852eae50
# ╠═b74929b4-80d8-11eb-08d6-f52ce81e4da4
# ╠═62e0efe8-8981-11eb-03db-9da240116edb


max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
using CatmullRom, Plots

b = 1/2;
a = 2*b;
fx(t) = a * cospi(t)^2;
fy(t) = (a^2/b) * cospi(t)^3 * sinpi(t);

n_points  =  7;
n_between = 17; # interpolants between each adjacent pair of points

xs = [fx(t) for t=range( 0.0, stop=1.0, length=n_points)];
ys = [fy(t) for t=range( 0.0, stop=1.0, length=n_points)];
points = collect(zip(xs, ys))
close_seq!(points)

crpoints = catmullrom(points, n_between);
cxs, cys = crpoints

plot(xs, ys, linecolor=:lightgreen,  linewidth=3, size=(600,600), legend=nothing, xaxis=nothing, yaxis=nothing)
plot!(cxs, cys, linecolor=:black, linewidth=2, size=(600,600), legend=nothing, xaxis=nothing, yaxis=nothing)

```

Here we use the same parameters on an open section of the Piriform curve.

```julia

n_points  =  9;
n_between = 24; # interpolants between each adjacent pair of points

xs = [fx(t) for t=range( -0.25, stop=0.25, length=n_points)];
ys = [fy(t) for t=range( -0.25, stop=0.25, length=n_points)];
points = collect(zip(xs,ys))

plot(points, linecolor=:blue,  size=(600,600), legend=nothing, xaxis=nothing, yaxis=nothing)
# overlay in fitted curve
cxs, cys = catmullrom(points, n_between);
plot!(cxs, cys, linecolor=:black, size=(600,600), legend=nothing, xaxis=nothing, yaxis=nothing)
```


text:
"""
    replace_all_caps(String)

Replace tokens in ALL CAPS by their lower version and add xxup before.
"""

function replace_all_caps(t)
    t = replace(t, r"([A-Z]+[^a-z\s]*)(?=(\s|$))" => s"xxup \1")
    return replace(t, r"([A-Z]*[^a-z\s]+)(?=(\s|$))" => lowercase)
end

"""
    replace_sentence_case(String)

Replace tokens in Sentence Case by their lower verions and add xxmaj before.
"""
function replace_sentence_case(t)
    t = replace(t, r"(?<!\w)([A-Z][A-Z0-9]*[a-z0-9]+)(?!\w)" => s"xxmaj \1")
    return replace(t, r"(?<!\w)([A-Z][A-Z0-9]*[a-z0-9]+)(?!\w)" => lowercase)
end

convert_lowercase(t) = string("xxbos ", lowercase(t))

## Tests

@testset "Text Transforms" begin
    str1 = "Hello WORLD CAPITAL Sentence Case"

    @test replace_all_caps(str1) == "Hello xxup world xxup capital Sentence Case"
    @test replace_sentence_case(str1) == "xxmaj hello WORLD CAPITAL xxmaj sentence xxmaj case"
    @test convert_lowercase(str1) == "xxbos hello world capital sentence case"
end


text:
"""
    ww = rfweights(orofile, reffile, nf; weightsfn="", varname="", fsmooth=false)

Compute orographic weights from a fine-scale precipitation climatology file.

#Arguments
* `orofile`  : filename of input climatology
* `reffile`  : filename of reference file (for metadata, e.g. the file to downscale)
* `nf`       : refinement factor for spatial downscaling
* `weightsfn`: write weights to file weightsfn
* `varname`  : variable name in climatology
* `fsmooth`  : use smoothing instead of gp conservation

#Returns
* `ww`       : a weight matrix also saved to weightsfn

#Depends

This function uses external system calls using the "cdo" command (https://code.mpimet.mpg.de/projects/cdo/wiki/Cdo) which needs to be available on your system.
"""
function rfweights(orofile, reffile, nf; weightsfn="", varname="", fsmooth=false)

    # Create a reference gridrf.nc file (same grid as rainfarm output files)
    (pr, lon_mat, lat_mat) = read_netcdf2d(reffile, varname)
    # Create fine scale grid
    nss = size(pr)
    if (length(nss)>=3)
        pr = pr[:,:,1]
    end
    println(nss)
    ns = nss[1]
    (lon_f, lat_f) = lon_lat_fine(lon_mat, lat_mat,nf)

    rr = round.(Int, rand(1)*100000)[1]

    println("Output size: ", size(lon_f))
    if (varname=="")
        varname = "pr"
        run(`cdo -s  setname,pr $reffile reffile_rr.nc`)
        reffile = "reffile_rr.nc"
    end
    if (weightsfn=="")
        weightsfn = "weights$rr.nc"
        fsavew = false
    else
        fsavew = true
    end

    # The rest is done in CDO
    println("Computing weights")
    write_netcdf2d("gridrf.nc", reshape(pr, ns, ns, 1), lon_f, lat_f, varname, reffile)
    run(`cdo -s -b F32 timmean $orofile pr_orofile_$rr.nc`)
    run(`cdo -s -f nc copy gridrf.nc gridrf_2_$rr.nc`)
    run(`cdo -s -f nc remapbil,gridrf_2_$rr.nc pr_orofile_$rr.nc pr_remap_rr.nc`)
    if (fsmooth)
        (prr, lon, lat) = read_netcdf2d("pr_remap_rr.nc", "")
        ww = prr./smoothconv(prr, ns)
        write_netcdf2d(weightsfn, ww, lon_f, lat_f, varname, reffile)
        run(`rm -f pr_remap_rr.nc pr_orofile_$rr.nc gridrf.nc gridrf_2_$rr.nc reffile_rr.nc`)
    else
        run(`cdo -s gridboxmean,$nf,$nf pr_remap_rr.nc pr_remap_gbm_$rr.nc`)
        run(`cdo -s remapnn,pr_remap_rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc`)
        run(`cdo -s div pr_remap_rr.nc pr_remap_nn_$rr.nc $weightsfn`)
        run(`rm -f pr_remap_rr.nc pr_remap_gbm_$rr.nc pr_remap_nn_$rr.nc pr_orofile_$rr.nc gridrf.nc gridrf_2_$rr.nc reffile_rr.nc`)
    end
    (ww, lon_f, lat_f) = read_netcdf2d(weightsfn, "")
    if (!fsavew)
        run(`rm -f $weightsfn`)
    end
    return ww
end


text:
result = [2.32551,
-0.155293,
0.469214,
-0.0870457,
0.271364,
-0.222398,
1.01423,
-0.179009,
1.01724,
-0.152972,
0.216039,
-0.251637,
1.05034,
-0.231289,
0.245388,
-0.231289,
1.03371,
-0.258872,
1.03256,
-0.262519,
1.04748,
-0.259143,
1.0535,
-0.266484,
1.04711,
-0.267177,
1.02131,
-0.280381]

include(joinpath(BASE_DIR,"data/data_14bus_pu.jl"))
c_sys14 = System(nodes14, thermal_generators14, loads14, branches14, nothing, 100.0, nothing, nothing, nothing);

include(joinpath(BASE_DIR,"data/data_5bus_pu.jl"))
c_sys5_re = System(nodes5, vcat(thermal_generators5, renewable_generators5), loads5,
                nothing, nothing,  100.0, nothing, nothing, nothing)

@testset begin
    using NLsolve
    pf!, x0 = make_pf(c_sys14);
    res = NLsolve.nlsolve(pf!, x0)
    @test res.zero ≈ result rtol=1e-3

    solve_powerflow!(c_sys14, nlsolve, method = :newton)
    @test_throws PowerSystems.DataFormatError solve_powerflow!(c_sys5_re, nlsolve)
end


text:
abstract type InnerProduct end

dot_arg(f,p::InnerProduct) = evaluate(f,p)
dot_arg(f::AbstractArray,p::InnerProduct) = [f(x) for x in evalpoints(p), f in f]
dot_arg(f::AbstractBasis,p::InnerProduct) = Matrix(f,evalpoints(p))

LinearAlgebra.norm(f, p::InnerProduct) = norm(dot_arg(f,p), p)
LinearAlgebra.dot(f,g,p::InnerProduct) = dot(dot_arg.((f,g),Ref(p))...,p)

"""
    WeightedL2(x,w) <: InnerProduct

# Example
```
p = WeightedL2([0,1],[2,3])
norm(identity,p) -> sqrt(0^2*2 + 1^2*3)
```
"""
struct WeightedL2{X,W} <: InnerProduct
    points::X
    weights::W
end

evalpoints(p::WeightedL2) = p.points
evaluate(f,p::WeightedL2) = f.(evalpoints(p))
dot_arg(f::AbstractArray{<:Number},p::WeightedL2) = f

LinearAlgebra.norm(f::AbstractVecOrMat{<:Number}, p::WeightedL2) = sqrt(dot(f,f,p))

function LinearAlgebra.dot(
    f::AbstractVecOrMat{<:Number},
    g::AbstractVecOrMat{<:Number},
    p::WeightedL2
)
    w = p.weights
    # Workaround for https://github.com/JuliaLang/julia/issues/35424
    length(w) == 0 && return zero(eltype(f))*zero(eltype(w))*zero(eltype(g))
    return f'*Diagonal(w)*g
end

using FastGaussQuadrature

LegendreInnerProduct(n) = WeightedL2(gausslegendre(n)...)
JacobiInnerProduct(n,a,b) = WeightedL2(gaussjacobi(n,a,b)...)


text:
function isvow(a)
  vowels = map(codepoint, ['a','e','i','o','u'])
  map(it -> it in vowels ? Char(it) : it, a)
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
if is_unix()
    cd(joinpath(pwd(), "src"))

    run(`cmake .`)
    run(`make`)
end


text:
module AbstractCircuitComponentRounding

abstract type AbstractPrefNumbSys end

include("UtilityFunctions.jl")

end


text:
module TestMOIHighs

import HiGHS
using Test

const MOI = HiGHS.MOI

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$name" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

const _EXCLUDE = [
    # TODO(odow): investigate test failures
    "test_model_LowerBoundAlreadySet",
    "test_model_UpperBoundAlreadySet",
    "test_model_copy_to_UnsupportedAttribute",
]

function test_runtests()
    model = MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(model, MOI.Test.Config(), exclude = _EXCLUDE)
    return
end

function test_runtests_cache()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            HiGHS.Optimizer(),
        ),
        Float64,
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(),
        exclude = vcat(
            _EXCLUDE,
            "test_conic_NormInfinityCone_VectorAffineFunction",
            "test_conic_NormInfinityCone_VectorOfVariables",
            "test_conic_NormOneCone_VectorAffineFunction",
            "test_conic_NormOneCone_VectorOfVariables",
            # TODO(odow): investigate. These all seem related.
            "test_linear_integration",
            "test_linear_integration_Interval",
            "test_linear_integration_delete_variables",
            "test_modification_delete_variables_in_a_batch",
        ),
    )
    return
end

function test_runtests_simplex()
    model = MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawOptimizerAttribute("solver"), "simplex")
    for presolve in ("on", "off")
        MOI.set(model, MOI.RawOptimizerAttribute("presolve"), presolve)
        MOI.Test.runtests(model, MOI.Test.Config(), exclude = _EXCLUDE)
    end
    return
end

function test_runtests_ipm()
    model = MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawOptimizerAttribute("solver"), "ipm")
    return MOI.Test.runtests(model, MOI.Test.Config(), exclude = _EXCLUDE)
end

function test_runtests_ipm_no_presolve()
    model = MOI.Bridges.full_bridge_optimizer(HiGHS.Optimizer(), Float64)
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawOptimizerAttribute("solver"), "ipm")
    MOI.set(model, MOI.RawOptimizerAttribute("presolve"), "off")
    MOI.Test.runtests(
        model,
        MOI.Test.Config(),
        exclude = vcat(
            _EXCLUDE,
            ["test_conic_linear_INFEASIBLE", "test_conic_linear_INFEASIBLE_2"],
        ),
    )
    return
end

function test_SolverName()
    @test MOI.get(HiGHS.Optimizer(), MOI.SolverName()) == "HiGHS"
    return
end

function test_attributes()
    model = HiGHS.Optimizer()
    @test MOI.get(model, MOI.SolverName()) == "HiGHS"
    @test MOI.get(model, MOI.TimeLimitSec()) > 10000
    MOI.set(model, MOI.TimeLimitSec(), 500)
    @test MOI.get(model, MOI.TimeLimitSec()) == 500.0
    @test MOI.get(model, MOI.RawSolver()) == model
    return
end

function test_MOI_variable_count_and_empty()
    model = HiGHS.Optimizer()
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
    x1 = MOI.add_variable(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 1
    @test MOI.supports_constraint(
        model,
        MOI.VariableIndex,
        MOI.Interval{Float64},
    )
    x2, _ = MOI.add_constrained_variable(model, MOI.Interval(0.0, 1.0))
    @test MOI.get(model, MOI.NumberOfVariables()) == 2
    MOI.empty!(model)
    @test MOI.get(model, MOI.NumberOfVariables()) == 0
end

function test_HiGHS_custom_options()
    model = HiGHS.Optimizer()
    @test MOI.supports(model, MOI.RawOptimizerAttribute("solver"))
    @test MOI.get(model, MOI.RawOptimizerAttribute("solver")) == "choose"
    MOI.set(model, MOI.RawOptimizerAttribute("solver"), "simplex")
    @test MOI.get(model, MOI.RawOptimizerAttribute("solver")) == "simplex"
    @test MOI.get(model, MOI.RawOptimizerAttribute("output_flag")) == true
    MOI.set(model, MOI.RawOptimizerAttribute("output_flag"), false)
    @test MOI.get(model, MOI.RawOptimizerAttribute("output_flag")) == false
    @test MOI.get(model, MOI.RawOptimizerAttribute("time_limit")) > 1000
    MOI.set(model, MOI.RawOptimizerAttribute("time_limit"), 1000.0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("time_limit")) == 1000.0
    # unsupported test
    @test MOI.supports(model, MOI.RawOptimizerAttribute("wrong_param")) == false
    @test_throws(
        MOI.UnsupportedAttribute,
        MOI.get(model, MOI.RawOptimizerAttribute("wrong_param")),
    )
    for v in [false, 1, 1.0, "A"]
        @test_throws(
            MOI.UnsupportedAttribute,
            MOI.set(model, MOI.RawOptimizerAttribute("wrong_param"), v)
        )
    end
    return
end

function test_show()
    model = HiGHS.Optimizer()
    @test sprint(show, model) == "A HiGHS model with 0 columns and 0 rows."
    return
end

function test_options()
    model = HiGHS.Optimizer()
    options = [
        "write_solution_to_file",   # Bool
        "simplex_strategy",         # Cint
        "time_limit",               # Cdouble
        "presolve",                 # String
    ]
    for key in options
        v = MOI.get(model, MOI.RawOptimizerAttribute(key))
        MOI.set(model, MOI.RawOptimizerAttribute(key), v)
        v2 = MOI.get(model, MOI.RawOptimizerAttribute(key))
        @test v == v2
    end
    return
end

function test_option_invalid()
    model = HiGHS.Optimizer()
    @test_throws(
        ErrorException(
            "Encountered an error in HiGHS: Check the log for details.",
        ),
        MOI.set(model, MOI.RawOptimizerAttribute("time_limit"), -1.0),
    )
    return
end

function test_option_unknown_option()
    model = HiGHS.Optimizer()
    err = ErrorException(
        "Encountered an error in HiGHS: Check the log for details.",
    )
    @test_throws(
        err,
        MOI.set(model, MOI.RawOptimizerAttribute("write_solution_to_file"), 1),
    )
    @test_throws(
        err,
        MOI.set(model, MOI.RawOptimizerAttribute("simplex_strategy"), "on"),
    )
    @test_throws err MOI.set(model, MOI.RawOptimizerAttribute("time_limit"), 1)
    @test_throws err MOI.set(model, MOI.RawOptimizerAttribute("presolve"), 1)
    return
end

function test_copy_to()
    src = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(
        src,
        """
variables: w, x, y, z
minobjective: w + x + y + z
c1: w <= 1.0
c2: w >= 0.5
c3: x in MathOptInterface.Interval(1.0, 2.0)
c4: y == 2.0
c5: y + z >= 4.5
c6: x + y <= 3.0
c7: w + x == 1.5
c8: w + x in MathOptInterface.Interval(1.0, 2.0)
""",
    )
    dest = HiGHS.Optimizer()
    MOI.copy_to(dest, src)
    @test MOI.get(dest, MOI.NumberOfVariables()) == 4
    list = MOI.get(dest, MOI.ListOfConstraintTypesPresent())
    @test length(list) == 8
    for S in (
        MOI.GreaterThan{Float64},
        MOI.LessThan{Float64},
        MOI.EqualTo{Float64},
        MOI.Interval{Float64},
    )
        @test (MOI.VariableIndex, S) in list
        @test (MOI.ScalarAffineFunction{Float64}, S) in list
    end
    MOI.optimize!(dest)
    @test MOI.get(dest, MOI.ObjectiveValue()) == 0.5 + 1.0 + 2.0 + 2.5
    return
end

end

TestMOIHighs.runtests()


text:
function hyperbolic_cross_derivative(weights::Array{T,1},node::Array{R,1},multi_index::Array{S,2},pos::S) where {T<:AbstractFloat,R<:Number,S<:Integer}

  n = 2*maximum(multi_index,dims=1).+1
  d = length(n) # d is the number of dimensions

  base_polynomials                   = Array{Array{R,2},1}(undef,d)
  base_polynomial_derivatives        = Array{Array{R,2},1}(undef,d)
  unique_base_polynomials            = Array{Array{R,1},2}(undef,size(multi_index))
  unique_base_polynomial_derivatives = Array{Array{R,1},2}(undef,size(multi_index))

  # Construct the base polynomials

  for i = 1:d
    base_polynomials[i] = chebyshev_polynomial(n[i]-1,node[i])
    base_polynomial_derivatives[i] = derivative_of_chebyshev_polynomial(n[i],node[i])
  end

  # Compute the unique polynomial terms from the base polynomials

  for i in axes(multi_index,1)
   for j = 1:d
      if multi_index[i,j] == 0
        unique_base_polynomials[i,j]            = [base_polynomials[j][1]]
        unique_base_polynomial_derivatives[i,j] = [base_polynomial_derivatives[j][1]]
      else
        unique_base_polynomials[i,j]            = [base_polynomials[j][2*multi_index[i,j]], base_polynomials[j][2*multi_index[i,j]+1]]
        unique_base_polynomial_derivatives[i,j] = [base_polynomial_derivatives[j][2*multi_index[i,j]], base_polynomial_derivatives[j][2*multi_index[i,j]+1]]
      end
    end
  end

  # Construct the first row of the interplation matrix

  polynomials = Array{R,1}(undef,length(weights))

  l = 1
  @inbounds for j in axes(multi_index,1)
    if pos == 1
      new_polynomials = unique_base_polynomial_derivatives[j,1]
    else
      new_polynomials = unique_base_polynomials[j,1]
    end
    for i = 2:size(multi_index,2)
      if pos == i
        new_polynomials = kron(new_polynomials,unique_base_polynomial_derivatives[j,i])
      else
        new_polynomials = kron(new_polynomials,unique_base_polynomials[j,i])
      end
    end
    m = length(new_polynomials)
    polynomials[l:l+m-1] = new_polynomials
    l += m
  end

  evaluated_derivative = zero(T)

  for i in eachindex(polynomials)
    evaluated_derivative += polynomials[i]*weights[i]
  end

  return evaluated_derivative

end

function hyperbolic_cross_derivative(weights::Array{T,1},node::Array{R,1},multi_index::Array{S,2},domain::Union{Array{T,1},Array{T,2}},pos::S) where {T<:AbstractFloat,R<:Number,S<:Integer}

  node = copy(node)

  if size(domain,2) != length(node)
    error("domain is inconsistent with the number of dimensions")
  end

  d = length(node)
  for i = 1:d
    node[i] = normalize_node(node[i],domain[:,i])
  end

  evaluated_derivative = hyperbolic_cross_derivative(weights,node,multi_index,pos)

  return evaluated_derivative

end


text:
#!/usr/bin/env julia
using Pkg;

# FIXME: Specify versions here :'(
Pkg.add.([
    "PyPlot",
    "PyCall",
    "Roots",
    "Polynomials",
    "DifferentialEquations",
    "Triangle",
    "Optim",
    "Images",
    "Plots",
    "Mimi"
]);

# Add support for various MIMI models
Pkg.Registry.add(RegistrySpec(url = "https://github.com/mimiframework/MimiRegistry.git"))

Pkg.add.([
    "MimiFUND"
])
# Precompiling installed packages
Pkg.precompile()

max_stars_repo_name: ryanlovett/datahub

max_issues_repo_name: ryanlovett/datahub

max_forks_repo_name: ryanlovett/datahub


text:
### fixes

using CUDArt
using CUSPARSE
using CUSPARSE: CudaSparseMatrix

Base.convert(::Type{CudaSparseMatrixCSC}, x::SparseMatrixCSC)=CudaSparseMatrixCSC(x)
Base.convert{T<:Array}(::Type{T},a::CudaSparseMatrix)=full(to_host(a))
Base.isempty(a::CudaSparseMatrix) = (length(a) == 0)
Base.issparse(a::CudaSparseMatrix) = true
# Base.ndims(::CudaSparseMatrix) = 2
Base.nnz(x::CudaSparseMatrix)=x.nnz
Base.scale!(c,a::CudaSparseMatrix) = (scale!(c,a.nzVal); a)
Base.stride(g::CudaSparseMatrix,i)=(i==1 ? 1 : i==2 ? g.dims[1] : length(g))
Base.strides(g::CudaSparseMatrix)=(1,g.dims[1])
# Base.summary(a::CudaSparseMatrix) = string(Base.dims2string(size(a)), " ", typeof(a))
Base.vecnorm(a::CudaSparseMatrix) = vecnorm(a.nzVal)
Base.LinAlg.BLAS.nrm2(a::CudaSparseMatrix) = Base.LinAlg.BLAS.nrm2(a.nzVal)

function copysync!{T}(a::CudaSparseMatrixCSC{T}, b::SparseMatrixCSC{T})
    a.dims = (b.m,b.n)
    a.nnz = convert(Cint, length(b.nzval))
    resizecopy!(a.colPtr, convert(Vector{Cint},b.colptr))
    resizecopy!(a.rowVal, convert(Vector{Cint},b.rowval))
    resizecopy!(a.nzVal, b.nzval)
    gpusync(); return a
end

function copysync!{T}(a::CudaSparseMatrixCSC{T}, b::CudaSparseMatrixCSC{T})
    a.dims = b.dims
    a.nnz = b.nnz
    resizecopy!(a.colPtr, convert(Vector{Cint},b.colPtr))
    resizecopy!(a.rowVal, convert(Vector{Cint},b.rowVal))
    resizecopy!(a.nzVal, b.nzVal)
    gpusync(); return a
end

function copysync!{T}(a::CudaSparseMatrixCSR{T}, b::CudaSparseMatrixCSR{T})
    a.dims = b.dims
    a.nnz = b.nnz
    resizecopy!(a.rowPtr, convert(Vector{Cint},b.rowPtr))
    resizecopy!(a.colVal, convert(Vector{Cint},b.colVal))
    resizecopy!(a.nzVal, b.nzVal)
    gpusync(); return a
end

function resizecopy!{T}(a::CudaVector{T}, b::Vector{T})
    resize!(a, length(b))       # TODO: is this efficient?
    copysync!(a, b)
end

function resizecopy!{T}(a::CudaVector{T}, b::CudaVector{T})
    resize!(a, length(b))       # TODO: is this efficient?
    copysync!(a, b)
end

CUDArt.free(x::CudaSparseMatrixCSR)=(free(x.rowPtr);free(x.colVal);free(x.nzVal))

function fillsync!(x::CudaSparseMatrixCSR,n)
    n == 0 || error("Only 0 fill for sparse")
    fillsync!(x.rowPtr,1)
    resize!(x.colVal,0)
    resize!(x.nzVal,0)
    x.nnz = 0
    return x
end

function fillsync!(x::CudaSparseMatrixCSC,n)
    n == 0 || error("Only 0 fill for sparse")
    fillsync!(x.colPtr,1)
    resize!(x.rowVal,0)
    resize!(x.nzVal,0)
    x.nnz = 0
    return x
end

# A_mul_Bt!(csr,dns,csc) is more efficient if we do not insist on
# unique sorted csr entries.  To signal unsorted csr, we introduce a
# new type so we don't accidentally pass it to an unprepared function.
# The following adapted from CUSPARSE.jl:

# For dw = r' * du we need the CSC versions below

if !isdefined(:CudaSparseMatrixCSRU)
    type CudaSparseMatrixCSRU{T}
        rowPtr::CudaVector{Cint}
        colVal::CudaVector{Cint}
        nzVal::CudaVector{T}
        dims::NTuple{2,Int}
        nnz::Cint
        dev::Int
    end
end

if !isdefined(:CudaSparseMatrixCSCU)
    type CudaSparseMatrixCSCU{T}
        colPtr::CudaVector{Cint}
        rowVal::CudaVector{Cint}
        nzVal::CudaVector{T}
        dims::NTuple{2,Int}
        nnz::Cint
        dev::Int
    end
end

CudaSparseMatrixCSRU(T::Type, m::Integer, n::Integer)=CudaSparseMatrixCSRU{T}(fillsync!(CudaArray(Cint,m+1),1), CudaArray(Cint,0), CudaArray(T,0), (convert(Int,m),convert(Int,n)), convert(Cint,0), convert(Int,device()))
CudaSparseMatrixCSRU(T::Type, rowPtr::CudaArray, colVal::CudaArray, nzVal::CudaArray, dims::NTuple{2,Int}) = CudaSparseMatrixCSRU{T}(rowPtr, colVal, nzVal, dims, convert(Cint,length(nzVal)), device())
CudaSparseMatrixCSRU(T::Type, rowPtr::CudaArray, colVal::CudaArray, nzVal::CudaArray, nnz, dims::NTuple{2,Int}) = CudaSparseMatrixCSRU{T}(rowPtr, colVal, nzVal, dims, nnz, device())

CudaSparseMatrixCSCU(T::Type, m::Integer, n::Integer)=CudaSparseMatrixCSCU{T}(fillsync!(CudaArray(Cint,n+1),1), CudaArray(Cint,0), CudaArray(T,0), (convert(Int,m),convert(Int,n)), convert(Cint,0), convert(Int,device()))
CudaSparseMatrixCSCU(T::Type, colPtr::CudaArray, rowVal::CudaArray, nzVal::CudaArray, dims::NTuple{2,Int}) = CudaSparseMatrixCSCU{T}(colPtr, rowVal, nzVal, dims, convert(Cint,length(nzVal)), device())
CudaSparseMatrixCSCU(T::Type, colPtr::CudaArray, rowVal::CudaArray, nzVal::CudaArray, nnz, dims::NTuple{2,Int}) = CudaSparseMatrixCSCU{T}(colPtr, rowVal, nzVal, dims, nnz, device())
CudaSparseMatrixCSCU(T::Type, colPtr::Vector, rowVal::Vector, nzVal::Vector, dims::NTuple{2,Int}) = CudaSparseMatrixCSCU{T}(CudaArray(convert(Vector{Cint},colPtr)), CudaArray(convert(Vector{Cint},rowVal)), CudaArray(nzVal), dims, convert(Cint,length(nzVal)), device())
CudaSparseMatrixCSCU(Mat::SparseMatrixCSC) = CudaSparseMatrixCSCU(eltype(Mat), Mat.colptr, Mat.rowval, Mat.nzval, size(Mat))

Base.eltype{T}(x::CudaSparseMatrixCSRU{T})=T
Base.size(x::CudaSparseMatrixCSRU)=x.dims
Base.issparse(x::CudaSparseMatrixCSRU)=true
Base.scale!(s,x::CudaSparseMatrixCSRU)=scale!(s,x.nzVal)
Base.similar(Mat::CudaSparseMatrixCSRU) = CudaSparseMatrixCSRU(eltype(Mat), copy(Mat.rowPtr), copy(Mat.colVal), similar(Mat.nzVal), Mat.nnz, Mat.dims)
Base.copy(Mat::CudaSparseMatrixCSRU; stream=null_stream) = copysync!(similar(Mat),Mat;stream=null_stream)

Base.eltype{T}(x::CudaSparseMatrixCSCU{T})=T
Base.size(x::CudaSparseMatrixCSCU)=x.dims
Base.issparse(x::CudaSparseMatrixCSCU)=true
Base.scale!(s,x::CudaSparseMatrixCSCU)=scale!(s,x.nzVal)
Base.similar(Mat::CudaSparseMatrixCSCU) = CudaSparseMatrixCSCU(eltype(Mat), copy(Mat.colPtr), copy(Mat.rowVal), similar(Mat.nzVal), Mat.nnz, Mat.dims)
Base.copy(Mat::CudaSparseMatrixCSCU; stream=null_stream) = copysync!(similar(Mat),Mat;stream=null_stream)

function copysync!(dst::CudaSparseMatrixCSRU, src::CudaSparseMatrixCSRU; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Matrix size"))
    end
    resizecopy!( dst.rowPtr, src.rowPtr )
    resizecopy!( dst.colVal, src.colVal )
    resizecopy!( dst.nzVal, src.nzVal )
    dst.nnz = src.nnz
    gpusync(); return dst
end

function copysync!(dst::CudaSparseMatrixCSCU, src::CudaSparseMatrixCSCU; stream=null_stream)
    if dst.dims != src.dims
        throw(ArgumentError("Inconsistent Sparse Matrix size"))
    end
    resizecopy!( dst.colPtr, src.colPtr )
    resizecopy!( dst.rowVal, src.rowVal )
    resizecopy!( dst.nzVal, src.nzVal )
    dst.nnz = src.nnz
    gpusync(); return dst
end

function CUDArt.to_host{T}(Mat::CudaSparseMatrixCSRU{T})
    rowPtr = to_host(Mat.rowPtr)
    colVal = to_host(Mat.colVal)
    nzVal = to_host(Mat.nzVal)
    #construct Is
    I = similar(colVal)
    counter = 1
    for row = 1 : size(Mat)[1], k = rowPtr[row] : (rowPtr[row+1]-1)
        I[counter] = row
        counter += 1
    end
    return sparse(I,colVal,nzVal,Mat.dims[1],Mat.dims[2])
end

function CUDArt.to_host{T}(Mat::CudaSparseMatrixCSCU{T})
    (m,n) = Mat.dims
    colptr = to_host(Mat.colPtr)
    rowval = to_host(Mat.rowVal)
    nzval = to_host(Mat.nzVal)
    colval = similar(rowval)
    for i=1:n
        colval[colptr[i]:colptr[i+1]-1] = i
    end
    return sparse(rowval, colval, nzval, m, n)
end

deepcopy_internal(x::CudaSparseMatrixCSRU, s::ObjectIdDict)=(haskey(s,x)||(s[x]=copy(x));s[x])
gpucopy_internal(x::CudaSparseMatrixCSRU, s::ObjectIdDict)=deepcopy_internal(x,s)

deepcopy_internal(x::CudaSparseMatrixCSCU, s::ObjectIdDict)=(haskey(s,x)||(s[x]=copy(x));s[x])
gpucopy_internal(x::CudaSparseMatrixCSCU, s::ObjectIdDict)=deepcopy_internal(x,s)

deepcopy_internal(x::CudaSparseMatrix, s::ObjectIdDict)=(haskey(s,x)||(s[x]=copy(x));s[x])
gpucopy_internal(x::CudaSparseMatrix, s::ObjectIdDict)=deepcopy_internal(x,s)
gpucopy_internal{T<:Number}(x::SparseMatrixCSC{T}, s::ObjectIdDict)=(haskey(s,x)||(s[x]=CudaSparseMatrixCSC(x));s[x])

# We need cpu versions of CudaSparseMatrices for cpu/gpucopy to work
if !isdefined(:SparseMatrixCSR0)
    type SparseMatrixCSR0{T}
        rowPtr::Array{Cint,1}
        colVal::Array{Cint,1}
        nzVal::Array{T,1}
        dims::NTuple{2,Int}
        nnz::Cint
        dev::Int
    end
end

if !isdefined(:SparseMatrixCSRU)
    type SparseMatrixCSRU{T}
        rowPtr::Array{Cint,1}
        colVal::Array{Cint,1}
        nzVal::Array{T,1}
        dims::NTuple{2,Int}
        nnz::Cint
        dev::Int
    end
end

if !isdefined(:SparseMatrixCSC0)
    type SparseMatrixCSC0{T}
        colPtr::Array{Cint,1}
        rowVal::Array{Cint,1}
        nzVal::Array{T,1}
        dims::NTuple{2,Int}
        nnz::Cint
        dev::Int
    end
end

if !isdefined(:SparseMatrixCSCU)
    type SparseMatrixCSCU{T}
        colPtr::Array{Cint,1}
        rowVal::Array{Cint,1}
        nzVal::Array{T,1}
        dims::NTuple{2,Int}
        nnz::Cint
        dev::Int
    end
end

# preserve matrix format in cpu
isdefined(:SparseMatrix0) || (typealias SparseMatrix0 Union{SparseMatrixCSC0,SparseMatrixCSR0,SparseMatrixCSCU,SparseMatrixCSRU})
isdefined(:CudaSparseMatrix0) || typealias CudaSparseMatrix0 Union{CudaSparseMatrixCSC,CudaSparseMatrixCSR,CudaSparseMatrixCSCU,CudaSparseMatrixCSRU}

cpucopy_internal(x::CudaSparseMatrix0, s::ObjectIdDict)=(haskey(s,x)||(s[x]=cpucopy_sparse(x));s[x])
gpucopy_internal(x::SparseMatrix0, s::ObjectIdDict)=(haskey(s,x)||(s[x]=gpucopy_sparse(x));s[x])

cpucopy_sparse{T}(x::CudaSparseMatrixCSC{T})=SparseMatrixCSC0{T}(to_host(x.colPtr),to_host(x.rowVal),to_host(x.nzVal),x.dims,x.nnz,x.dev)
cpucopy_sparse{T}(x::CudaSparseMatrixCSCU{T})=SparseMatrixCSCU{T}(to_host(x.colPtr),to_host(x.rowVal),to_host(x.nzVal),x.dims,x.nnz,x.dev)
cpucopy_sparse{T}(x::CudaSparseMatrixCSR{T})=SparseMatrixCSR0{T}(to_host(x.rowPtr),to_host(x.colVal),to_host(x.nzVal),x.dims,x.nnz,x.dev)
cpucopy_sparse{T}(x::CudaSparseMatrixCSRU{T})=SparseMatrixCSRU{T}(to_host(x.rowPtr),to_host(x.colVal),to_host(x.nzVal),x.dims,x.nnz,x.dev)

gpucopy_sparse{T}(x::SparseMatrixCSC0{T})=CudaSparseMatrixCSC{T}(CudaArray(x.colPtr),CudaArray(x.rowVal),CudaArray(x.nzVal),x.dims,x.nnz,x.dev)
gpucopy_sparse{T}(x::SparseMatrixCSCU{T})=CudaSparseMatrixCSCU{T}(CudaArray(x.colPtr),CudaArray(x.rowVal),CudaArray(x.nzVal),x.dims,x.nnz,x.dev)
gpucopy_sparse{T}(x::SparseMatrixCSR0{T})=CudaSparseMatrixCSR{T}(CudaArray(x.rowPtr),CudaArray(x.colVal),CudaArray(x.nzVal),x.dims,x.nnz,x.dev)
gpucopy_sparse{T}(x::SparseMatrixCSRU{T})=CudaSparseMatrixCSRU{T}(CudaArray(x.rowPtr),CudaArray(x.colVal),CudaArray(x.nzVal),x.dims,x.nnz,x.dev)

# equality testing
function Base.isequal(a::Union{CudaSparseMatrix0,SparseMatrix0},b::Union{CudaSparseMatrix0,SparseMatrix0})
    typeof(a)==typeof(b) || return false
    for n in fieldnames(a)
        if isdefined(a,n) && isdefined(b,n)
            isequal(a.(n), b.(n)) || return false
        elseif isdefined(a,n) || isdefined(b,n)
            return false
        end
    end
    return true
end

import Base.LinAlg.BLAS: gemm!
using CUSPARSE: SparseChar, cusparseop, cusparseindex,
    cusparseMatDescr_t, CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT,
    statuscheck, libcusparse, cusparseStatus_t, cusparseHandle_t, cusparseOperation_t, cusparsehandle

for (fname,elty) in ((:cusparseScsrgemm, :Float32),
                     (:cusparseDcsrgemm, :Float64),
                     (:cusparseCcsrgemm, :Complex64),
                     (:cusparseZcsrgemm, :Complex128))
    @eval begin
        function gemm!(transa::SparseChar,
                       transb::SparseChar,
                       A::CudaSparseMatrixCSR{$elty},
                       B::CudaSparseMatrixCSR{$elty},
                       C::CudaSparseMatrixCSR{$elty})
            cutransa = cusparseop(transb)
            cutransb = cusparseop(transa)
            cudesc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cusparseindex('O'))
            m,k  = transa == 'N' ? A.dims : (A.dims[2],A.dims[1])
            kB,n = transb == 'N' ? B.dims : (B.dims[2],B.dims[1])
            if k != kB
                throw(DimensionMismatch("Interior dimension of A, $k, and B, $kB, must match"))
            end
            nnzC = Array(Cint,1)
            resize!(C.rowPtr, m+1)
            statuscheck(ccall((:cusparseXcsrgemmNnz,libcusparse), cusparseStatus_t, # t:186/868
                              (cusparseHandle_t, cusparseOperation_t,
                               cusparseOperation_t, Cint, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Cint, Ptr{Cint},
                               Ptr{Cint}, Ptr{cusparseMatDescr_t}, Cint, Ptr{Cint},
                               Ptr{Cint}, Ptr{cusparseMatDescr_t}, Ptr{Cint},
                               Ptr{Cint}), cusparsehandle[1], cutransa, cutransb,
                              m, n, k, &cudesc, A.nnz, A.rowPtr, A.colVal,
                              &cudesc, B.nnz, B.rowPtr, B.colVal, &cudesc,
                              C.rowPtr, nnzC))
            C.nnz = nnzC[1]
            C.dims = (m,n)
            resize!(C.nzVal, C.nnz)
            resize!(C.colVal, C.nnz)
            statuscheck(ccall(($(string(fname)),libcusparse), cusparseStatus_t, # t:659/868
                              (cusparseHandle_t, cusparseOperation_t,
                               cusparseOperation_t, Cint, Cint, Cint,
                               Ptr{cusparseMatDescr_t}, Cint, Ptr{$elty},
                               Ptr{Cint}, Ptr{Cint}, Ptr{cusparseMatDescr_t},
                               Cint, Ptr{$elty}, Ptr{Cint}, Ptr{Cint},
                               Ptr{cusparseMatDescr_t}, Ptr{$elty}, Ptr{Cint},
                               Ptr{Cint}), cusparsehandle[1], cutransa,
                              cutransb, m, n, k, &cudesc, A.nnz, A.nzVal,
                              A.rowPtr, A.colVal, &cudesc, B.nnz, B.nzVal,
                              B.rowPtr, B.colVal, &cudesc, C.nzVal,
                              C.rowPtr, C.colVal))
            C
        end
    end
end

        # function geam!(alpha::$elty,
        #                A::CudaSparseMatrixCSR{$elty},
        #                beta::$elty,
        #                B::CudaSparseMatrixCSR{$elty},
        #                C::CudaSparseMatrixCSR{$elty},
        #                indexA::SparseChar='O',
        #                indexB::SparseChar='O',
        #                indexC::SparseChar='O')
        #     cuinda = cusparseindex(indexA)
        #     cuindb = cusparseindex(indexB)
        #     cuindc = cusparseindex(indexB)
        #     cudesca = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuinda)
        #     cudescb = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuindb)
        #     cudescc = cusparseMatDescr_t(CUSPARSE_MATRIX_TYPE_GENERAL, CUSPARSE_FILL_MODE_LOWER, CUSPARSE_DIAG_TYPE_NON_UNIT, cuindc)
        #     mA,nA = A.dims
        #     mB,nB = B.dims
        #     if ( (mA != mB) || (nA != nB) )
        #         throw(DimensionMismatch(""))
        #     end
        #     nnzC = Array(Cint,1)
        #     length(C.rowPtr) < mA+1 && (C.rowPtr = CudaArray(zeros(Cint, mA+1)))
        #     statuscheck(ccall((:cusparseXcsrgeamNnz,libcusparse), cusparseStatus_t,
        #                       (cusparseHandle_t, Cint, Cint,
        #                        Ptr{cusparseMatDescr_t}, Cint, Ptr{Cint},
        #                        Ptr{Cint}, Ptr{cusparseMatDescr_t}, Cint, Ptr{Cint},
        #                        Ptr{Cint}, Ptr{cusparseMatDescr_t}, Ptr{Cint},
        #                        Ptr{Cint}), cusparsehandle[1], mA, nA, &cudesca,
        #                       A.nnz, A.rowPtr, A.colVal, &cudescb, B.nnz,
        #                       B.rowPtr, B.colVal, &cudescc, C.rowPtr, nnzC))
        #     C.nnz = nnzC[1]
        #     C.dims = A.dims
        #     length(C.nzVal)  < C.nnz && (C.nzVal = CudaArray($elty,Int(C.nnz)))
        #     length(C.colVal) < C.nnz && (C.colVal = CudaArray(Cint,Int(C.nnz)))
        #     statuscheck(ccall(($(string(fname)),libcusparse), cusparseStatus_t,
        #                       (cusparseHandle_t, Cint, Cint, Ptr{$elty},
        #                        Ptr{cusparseMatDescr_t}, Cint, Ptr{$elty},
        #                        Ptr{Cint}, Ptr{Cint}, Ptr{$elty},
        #                        Ptr{cusparseMatDescr_t}, Cint, Ptr{$elty},
        #                        Ptr{Cint}, Ptr{Cint}, Ptr{cusparseMatDescr_t},
        #                        Ptr{$elty}, Ptr{Cint}, Ptr{Cint}),
        #                       cusparsehandle[1], mA, nA, [alpha], &cudesca,
        #                       A.nnz, A.nzVal, A.rowPtr, A.colVal, [beta],
        #                       &cudescb, B.nnz, B.nzVal, B.rowPtr, B.colVal,
        #                       &cudescc, C.nzVal, C.rowPtr, C.colVal))
        #     C
        # end
# # This is buggy in master:util.jl:75, submitted pull request:
# function to_host2{T}(Mat::CudaSparseMatrixCSR{T})
#     rowPtr = to_host(Mat.rowPtr)
#     colVal = to_host(Mat.colVal)
#     nzVal = to_host(Mat.nzVal)
#     #construct Is
#     I = similar(colVal)
#     counter = 1
#     for row = 1 : size(Mat)[1], k = rowPtr[row] : (rowPtr[row+1]-1)
#         I[counter] = row
#         counter += 1
#     end
#     return sparse(I,colVal,nzVal,Mat.dims[1],Mat.dims[2])
# end


text:
mutable struct Descriptor <: AbstractDescriptor
    name::String
    p::libgb.GrB_Descriptor
    function Descriptor(name, p::libgb.GrB_Descriptor)
        d = new(name, p)
        function f(descriptor)
            libgb.GrB_Descriptor_free(Ref(descriptor.p))
        end
        return finalizer(f, d)
    end
end

const DEFAULT = libgb.GxB_DEFAULT
const REPLACE = libgb.GrB_REPLACE
const COMPLEMENT = libgb.GrB_COMP
const STRUCTURE = libgb.GrB_STRUCTURE
const STRUCT_COMP = libgb.GrB_STRUCT_COMP
const TRANSPOSE = libgb.GrB_TRAN
const GUSTAVSON = libgb.GxB_AxB_GUSTAVSON
const DOT = libgb.GxB_AxB_DOT
const HASH = libgb.GxB_AxB_HASH
const SAXPY = libgb.GxB_AxB_SAXPY

const T1 = Descriptor("T1", libgb.GrB_Descriptor(C_NULL))
const T0 = Descriptor("T0", libgb.GrB_Descriptor(C_NULL))
const T0T1 = Descriptor("T0T1", libgb.GrB_Descriptor(C_NULL))
const C = Descriptor("C", libgb.GrB_Descriptor(C_NULL))
const CT1 = Descriptor("CT1", libgb.GrB_Descriptor(C_NULL))
const CT0 = Descriptor("CT0", libgb.GrB_Descriptor(C_NULL))
const CT0T1 = Descriptor("CT0T1", libgb.GrB_Descriptor(C_NULL))
const S = Descriptor("S", libgb.GrB_Descriptor(C_NULL))
const ST1 = Descriptor("ST1", libgb.GrB_Descriptor(C_NULL))
const ST0 = Descriptor("ST0", libgb.GrB_Descriptor(C_NULL))
const ST0T1 = Descriptor("ST0T1", libgb.GrB_Descriptor(C_NULL))
const SC = Descriptor("SC", libgb.GrB_Descriptor(C_NULL))
const SCT1 = Descriptor("SCT1", libgb.GrB_Descriptor(C_NULL))
const SCT0 = Descriptor("SCT0", libgb.GrB_Descriptor(C_NULL))
const SCT0T1 = Descriptor("SCT0T1", libgb.GrB_Descriptor(C_NULL))
const R = Descriptor("R", libgb.GrB_Descriptor(C_NULL))
const RT1 = Descriptor("RT1", libgb.GrB_Descriptor(C_NULL))
const RT0 = Descriptor("RT0", libgb.GrB_Descriptor(C_NULL))
const RT0T1 = Descriptor("RT0T1", libgb.GrB_Descriptor(C_NULL))
const RC = Descriptor("RC", libgb.GrB_Descriptor(C_NULL))
const RCT1 = Descriptor("RCT1", libgb.GrB_Descriptor(C_NULL))
const RCT0 = Descriptor("RCT0", libgb.GrB_Descriptor(C_NULL))
const RCT0T1 = Descriptor("RCT0T1", libgb.GrB_Descriptor(C_NULL))
const RS = Descriptor("RS", libgb.GrB_Descriptor(C_NULL))
const RST1 = Descriptor("RST1", libgb.GrB_Descriptor(C_NULL))
const RST0 = Descriptor("RST0", libgb.GrB_Descriptor(C_NULL))
const RST0T1 = Descriptor("RST0T1", libgb.GrB_Descriptor(C_NULL))
const RSC = Descriptor("RSC", libgb.GrB_Descriptor(C_NULL))
const RSCT1 = Descriptor("RSCT1", libgb.GrB_Descriptor(C_NULL))
const RSCT0 = Descriptor("RSCT0", libgb.GrB_Descriptor(C_NULL))
const RSCT0T1 = Descriptor("RSCT0T1", libgb.GrB_Descriptor(C_NULL))
const DEFAULTDESC = Descriptor("DEFAULTDESC", libgb.GrB_Descriptor(C_NULL))

function Descriptor()
    return Descriptor("", libgb.GrB_Descriptor_new())
end

Base.unsafe_convert(::Type{libgb.GrB_Descriptor}, d::Descriptor) = d.p

function Base.getproperty(d::Descriptor, s::Symbol)
    if s == :p
        return getfield(d, s)
    elseif s == :output
        f = libgb.GrB_OUTP
    elseif s == :mask
        f = libgb.GrB_MASK
    elseif s == :input1
        f = libgb.GrB_INP0
    elseif s == :input2
        f = libgb.GrB_INP1
    #elseif s == :nthreads
        #f = libgb.GxB_DESCRIPTOR_NTHREADS
    #elseif s == :chunk
        #f = libgb.GxB_DESCRIPTOR_CHUNK
    #elseif s == :axb_method
        #f = libgb.GxB_AxB_METHOD
    #elseif s == :sort
        #f = libgb.GxB_SORT
    else
        return getfield(d, s)
    end
    return libgb.GxB_Descriptor_get(d, f)
end

function Base.setproperty!(d::Descriptor, s::Symbol, x)
    if s == :p
        setfield!(d, s, x)
        return nothing
    elseif s == :output
        f = libgb.GrB_OUTP
    elseif s == :mask
        f = libgb.GrB_MASK
    elseif s == :input1
        f = libgb.GrB_INP0
    elseif s == :input2
        f = libgb.GrB_INP1
    #elseif s == :nthreads
        #f = libgb.GxB_DESCRIPTOR_NTHREADS
    #elseif s == :chunk
        #f = libgb.GxB_DESCRIPTOR_CHUNK
    #elseif s == :axb_method
        #f = libgb.GxB_AxB_METHOD
    #elseif s == :sort
        #f = libgb.GxB_SORT
    else
        setfield!(d, s, x)
        return nothing
    end
    libgb.GrB_Descriptor_set(d, f, x)
end

function Base.:+(d1::Descriptor, d2::Descriptor)
    d = Descriptor()
    for f ∈ propertynames(d)
        if f == :input1
            if getproperty(d1, f) == TRANSPOSE && getproperty(d2, f) == TRANSPOSE
                setproperty!(d, f, DEFAULT)
            elseif getproperty(d1, f) == TRANSPOSE || getproperty(d2, f) == TRANSPOSE
                setproperty!(d, f, TRANSPOSE)
            end
        elseif f == :input2
            if getproperty(d1, f) == TRANSPOSE && getproperty(d2, f) == TRANSPOSE
                setproperty!(d, f, DEFAULT)
            elseif getproperty(d1, f) == TRANSPOSE || getproperty(d2, f) == TRANSPOSE
                setproperty!(d, f, TRANSPOSE)
            end
        else
            if getproperty(d1, f) != DEFAULT
                setproperty!(d, f, getproperty(d1, f))
            end
            if getproperty(d2, f) != DEFAULT
                setproperty!(d, f, getproperty(d2, f))
            end
        end
    end
    return d
end

Base.:+(d1::Descriptor, ::Nothing) = d1
Base.:+(::Nothing, d2::Descriptor) = d2
Base.:+(f1::libgb.GrB_Desc_Value, f2::libgb.GrB_Desc_Value) = libgb.GrB_Desc_Value(UInt32(f1) + UInt32(f2))
function Base.propertynames(::Descriptor)
    return (
    :output,
    :mask,
    :input1,
    :input2,
    #:nthreads,
    #:chunk,
    #:axb_method,
    #:sort,
    )
end

function _loaddescriptors()
    T1.p = load_global("GrB_DESC_T1", libgb.GrB_Descriptor)
    T0.p = load_global("GrB_DESC_T0", libgb.GrB_Descriptor)
    T0T1.p = load_global("GrB_DESC_T0T1", libgb.GrB_Descriptor)
    C.p = load_global("GrB_DESC_C", libgb.GrB_Descriptor)
    CT1.p = load_global("GrB_DESC_CT1", libgb.GrB_Descriptor)
    CT0.p = load_global("GrB_DESC_CT0", libgb.GrB_Descriptor)
    CT0T1.p = load_global("GrB_DESC_CT0T1", libgb.GrB_Descriptor)
    S.p = load_global("GrB_DESC_S", libgb.GrB_Descriptor)
    ST1.p = load_global("GrB_DESC_ST1", libgb.GrB_Descriptor)
    ST0.p = load_global("GrB_DESC_ST0", libgb.GrB_Descriptor)
    ST0T1.p = load_global("GrB_DESC_ST0T1", libgb.GrB_Descriptor)
    SC.p = load_global("GrB_DESC_SC", libgb.GrB_Descriptor)
    SCT1.p = load_global("GrB_DESC_SCT1", libgb.GrB_Descriptor)
    SCT0.p = load_global("GrB_DESC_SCT0", libgb.GrB_Descriptor)
    SCT0T1.p = load_global("GrB_DESC_SCT0T1", libgb.GrB_Descriptor)
    R.p = load_global("GrB_DESC_R", libgb.GrB_Descriptor)
    RT1.p = load_global("GrB_DESC_RT1", libgb.GrB_Descriptor)
    RT0.p = load_global("GrB_DESC_RT0", libgb.GrB_Descriptor)
    RT0T1.p = load_global("GrB_DESC_RT0T1", libgb.GrB_Descriptor)
    RC.p = load_global("GrB_DESC_RC", libgb.GrB_Descriptor)
    RCT1.p = load_global("GrB_DESC_RCT1", libgb.GrB_Descriptor)
    RCT0.p = load_global("GrB_DESC_RCT0", libgb.GrB_Descriptor)
    RCT0T1.p = load_global("GrB_DESC_RCT0T1", libgb.GrB_Descriptor)
    RS.p = load_global("GrB_DESC_RS", libgb.GrB_Descriptor)
    RST1.p = load_global("GrB_DESC_RST1", libgb.GrB_Descriptor)
    RST0.p = load_global("GrB_DESC_RST0", libgb.GrB_Descriptor)
    RST0T1.p = load_global("GrB_DESC_RST0T1", libgb.GrB_Descriptor)
    RSC.p = load_global("GrB_DESC_RSC", libgb.GrB_Descriptor)
    RSCT1.p = load_global("GrB_DESC_RSCT1", libgb.GrB_Descriptor)
    RSCT0.p = load_global("GrB_DESC_RSCT0", libgb.GrB_Descriptor)
    RSCT0T1.p = load_global("GrB_DESC_RSCT0T1", libgb.GrB_Descriptor)
    DEFAULTDESC.p = libgb.GrB_Descriptor_new()
end

Base.show(io::IO, ::MIME"text/plain", d::Descriptor) = gxbprint(io, d)
Base.print(io::IO, d::Descriptor) = gxbprint(io, d)


text:
@testset "Result" begin
    @testset "Result+PathResult" begin
        R = solve(equations(heart()), save_all_paths=true, show_progress=false)
        @test R isa Result
        @test natinfinity(R) ≤ 572
        @test nfinite(R) == 4
        @test length(collect(R)) == 576
        @test finite(R) isa Vector{<:PathResult}

        @test length(finite(R, onlynonsingular=false)) == 4
        @test length(finite(R, onlynonsingular=true)) == 4
        @test length(finite(R, onlysingular=true)) == 0
        @test 572 - length(failed(R)) == natinfinity(R)
        @test length(real(R, tol=1e-6)) == 2
        @test nreal(R, tol=1e-6) == 2
        @test length(atinfinity(R)) ≤ 572
        @test length(results(R, onlyreal=true, realtol=1e-8)) == 2
        @test length(results(R, onlynonsingular=true, singulartol=1e9)) == 4
        @test length(finite(results(R, onlyreal=true))) == 2
        @test nresults(R, onlynonsingular=true, singulartol=1e9) == 4
        @test length(results(R, onlyfinite=false)) == 576
        @test nresults(R, onlyfinite=false) == 576
        @test nnonsingular(R) == 4
        @test length(nonsingular(R)) == 4
        @test length(singular(R)) == 0
        @test length(singular(real(R))) == 0

        @test_nowarn mapresults(solution, R)
        @test_nowarn mapresults(start_solution, R)
        @test_nowarn mapresults(residual, R)
        @test_nowarn mapresults(issingular, R)
        @test count(mapresults(isreal, R)) == 2
        # test fallback
        @test count(results(isreal, R)) == 2

        @test length(solutions(R)) == 4
        @test solutions(R) isa Vector{Vector{ComplexF64}}
        @test realsolutions(R, realtol=1e-8) isa Vector{Vector{Float64}}
        @test length(realsolutions(R, realtol=1e-8)) == 2

        @test_nowarn string(R)
        @test_nowarn string(R[end])
        @test_nowarn show(IOContext(devnull, :compact=>true), R)
        @test_nowarn show(IOContext(devnull, :compact=>true), R[end])

        test_treeviews(R)

        @polyvar x y z
        R = solve([(x-3)^3,(y-2)], affine_tracking=true)
        @test R isa Result
        test_treeviews(R)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 3
        @test_nowarn sprint(show, R)

        R = solve([(x-3)^3,(y-2)], affine_tracking=true, system_scaling=nothing)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 3
        @test_nowarn sprint(show, R)

        @polyvar x y z
        R = solve([(x-3z),(y-2z)])
        @test R isa Result{<:ProjectiveVectors.PVector}
        @test_nowarn string(R)
        test_treeviews(R)
        @test nnonsingular(R) == 1

        @polyvar x y z
        R = solve([(x-3z)^3,(y-2z)])
        @test R isa Result{<:ProjectiveVectors.PVector}
        test_treeviews(R)
        @test nnonsingular(R) == 0
        @test nsingular(R) == 3

        @polyvar x y z
        V = x^2 + y^2 - z^2
        L = randn(1,3) * [x, y, z]
        R = solve([V; L])
        @test_nowarn sprint(show, R)
    end

    @testset "multiplicities" begin
        @polyvar x y z
        f = (x-1*y)*(x-1.01*y)
        g = x - y - z
        S = solve([f,g], threading=false)
        @test length(multiplicities(S, tol=1e-1)) == 1
        @test length(multiplicities(S, tol=1e-5)) == 0
        @test length(multiplicities(S.pathresults, tol=1e-1)) == 1
        @test length(multiplicities(S.pathresults, tol=1e-5)) == 0

        @polyvar x z
        f = (x-1)*(x-1.01)
        g = x - 1 - z
        S = solve([f,g], threading=false)
        @test length(multiplicities(S, tol=1e-1)) == 1
        @test length(multiplicities(S, tol=1e-5)) == 0
        @test length(multiplicities(S.pathresults, tol=1e-1)) == 1
        @test length(multiplicities(S.pathresults, tol=1e-5)) == 0
        @test norm(solution(S[1]) - solution(S[2])) < 1e-1
    end

    @testset "uniquesolutions" begin
        @polyvar x
        f = (x-3)^3*(x-2)
        R = solve([f])
        @test length(uniquesolutions(R)) == 2
        @test uniquesolutions(R) isa Vector{Vector{Complex{Float64}}}
        @test length(uniquesolutions(R, multiplicities=true)) == 2
        a, b = uniquesolutions(R, multiplicities=true)
        @test (a[2] == 3 && b[2] == 1) || (b[2] == 3 && a[2] == 1)

        @polyvar x y
        f = (x-3y)^3*(x-2y)
        R = solve([f])
        @test length(uniquesolutions(R)) == 2
        @test uniquesolutions(R) isa Vector{ProjectiveVectors.PVector{Complex{Float64}, 1}}
        @test length(uniquesolutions(R, multiplicities=true)) == 2
        a, b = uniquesolutions(R, multiplicities=true)
        @test (a[2] == 3 && b[2] == 1) || (b[2] == 3 && a[2] == 1)
    end
end


text:
# Merge with EventTraceData
printstyled(stdout,"  merge! on EventTraceData\n", color=:light_green)

# Test case where we have to merge phase catalogs
(S,T) = mktestseis()
S = convert(EventTraceData, S)
T = convert(EventTraceData, T)
for i = 1:S.n
  S.az[i] = (rand()-0.5)*360.0
  S.baz[i] = (rand()-0.5)*360.0
  S.dist[i] = (rand()-0.5)*360.0
  S.pha[i] = randPhaseCat()
end

for i = 1:T.n
  T.az[i] = (rand()-0.5)*360.0
  T.baz[i] = (rand()-0.5)*360.0
  T.dist[i] = (rand()-0.5)*360.0
end

# Force a phase pick mismatch
P_true = SeisPha(rand(Float64,8)..., 'U', '0')
P_false = SeisPha(rand(Float64,8)..., '+', '1')
i = findid("CC.LON..BHZ", S)
j = findid("CC.LON..BHZ", T)

# T always begins later and should thus preserve P_true in the merge
S.pha[i]["P"] = deepcopy(P_false)
T.pha[j]["P"] = deepcopy(P_true)
merge!(S,T)
i = findid("CC.LON..BHZ", S)
@test S.pha[i]["P"] == P_true

# Check that purge works
V = purge(S)
purge!(S)
@test S == V

printstyled(stdout,"  merge! extensions to EventTraceData, EventChannel\n", color=:light_green)
S = convert(EventTraceData, randSeisData())
T = deepcopy(S)
merge!(S)
@test merge(T) == S
merge!(S,T)
sort!(T)
@test S == T

C = convert(EventChannel, randSeisChannel())
T = merge(S, C)
merge!(S, C)
@test S == T

C = convert(EventChannel, randSeisChannel())
S = convert(EventTraceData, randSeisData())
@test merge(C, S) == merge(S, C)

A = EventTraceData[convert(EventTraceData, randSeisData()),
convert(EventTraceData, randSeisData()),
convert(EventTraceData, randSeisData())]
merge(A)

S = convert(EventTraceData, randSeisData())
T = convert(EventTraceData, randSeisData())
@test S*T == T*S

@test S*C == merge(S, EventTraceData(C))

C = convert(EventChannel, randSeisChannel())
D = convert(EventChannel, randSeisChannel())
S = merge(C,D)
@test typeof(S) == EventTraceData
@test C*D == S

# Check that dist, az, baz merge correctly
(S,T) = mktestseis()
S = convert(EventTraceData, S)
T = convert(EventTraceData, T)
T.az[3:4] .= 0.0
T.baz[3:4] .= 0.0
T.dist[3:4] .= 0.0
U = deepcopy(T)

for i = 1:S.n
  S.az[i] = (rand()-0.5)*360.0
  S.baz[i] = (rand()-0.5)*360.0
  S.dist[i] = (rand()-0.5)*360.0
end
merge!(T,S)
for i in 1:S.n
  n = findid(S.id[i], T)
  if n > 0
    @test S.az[i] == T.az[n]
    @test S.baz[i] == T.baz[n]
    @test S.dist[i] == T.dist[n]
  end
end


text:
using Test
using QBase

@testset "./src/measure.jl" begin

@testset "measure()" begin
    @testset "State" begin
        ρ = State([0.5 0.5im;-0.5im 0.5])
        Π = POVM([[1 0;0 0],[0 0;0 1]])

        marginals = measure(Π, ρ)

        @test marginals isa Probabilities
        @test marginals == [0.5,0.5]
    end

    @testset "Ket" begin
        ψ = Ket(1/sqrt(2)*[1,im])
        Π = POVM([[0.5 0.5im;-0.5im 0.5],[0.5 -0.5im;0.5im 0.5]])

        marginals = measure(Π, ψ)

        @test marginals isa Probabilities
        @test marginals ≈ [0,1]
    end

    @testset "Vector{State}" begin
        @testset "trine_measurements" begin
           conditionals = measure(
               mirror_symmetric_qubit_3povm(π/3), trine_qubit_states()
           )
           @test conditionals ≈ [2/3 1/6 1/6;1/6 2/3 1/6;1/6 1/6 2/3]
           @test conditionals isa Conditionals
       end

       @testset "classical_measurements" begin
           Π = POVM([[1 0;0 0],[0 0;0 1],[0 0;0 0]])
           ρ_set = State.([[0.5 0.5;0.5 0.5],[1 0;0 0],[0 0;0 1]])

           @test measure(Π, ρ_set) == [0.5 1 0;0.5 0 1;0 0 0]

           Π = POVM([[1 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]])
           ρ_set = State.([[1 0 0;0 0 0;0 0 0],[0 0 0;0 1 0;0 0 0],[0 0 0;0 0 0;0 0 1]])

           @test measure(Π, ρ_set) == [1 1 0;0 0 1]
       end
    end

    @testset "Vector{Ket}" begin
        ψ_kets = Ket.([[1,0],[0,1]])
        Π = POVM([[0 0;0 1],[1 0;0 0],[0 0;0 0]])

        conditionals = measure(Π, ψ_kets)

        @test conditionals isa Conditionals
        @test conditionals == [0 1;1 0;0 0]
    end
end

end