text:
using MultiplesOfPi
using Test
using Aqua

@testset "project quality" begin
    Aqua.test_all(MultiplesOfPi, ambiguities=false)
end

@testset "Constructor" begin
    @test (@inferred (1)*Pi) ≈ π
    @test (@inferred (1.5)*Pi) ≈ 1.5*pi
    @test (@inferred (1//2)*Pi) ≈ π/2

    @test (1)*Pi === Pi
    @test Pi == π
    @test π == Pi
    @test 2π ≈ 2Pi

    @test float(Pi^0.5) == pi^0.5

    @testset "nested" begin
        @test PiExpTimes{PiExpTimes{Int}}(3, 0) == 3
        @test PiExpTimes{PiExpTimes{Float64}}(3, 0) == PiExpTimes{Float64}(3)
        @test PiExpTimes{PiExpTimes{Float64}}(Pi, 0) == Pi

        @test PiExpTimes{PiExpTimes}(Pi) isa PiExpTimes{PiExpTimes}
        @test PiExpTimes{PiExpTimes}(Pi) == Pi
        @test PiExpTimes{PiExpTimes{Int}}(Pi) isa PiExpTimes{PiExpTimes{Int}}
        @test PiExpTimes{PiExpTimes{Int}}(Pi) == Pi
    end

    @testset "Irrational{:π}" begin
        @test (π)*Pi == PiExpTimes(1, 2)
        @test PiExpTimes{Irrational{:π}}(π) isa PiExpTimes{Irrational{:π}}
        @test PiExpTimes{Irrational{:π}}(π).x === π

        @test PiExpTimes{Real}(π) isa PiExpTimes{Real}
        @test PiExpTimes{Real}(π).x === 1

        @test PiExpTimes{Irrational{:π}}(π, 2) == PiExpTimes(1, 3)
    end

    @test PiExpTimes(1, -1) ≈ 1/π
    @test PiExpTimes(2, -2) ≈ 2/π^2

    @test PiExpTimes(Pi^2, 2) === Pi^4
    @test PiExpTimes(Pi^2, -2) == 1
    @test PiExpTimes(1.0Pi^2, -2) == 1.0
    @test PiExpTimes(Pi^-2, 2) == 1.0
    @test PiExpTimes{Int}(Pi^2, 2) === Pi^4
    @test PiExpTimes{Float64}(Pi^2, 2) === PiExpTimes(1.0, 4)

    @test PiExpTimes{Int}(2.0, 2) === PiExpTimes(2, 2)

    @test (pi)*Pi === PiExpTimes(1, 2)
    @test PiExpTimes(pi, 2) === PiExpTimes(1, 3)

    @test PiExpTimes(π, 2) === PiExpTimes(1, 3)

    @test PiExpTimes(PiExpTimes(1, 3), 2) === Pi^5

    @test float(PiExpTimes(2, 2)) === 2*pi*pi

    @test float(PiExpTimes(BigFloat(2)^10000)) isa BigFloat
end

@testset "Ordering" begin
    @testset "inequality" begin
        @testset "PiExpTimes" begin
            @test PiExpTimes(3, 2) < PiExpTimes(4, 2)
            @test PiExpTimes(3, 1) < PiExpTimes(3, 2)
            @test PiExpTimes(3, 1) < PiExpTimes(3.5, 2)
            @test (3)*Pi < PiExpTimes(3, 2)
            @test π < PiExpTimes(1, 2)
            @test !(PiExpTimes(1, 2) < π)
            @test 1 < PiExpTimes(1, 2)
            @test PiExpTimes(1, 2) > π
            @test PiExpTimes(1, 2) > 1
            @test PiExpTimes(1, 0) < 2
            @test PiExpTimes(1, 0) < π
        end
    end
    @testset "equality" begin
        @testset "PiTimes" begin
            @test Pi == π
            @test π == Pi
            @test (2)*Pi == (2)*Pi
            @test (0)*Pi == 0
            @test 0 == (0)*Pi
            @test Pi != 2Pi
            @test Pi != 4
        end
        @testset "PiExpTimes" begin
            Pi² = PiExpTimes(1, 2)
            Pi³ = PiExpTimes(1, 3)

            z² = PiExpTimes(0, 2)
            z³ = PiExpTimes(0, 3)

            @test Pi² == Pi²

            @test Pi² != π^2
            @test π^2 != Pi²

            @test Pi² != 3
            @test 3 != Pi²

            @test Pi² != Pi
            @test Pi  != Pi²

            @test Pi² != Pi³
            @test Pi³ != Pi²

            @test !(Pi² == π)
            @test !(π == Pi²)

            @test z² == z³
            @test z² == 0
        end
        @testset "nested" begin
            @test π == PiExpTimes{PiExpTimes{Int}}(Pi^2, -1)
            @test PiExpTimes{PiExpTimes{Int}}(Pi^2, -1) == π
        end
    end
end

@testset "sign" begin
    @testset "PiExpTimes" begin
        @testset "sign" begin
            @test sign(PiExpTimes(1, 2)) == 1
            @test sign(PiExpTimes(3.4, 2)) == 1.0
            @test sign(PiExpTimes(0.0, 2)) == 0.0
            @test sign(PiExpTimes(-0.0, 2)) == -0.0
            @test sign(PiExpTimes(1//2, 2)) == 1//1
        end
        @testset "signbit" begin
            @test !signbit(PiExpTimes(1, 2))
            @test !signbit(PiExpTimes(0, 2))
            @test !signbit(PiExpTimes(-0, 2))
            @test signbit(PiExpTimes(-1, 2))
            @test signbit(PiExpTimes(-0.0, 2))
            @test !signbit(PiExpTimes(0.0, 2))
        end
    end
end

@testset "zero and one" begin
    @testset "one" begin
        for T in (Float16,Float32,Float64,Int8,Int16,Int32,Int64,Int128,BigInt,BigFloat)
            @test one(PiExpTimes{T,Int}) == 1
            @test one(PiExpTimes{T,Float64}) == 1
            @test one(PiExpTimes{T,T}) == 1
            @test one(PiExpTimes(zero(T))) == 1
        end

        @test one(Pi) * Pi == Pi
        @test one(Pi^2) * Pi^2 == Pi^2

        @test one(PiExpTimes{PiExpTimes{Int}}(1)) == 1
    end
    @testset "zero" begin
        @test zero(Pi) == 0
        @test zero(typeof(Pi)) == 0

        @test zero(Pi) + Pi == Pi
        @test zero(Pi^2) + Pi^2 == Pi^2

        for T in (Float16,Float32,Float64,Int8,Int16,Int32,Int64,Int128)
            @test zero(PiExpTimes{T,Int}) isa PiExpTimes{T,Int}
            @test zero(PiExpTimes{T,Float64}) isa PiExpTimes{T,Float64}
            @test zero(PiExpTimes{T,T}) isa PiExpTimes{T,T}
            @test zero(PiExpTimes{T}(1)) isa PiExpTimes{T}
            @test zero(PiExpTimes{T}(1)) == 0
        end
        for T in (BigInt,BigFloat)
            @test zero(PiExpTimes{T,Int}) isa PiExpTimes{T,Int}
            @test zero(PiExpTimes{T,BigInt}) isa PiExpTimes{T,BigInt}
            @test zero(PiExpTimes{T}(1)) isa PiExpTimes{T}
            @test zero(PiExpTimes{T}(1)) == 0
        end
    end
end

@testset "Range" begin
    @testset "LinRange" begin
        l = LinRange(0.0Pi,1.0Pi,11)
        lpi = LinRange(0,π,11)
        @test l[1] === 0.0Pi
        @test l[end] === 1.0Pi

        for i in eachindex(l,lpi)
            @test l[i] ≈ lpi[i]
        end
    end
    @testset "range" begin
        @testset "step" begin
            r = range(0Pi, step=Pi, length = 3)
            @test length(r) == 3
            @test r[1] == 0Pi
            @test r[end] == 2Pi
            @test step(r) == 1Pi
        end
    end
end

@testset "arithmetic" begin
    p = 3Pi
    q = 4Pi
    r = 4Pi^2

    @testset "negation" begin
        @test -p isa PiExpTimes{Int}
        @test -r isa PiExpTimes{Int}
        @test -p === -p.x*Pi
        @test -r === PiExpTimes(-r.x, 2)
    end

    @testset "addition" begin
        @test p + p === (2p.x)*Pi
        @test p + q === (p.x + q.x)*Pi
        @test r + r === PiExpTimes(2r.x, 2)
        @test p + r === (p.x + r.x*π)*Pi
        @test p + r ≈ float(p) + float(r)
    end

    @testset "subtraction" begin
        @test p - p === (0)*Pi == 0
        @test p - q === (p.x - q.x)*Pi
        @test r - r === PiExpTimes(zero(r.x), 2) == 0
        @test p - r === (p.x - r.x*π)*Pi
        @test p - r ≈ float(p) - float(r)
    end

    @testset "multiplication" begin
        @testset "PiTimes" begin
            @test (@inferred p*p) === PiExpTimes(p.x^2, 2)
            @test (@inferred p*p) ≈ (p.x*π)^2
            @test (@inferred p*q) === PiExpTimes(p.x*q.x, 2)
            @test (@inferred q*p) === PiExpTimes(p.x*q.x, 2)
            @test (@inferred p*q) ≈ (p.x*q.x)*π^2
            @test (@inferred p*r) === PiExpTimes(p.x*r.x, 3) ≈ (p.x*r.x)*π^3
            @test (@inferred r*p) === PiExpTimes(p.x*r.x, 3) ≈ (p.x*r.x)*π^3
        end
        @testset "Real" begin
            @testset "Int" begin
                @test p*2 === (p.x*2)*Pi
                @test 2*p === (p.x*2)*Pi
                @test 2p === (p.x*2)*Pi
                @test 2r === PiExpTimes(2r.x, 2)

                @test PiExpTimes(3, 0) * 2 == 6
                @test 2 * PiExpTimes(3, 0) == 6
            end
            @testset "Bool" begin
                @test (@inferred true*p) === p
                @test (@inferred p*true) === p
                @test (@inferred false*p) === zero(p)
                @test (@inferred p*false) === zero(p)

                @test PiExpTimes(3, 0) * true == 3
                @test PiExpTimes(3, 0) * false == 0
                @test PiExpTimes(3.0, 0) * true == 3.0
                @test PiExpTimes(3.0, 0) * false == 0.0
                @test PiExpTimes(-3.0, 0) * true == -3.0
                @test PiExpTimes(-3.0, 0) * false == -0.0
                @test sign(PiExpTimes(-3.0, 0) * false) === -0.0
                @test -3.0Pi * false === false * (-3.0Pi)
            end
            @testset "float" begin
                @test (@inferred p*1/p) == 1.0
                @test (@inferred r*1/r) == 1.0

                @test 2.3*p === (2.3*p.x)*Pi

                @test PiExpTimes(3.0, 0) * 2 == 6.0
                @test 2 * PiExpTimes(3.0, 0) == 6.0

                @test 2 * PiExpTimes(-3.0, 0) == -6.0
            end
        end

        @testset "Complex" begin
            @testset "Complex{Bool}" begin
                @test (@inferred p * im) === (0)*Pi + (p.x)*Pi*im
                @test (@inferred im* p) === (0)*Pi + (p.x)*Pi*im
                @test Complex(true,false)*Pi === Pi + im*zero(Pi)
                @test Complex(false,true)*Pi === zero(Pi) + im*Pi
                @test Pi*Complex(false,true) === zero(Pi) + im*Pi
            end
            @testset "Complex{Int}" begin
                @test (1+im)*p === (p.x)*Pi + (p.x)*Pi*im
                @test p*(1+im) === (p.x)*Pi + (p.x)*Pi*im
            end
        end
    end

    @testset "division" begin
        @testset "PiTimes" begin
            @test p/p == p.x/p.x
            @test p/q == p.x/q.x == 3/4 # works as 1/4 can be stored exactly
            @test q/p == q.x/p.x
        end
        @testset "PiExpTimes" begin
            @test r/r == 1.0
            @test r/p == (r.x/p.x)*Pi
            @test p/r == PiExpTimes(p.x/r.x, -1)
        end
        @testset "Real" begin
            @testset "Int" begin
                @test p/2 === (p.x/2)*Pi
                @test (3)*Pi/4 === (3/4)*Pi
            end
            @testset "Bool" begin
                @test p/true == (p.x/true)*Pi
                @test p/false == (Inf)*Pi
                @test true/Pi == PiExpTimes(1.0, -1)
                @test false/Pi == PiExpTimes(0.0, -1)
            end
        end
        @testset "pi" begin
            @test p/π == p.x
            @test (7.3)*Pi/π == 7.3
            @test π/((1/7.3)*Pi) ≈ 7.3 # within rounding errors
            @test π/((1//7)*Pi) == 7//1

            @test Pi^2/π == Pi
            @test π/Pi^2 == Pi^-1
        end
        @testset "Complex" begin
            @testset "Complex{Bool}" begin
                @test p/im === Complex(zero(Float64)*Pi,(-p.x/1)*Pi)
                @test p/Complex(true,false) == Complex((p.x/1),-zero(Float64))*Pi
                @test im/Pi === im*PiExpTimes(1.0, -1)
                @test Complex(true,false)/Pi == Complex(true,false)*PiExpTimes(1.0, -1)
                @test Pi*Pi/im == -1.0Pi*Pi*im
            end
            @testset "Complex{Int}" begin
                @test p/(1+im) === Complex(1.5Pi,-1.5Pi)
                @test Pi/2im === Complex(0.0Pi,-0.5Pi)
                @test (1+im)/Pi === (1+im)*PiExpTimes(1.0, -1)
            end
            @testset "Complex{PiTimes}" begin
                @testset "Real/Complex{PiExpTimes}" begin
                    @test Pi*Pi/(im*Pi*Pi) == Complex(0.0,-1.0)
                    resim = -PiExpTimes(1.0, -1)
                    @test Pi/(im*Pi*Pi) == Complex(zero(resim),resim)
                end

                @testset "Complex{PiExpTimes}/Real" begin
                    @testset "Complex{PiExpTimes}/Int" begin
                        @test Pi*Pi*im/1 == Complex(0.0*Pi^2,1.0*Pi^2)
                    end
                    @testset "Complex{PiExpTimes}/PiExpTimes" begin
                        @test (Pi*im)/Pi == 1.0im
                        @test Pi*Pi*im/Pi == 1.0Pi*im
                        @test Pi*Pi*im/(Pi*Pi) == Complex(0.0,1.0)
                    end
                end

                @testset "Complex{PiExpTimes}/Complex" begin
                    @testset "Complex{PiExpTimes}/Complex{Int}" begin
                        z = Pi*(1+2im)/(1+2im)
                        @test real(z) == (1.0)*Pi
                        @test iszero(imag(z))
                    end
                    @testset "Complex{PiExpTimes}/Complex{PiExpTimes}" begin
                        z = Pi*(1+2im)
                        @test isone(real(z/z))
                        @test iszero(imag(z/z))

                        @test Pi^2*im/(im*Pi^2) == Complex(1.0,0.0)
                        @test Pi^2*im/(im*Pi) == Complex(1.0,0.0)*Pi
                        @test Pi*im/(im*Pi^2) == Complex(1.0,0.0)*Pi^-1
                    end
                end
            end
        end
    end

    @testset "exponentiation" begin
        @testset "Int exponent" begin
            @test Pi^0 === PiExpTimes(1, 0) == 1
            @test Pi^1 === Pi == π
            @test Pi^2 === PiExpTimes(1, 2) ≈ π*π
            @test Pi^-1 === 1/Pi === PiExpTimes(1.0, -1) ≈ 1/π
            @test Pi^-2 === 1/Pi^2 === PiExpTimes(1.0, -2) ≈ 1/π^2

            for n=1:10
                @test Pi^n ≈ π^n
                @test Pi^-n == 1/Pi^n ≈ float(π)^-n
                @test float(Pi^-n) === float(1/Pi^n) ≈ 1/π^n
            end

            @test (Pi^2)^0 == 1

            @test float(Pi^big(10000)) == BigFloat(pi)^10000
        end

        @testset "Float64 exponent" begin
            @test float(@inferred Pi^2.0) === pi^2.0
            @test float(@inferred Pi^-2.0) === pi^-2.0
            @test (@inferred Pi^0.0) == 1.0
            @test float(@inferred Pi^0.5) == float(@inferred sqrt(Pi))
            @test BigFloat(@inferred Pi^0.5) == BigFloat(@inferred sqrt(Pi))
            @test float(@inferred Pi^(1//3)) == float(@inferred cbrt(Pi))
            @test BigFloat(@inferred Pi^(1//3)) == BigFloat(@inferred cbrt(Pi))
            for T in [Float64, BigFloat]
                @test T(sqrt(Pi)^2) == T(Pi)
                @test T(sqrt(Pi^2)) == T(Pi)
                @test T(sqrt(Pi^2.0)) == T(Pi)
                @test T(sqrt(Pi^(2//1))) == T(Pi)
                @test T(cbrt(Pi)^3) == T(Pi)
                @test T(cbrt(Pi^3)) == T(Pi)
                @test T(cbrt(Pi^3.0)) == T(Pi)
                @test T(cbrt(Pi^(3//1))) == T(Pi)
            end
        end

        @testset "Irrational exponent" begin
            for T in [Float64, BigFloat]
                @test T(@inferred Pi^Pi) == T(pi)^pi
                @test T(@inferred Pi^π) == T(pi)^pi
                @test T(@inferred π^Pi) == T(pi)^pi
            end
        end

        @testset "Rational exponent" begin
            @test Float64(@inferred Pi^(3//4)) == pi^(3//4)
        end

        @testset "Pi as exponent" begin
            @test (@inferred 1^Pi) == 1
            @test (@inferred 2^Pi) == 2^π
            @test (@inferred 2.0^Pi) == 2.0^π
            @test (@inferred 2.0^(1.0Pi)) == 2.0^(1.0*π)
            @test (@inferred big(2)^Pi) == big(2)^π
            @test (@inferred big(2.0)^Pi) == big(2.0)^π
        end

        @testset "reciprocal" begin
            @test (@inferred Pi * Pi^-1)== 1.0
            @test (@inferred Pi^-1 * Pi)== 1.0
            @test (@inferred Pi^-1.0 * Pi)== 1.0
            @test (@inferred 2Pi * Pi^-1)== 2.0
            @test (@inferred 2Pi^2 * Pi^-2) == 2.0
            @test (@inferred 2Pi^2.0 * Pi^-2) == 2.0
        end
    end

    @testset "inv" begin
        @test inv(Pi) === PiExpTimes(1.0, -1)
        @test inv(Pi^2) === PiExpTimes(1.0, -2)
        @test inv(PiExpTimes{Int}(1, 0)) == 1.0
    end

    @testset "Rational" begin
        @test p//2 === (p.x//2)*Pi
        @test Pi//2 === (1//2)*Pi
        @test Pi//Pi == 1//1
        @test 3Pi//4Pi == 3//4
        @test (1//2)Pi//3Pi == 1//6
        @test 2Pi//(3//2)Pi == 4//3
        @test (1//2)Pi//(1//3)Pi == 3//2
        @test (1//2 + 3//2 * im)Pi == (1+3im)*Pi // 2

        @test PiExpTimes{Int}(1, 0)//2 == 1//2

        @testset "pi" begin
            @test π//Pi == 1//1
            @test Pi//π == 1//1

            @test π//Pi^2 == (1//Pi)
            @test Pi^2//π == Pi//1
            @test PiExpTimes{Int}(1, 0)//Pi == 1//Pi
        end

        @testset "PiExpTimes // PiExpTimes" begin
            @test Pi^2//Pi === Pi//1
            @test Pi//Pi^2 === 1//Pi
            @test (3Pi^2//2)//(Pi//2) == 3Pi//1
            @test (3Pi//2)//(Pi^2//2) == (3//Pi)
        end

        @testset "Real // Complex" begin
            @testset "Complex{Int}" begin
                @test Pi//(2im) === (0//1 - 1//2*im)*Pi
                @test Pi^2//(2im) === (0//1 - 1//2*im)*Pi^2
                @test Pi^-1//(2im) === (0//1 - 1//2*im)*Pi^-1
            end
            @testset "Complex{PiTimes}" begin
                @test 2Pi // (2Pi + 0Pi*im) == 1//1 + 0//1*im
                @test 2Pi^2 // (2Pi^2 + 0Pi^2*im) == 1//1 + 0//1*im
                @test 2Pi^2 // (2Pi + 0Pi*im) === (1//1)*Pi + (0//1)*Pi*im
                @test 2Pi // (2Pi^2 + 0Pi^2*im) === (1//Pi) + (0//Pi)*im
            end
        end

        @testset "Complex // Real" begin
            @test (2//1 + 1//1*im)Pi // 2 === (1//1 + 1//2*im)Pi
            @test (2//1 + 1//1*im)Pi^2 // 2 === (1//1 + 1//2*im)Pi^2
            @test (2//1 + 1//1*im)Pi // Pi == (2//1 + 1//1*im)
            @test (2//1 + 1//1*im)Pi^2 // Pi^2 == (2//1 + 1//1*im)
            @test (2//1 + 1//1*im)Pi^2 // Pi === (2//1 + 1//1*im)*Pi
            @test (2//1 + 1//1*im)Pi // Pi^2 === (2//Pi + (1//Pi)*im)
            @test Complex(π,π)//Pi == Complex(1//1,1//1)
            @test Complex(π,π)//Pi^2 === Complex(1//Pi,1//Pi)
        end

        @testset "Complex // Complex" begin
            z1 = (1//2 + 3//2 * im)Pi
            z2 = (1//2 + 3//2 * im)Pi^2
            @test z1 // z1 == (1//1 + 0//1 *im)
            @test z2 // z2 == (1//1 + 0//1 *im)
            @test z2 // z1 === (1//1 + 0//1 * im)*Pi
            @test z1 // z2 === (1 + 0*im)//Pi
        end
    end

    @testset "Irrational" begin
        @testset "pi" begin
            @test Pi + π == 2Pi
            @test Pi^1.0 + π == 2Pi
            @test π + Pi == 2Pi
            @test Pi - π == zero(Pi) == 0
            @test π - Pi == zero(Pi) == 0
            @test π*Pi == Pi*Pi
            @test Pi*π == Pi*Pi
            @test 2Pi*π == (2Pi)*Pi
            @test π*2Pi == (2Pi)*Pi
            @test Pi/π == 1
            @test π/Pi == 1.0
            @test (Pi)*Pi/π == Pi
            @test Pi*Pi/π == Pi
            @test Pi*π/Pi == (one(Float64))*Pi
            @test Pi*π/Pi == Pi == pi
        end
        @testset "e" begin
            @test Pi + ℯ == π + ℯ
            @test ℯ + Pi == ℯ + π
            @test Pi - ℯ == π - ℯ
            @test ℯ - Pi == ℯ - π
            @test ℯ*Pi == ℯ*π
            @test Pi*ℯ == π*ℯ
            @test Pi/ℯ == π/ℯ
            @test ℯ/Pi == ℯ/π
            @test float(Pi^ℯ) == pi^ℯ
        end
    end
end

@testset "trigonometric functions" begin
    @testset "PiTimes" begin
        @testset "sin and cos" begin
            x = 1.5
            @test sin(x*Pi) == sinpi(x) == -1 ≈ sin(π*x)
            x = 1
            @test cos(x*Pi) == cospi(x) == -1 ≈ cos(π*x)

            a = sincos(x*Pi)
            @test a[1] == sinpi(x)
            @test a[2] == cospi(x)
        end

        @testset "tan and cot" begin
            @testset "tan" begin
                @test isinf(tan(Pi/2))
                @test isinf(tan(-Pi/2))
                @test (@inferred tan(Pi/4)) == 1
                @test (@inferred tan(-Pi/4)) == -1
                @test (@inferred tan(Pi*0)) == 0
                @test (@inferred tan(-Pi*0)) == 0
                @test (@inferred tan(Pi)) == 0
                @test (@inferred tan(-Pi)) == 0
                @test (@inferred tan(2Pi)) == 0
                @test (@inferred tan(-2Pi)) == 0
                @test (@inferred tan(Pi/8)) == tan(pi/8)

                @test tan(nextfloat(0.0)*π) == tan(nextfloat(0.0)*Pi)
                @test tan(prevfloat(0.0)*π) == tan(prevfloat(0.0)*Pi)
            end
            @testset "cot" begin
                @test cot(Pi/2) == 0
                @test cot(-Pi/2) == 0
                @test cot(Pi/4) == 1
                @test cot(-Pi/4) == -1
                @test isinf(cot(Pi*0))
                @test isinf(cot(-Pi*0))
                @test isinf(cot(Pi))
                @test isinf(cot(-Pi))
                @test isinf(cot(2Pi))
                @test isinf(cot(-2Pi))
                @test cot(Pi/8) == cot(pi/8)

                @test cot(nextfloat(0.0)*π) == cot(nextfloat(0.0)*Pi)
                @test cot(prevfloat(0.0)*π) == cot(prevfloat(0.0)*Pi)
            end
        end

        @testset "complex exp" begin
            x = 1.1
            @test cis(x*Pi) == cospi(x) + im*sinpi(x) ≈ cos(π*x) + im*sin(π*x)

            @test exp(im*(0)*Pi) == cis((0)*Pi) == 1
            @test exp(im*(1/2)*Pi) == cis((1/2)*Pi) == im
            @test exp(im*(1)*Pi) == cis((1)*Pi)== -1
            @test exp(im*(3/2)*Pi) == cis((3/2)*Pi) == -im

            @test exp(im*Pi/2) == im
            @test exp(im*Pi) == -1
            @test exp(im*3Pi/2) == -im
        end

        @testset "sec and csc" begin
            @test sec((1)*Pi) == -1
            @test sec(Pi) == -1
            @test sec((2)*Pi) == 1
            @test sec(2Pi) == 1
            @test !isfinite(sec((1/2)*Pi))
            @test !isfinite(sec(Pi/2))

            @test csc((1/2)*Pi) == 1
            @test csc(Pi/2) == 1
            @test csc((3/2)*Pi) == -1
            @test csc(3Pi/2) == -1
            @test !isfinite(csc((1)*Pi))
            @test !isfinite(csc(Pi))
        end

        @testset "sinc" begin
            @test sinc(Pi) === sinc(π)
        end
    end
    @testset "PiExpTimes" begin
        @testset "Pi^2" begin
            @test sin(Pi^2) ≈ sin(π^2)
            @test cos(Pi^2) ≈ cos(π^2)
            s1, c1 = sincos(Pi^2)
            s2, c2 = sincos(π^2)
            @test (s1 ≈ s2) && (c1 ≈ c2)
            @test tan(Pi^2) ≈ tan(π^2)
            @test cot(Pi^2) ≈ cot(π^2)
            @test sinc(Pi^2) ≈ sinc(π^2)
            @test sinc(Pi^-2) ≈ sinc((1/π)^2)
            @test sinc(1/Pi^2) ≈ sinc((1/π)^2)
            @test sec(Pi^2) ≈ sec(π^2)
            @test csc(Pi^2) ≈ csc(π^2)
        end
        @testset "Pi^0" begin
            x = 1
            p = PiExpTimes{Int}(x, 0)
            @test sin(p) == sin(x)
            @test cos(p) == cos(x)
            @test sincos(p) == sincos(x)
            @test tan(p) == tan(x)
            @test cot(p) == cot(x)
            @test sinc(p) == sinc(x)
            @test sinc(p) == sinc(x)
            @test sec(p) == sec(x)
            @test csc(p) == csc(x)
        end
    end
end

@testset "hyperbolic" begin
    @testset "cosh" begin
        @test cosh(im*(1/2)*Pi) == 0
    end
    @testset "sinh" begin
        @test sinh(im*(-1)*Pi) == 0
    end

    @testset "tanh" begin
        @test tanh(im*Pi/4) == im
        @test tanh(im*Pi/2) == im*Inf
        @test tanh(im*Pi/3) == im*tan(Pi/3)
        @test tanh((1+im)*Pi) == tanh((1+im)*π)
        @test tanh(2+im*Pi) == tanh(2+im*π)
    end

    @testset "PiExpTimes" begin
        @testset "Pi^2" begin
            @test sinh(Pi^2) === sinh(π^2)
            @test cosh(Pi^2) === cosh(π^2)
            @test tanh(Pi^2) === tanh(π^2)
        end
        @testset "Pi^0" begin
            x = 1
            p = PiExpTimes{Int}(x, 0)
            @test sinh(p) === sinh(x)
            @test cosh(p) === cosh(x)
            @test tanh(p) === tanh(x)
        end
    end
end

@testset "Inf NaN zero check" begin
    @test isinf((Inf)*Pi)
    @test !isfinite((Inf)*Pi)
    @test isfinite(Pi)
    @test !isfinite(PiExpTimes(Inf, 2))
    @test !isfinite(PiExpTimes(Inf, 0))
    @test !isfinite(PiExpTimes(Inf, -1))
    @test isinf(PiExpTimes(Inf, 2))
    @test isinf(PiExpTimes(Inf, 0))
    @test isinf(PiExpTimes(Inf, -1))
    @test isinf(PiExpTimes(1,Inf))
    @test isfinite(PiExpTimes(1,-Inf))

    @test isnan((NaN)*Pi)
    @test isnan(PiExpTimes(1,NaN))
    @test isnan(PiExpTimes(NaN,1))
    @test isnan(PiExpTimes(NaN,NaN))
    @test isnan(PiExpTimes(0, Inf))

    @test iszero((0)*Pi)
    @test iszero(PiExpTimes(0, 2))
    @test iszero(PiExpTimes(0, 0))
    @test iszero(Pi^-Inf)
    @test !iszero(PiExpTimes(0, NaN))
end

@testset "show" begin
    @testset "Int" begin
        @test repr(Pi) == "Pi"
        @test repr(2Pi) == "2Pi"
        @test repr(0*Pi) == "0"
    end

    @testset "Rational" begin
        @test repr((1//1)*Pi) == "(1//1)Pi"
        @test repr((2//3)*Pi) == "(2//3)Pi"
        @test repr((2//3)*Pi*Pi) == "(2//3)Pi^2"
    end

    @testset "float" begin
        @test repr(1.0Pi) == "Pi"
        @test repr(2.0Pi) == "2.0*Pi"
        @test repr(0.0Pi) == "0.0"
        @test repr(-0.0Pi) == "-0.0"

        @testset "inf nan" begin
            @test repr(Inf*Pi) == "Inf"
            @test repr(-Inf*Pi) == "-Inf"
            @test repr(NaN*Pi) == "NaN"
        end
    end

    @testset "PiExpTimes" begin
        @test repr(Pi*Pi) == "Pi^2"
        @test repr(Pi*Pi*Pi) == "Pi^3"

        @testset "nested" begin
            @test repr(PiExpTimes{PiExpTimes}(1)) == "Pi^0"
        end
    end

    @testset "complex" begin
        @test repr(im*Pi) == "0 + Pi*im"

        @testset "PiExpTimes" begin
            @test repr(im*Pi*Pi) == "0 + Pi^2*im"
        end

        @testset "float" begin
            @test repr(im*1.2*Pi) == "0.0 + 1.2*Pi*im"
            @test repr((-0.0)*Pi*im) == "-0.0 - 0.0*im"
        end

        @testset "Rational" begin
            @test repr(Pi + im*(2//3)Pi) == "(1//1)Pi + (2//3)Pi*im"
            @test repr(Pi*Pi*(2//3)*im) == "0//1 + (2//3)Pi^2*im"
        end
    end

    @testset "non-finite" begin
        @test repr(Pi^NaN) == repr(NaN)
        @test repr(Pi^Inf) == repr(Inf)
    end
end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# enumeration.jl: replicate Enumeration_C Spinnaker SDK example with low-level library

using Spinnaker
import Spinnaker: bool8_t, MAX_BUFFER_LEN

# Get system interface and library version
hSystem = Ref(spinSystem(C_NULL))
spinSystemGetInstance(hSystem)

hLibraryVersion = Ref(spinLibraryVersion(0,0,0,0))
spinSystemGetLibraryVersion(hSystem[], hLibraryVersion)

libver = VersionNumber(hLibraryVersion[].major, hLibraryVersion[].minor, hLibraryVersion[].type)
@info "Spinnaker library version: $libver, build $(hLibraryVersion[].build)"

# Retrieve list of interfaces
hInterfaceList = Ref(spinInterfaceList(C_NULL))
nInterfaces = Ref(Csize_t(0))
spinInterfaceListCreateEmpty(hInterfaceList)
spinSystemGetInterfaces(hSystem[], hInterfaceList[])
spinInterfaceListGetSize(hInterfaceList[], nInterfaces)

@info "Discovered $(nInterfaces[]) interfaces"

# Retrieve camera list (from system object)
hCameraList = Ref(spinCameraList(C_NULL))
nCameras = Ref(Csize_t(0))
spinCameraListCreateEmpty(hCameraList)
spinSystemGetCameras(hSystem[], hCameraList[])
spinCameraListGetSize(hCameraList[], nCameras)

@info "Discovered $(nCameras[]) cameras"

# Examine a given interface
function query_interface(hInterface)

  # Get interface and print details
  hNodeMapInterface = Ref(spinNodeMapHandle(C_NULL))
  hInterfaceDisplayName = Ref(spinNodeHandle(C_NULL))
  interfaceDisplayNameIsAvailable = Ref(bool8_t(false))
  interfaceDisplayNameIsReadable = Ref(bool8_t(false))
  spinInterfaceGetTLNodeMap(hInterface, hNodeMapInterface);
  spinNodeMapGetNode(hNodeMapInterface[], "InterfaceDisplayName", hInterfaceDisplayName);
  spinNodeIsAvailable(hInterfaceDisplayName[], interfaceDisplayNameIsAvailable);
  spinNodeIsReadable(hInterfaceDisplayName[], interfaceDisplayNameIsReadable);

  interfaceDisplayName = Vector{UInt8}(undef, MAX_BUFFER_LEN)
  lenInterfaceDisplayName = Ref(Csize_t(MAX_BUFFER_LEN))
  if (interfaceDisplayNameIsAvailable[] == true) && (interfaceDisplayNameIsReadable[] == true)
    spinStringGetValue(hInterfaceDisplayName[], interfaceDisplayName, lenInterfaceDisplayName);
    @info "$(unsafe_string(pointer(interfaceDisplayName)))"
  else
    @info "Interface display name not readable"
  end

  hCameraList = Ref(spinCameraList(C_NULL))
  nCameras = Ref(Csize_t(0))
  spinCameraListCreateEmpty(hCameraList)
  spinInterfaceGetCameras(hInterface, hCameraList[]);
  spinCameraListGetSize(hCameraList[], nCameras);

  # Get cameras on interface nd print details
  if nCameras[] > 0

    for i in 0:nCameras[]-1

      hCam = Ref(spinCamera(C_NULL))
      spinCameraListGet(hCameraList[], i, hCam);

      hNodeMapTLDevice = Ref(spinNodeMapHandle(C_NULL))
      spinCameraGetTLDeviceNodeMap(hCam[], hNodeMapTLDevice);

      #  Get camera vendor name
      hDeviceVendorName = Ref(spinNodeHandle(C_NULL))
      deviceVendorNameIsAvailable = Ref(bool8_t(false))
      deviceVendorNameIsReadable = Ref(bool8_t(false))
      spinNodeMapGetNode(hNodeMapTLDevice[], "DeviceVendorName", hDeviceVendorName)
      spinNodeIsAvailable(hDeviceVendorName[], deviceVendorNameIsAvailable)
      spinNodeIsReadable(hDeviceVendorName[], deviceVendorNameIsReadable)

      deviceVendorName = Vector{UInt8}(undef, MAX_BUFFER_LEN)
      lenDeviceVendorName = Ref(Csize_t(MAX_BUFFER_LEN))
      if (deviceVendorNameIsAvailable[] == true) && (deviceVendorNameIsReadable[] == true)
        spinStringGetValue(hDeviceVendorName[], deviceVendorName, lenDeviceVendorName);
        @info "$(unsafe_string(pointer(deviceVendorName)))"
      else
        @info "Camera vendor name not readable"
      end

      # Get camera model name
      hDeviceModelName = Ref(spinNodeHandle(C_NULL))
      deviceModelNameIsAvailable = Ref(bool8_t(false))
      deviceModelNameIsReadable = Ref(bool8_t(false))
      spinNodeMapGetNode(hNodeMapTLDevice[], "DeviceModelName", hDeviceModelName);
      spinNodeIsAvailable(hDeviceModelName[], deviceModelNameIsAvailable);
      spinNodeIsReadable(hDeviceModelName[], deviceModelNameIsReadable);

      deviceModelName = Vector{UInt8}(undef, MAX_BUFFER_LEN)
      lenDeviceModelName = Ref(Csize_t(MAX_BUFFER_LEN))
      if (deviceModelNameIsAvailable[] == true) && (deviceModelNameIsReadable[] == true)
        spinStringGetValue(hDeviceModelName[], deviceModelName, lenDeviceModelName);
        @info "$(unsafe_string(pointer(deviceModelName)))"
      else
        @info "Camera model name not readable"
      end

      spinCameraRelease(hCam[]);

    end
  else
    @info "No cameras detected"
  end

  spinCameraListClear(hCameraList[]);
  spinCameraListDestroy(hCameraList[]);
  return nothing

end

# Examine cameras on each interface
if nCameras[] > 0
  for i in 0:nInterfaces[]-1
    hInterface = Ref(spinInterface(C_NULL))
    spinInterfaceListGet(hInterfaceList[], i, hInterface)
    query_interface(hInterface[])
    spinInterfaceRelease(hInterface[])
  end
else
  # Finish if there are no cameras
  @info "No cameras, finishing"
end

# Cleanup
spinCameraListClear(hCameraList[])
spinCameraListDestroy(hCameraList[])
spinInterfaceListClear(hInterfaceList[])
spinInterfaceListDestroy(hInterfaceList[])
spinSystemReleaseInstance(hSystem[])

ext: jl


text:
# Find strongly connected components
function strongconnected(cfg; blocks = 1:length(cfg))
  preorder = zeros(Int, length(cfg))
  S, P = Int[], Int[]
  C = Ref(0)
  components = Vector{Int}[]
  let # avoid boxing `search`
    function search(v)
      preorder[v] = (C[] += 1)
      push!(S, v)
      push!(P, v)
      for w in cfg[v]
        w in blocks || continue
        if preorder[w] == 0
          search(w)
        elseif !any(c -> w in c, components)
          while preorder[P[end]] > preorder[w] pop!(P) end
        end
      end
      if P[end] == v
        pop!(P)
        pushfirst!(components, Int[])
        while true
          push!(components[1], pop!(S))
          components[1][end] == v && break
        end
      end
    end
    for b in blocks
      preorder[b] == 0 && search(b)
    end
  end
  return components
end

struct Component
  children::Vector{Union{Component,Int}}
end

blocks(c::Integer) = [c]
blocks(c::Component) = reduce(vcat, blocks.(c.children))

entries(c::Integer) = [c]
entries(c::Component) = entries(first(c.children))

function components(cfg::CFG; blocks = 1:length(cfg))
  # Assume the first block is the entry.
  cs = strongconnected(cfg, blocks = blocks[2:end])
  Component([blocks[1],
             [length(c) == 1 ? c[1] : components(cfg, blocks = sort(c)) for c in cs]...])
end

# Aaaargh
_union(xs...) = union(xs...)
_union() = []

function branchesto(cfg::CFG, cs, i)
  c = cs.children[i]
  valid = _union([entries(cs.children[j]) for j = 1:length(cs.children) if i != j]...)
  preds = _union([cfg'[c] for c in blocks(c)]...)
  forw = filter(next -> next < entries(c)[1] && next in valid, preds)
  back = filter(next -> next > entries(c)[1] && next in valid, preds)
  forw = isempty(forw) ? nothing : minimum(forw)
  back = isempty(back) ? nothing : maximum(back)
  forw, back
end

function stackify(cfg::CFG, cs::Component = components(cfg), forw = [], back = [])
  entry = cs.children[1]
  for (i, c) = enumerate(cs.children)
    target = minimum(entries(c))
    f, b = branchesto(cfg, cs, i)
    f == nothing || push!(forw, entry=>target)
    b == nothing || push!(back, b    =>target)
    c isa Component && stackify(cfg, c, forw, back)
  end
  return forw, back
end


text:
# # Examples

## load packages
using PyPlot
using WaterPhysics
FT = Float32;

## define preview_data function
function preview_data(
            xs::Array{FT,1},
            ys::Array{FT,1},
            title::String,
            xlab::String,
            ylab::String
) where {FT<:AbstractFloat}
    fig = figure(title, figsize=(4,3));
    ax0 = fig.add_subplot(1,1,1);
    ax0.plot(xs, ys);
    ax0.set_xlabel(xlab);
    ax0.set_ylabel(ylab);
    fig.set_tight_layout(true);

    return fig
end
#------------------------------------------------------------------------------

# ## Capilary pressure
# ### As a function of radius
Rs = collect(FT, 1:0.5:30) .* FT(1e-6);
Pc = capillary_pressure.(Rs, FT(298.15));

preview_data(Rs .* FT(1e6), Pc .* FT(1e-3),
             "Pc ~ Radius", "Radius (μm)", "Capillary pressure (kPa)")
#------------------------------------------------------------------------------

# ### As a function of contact angle
αs = collect(FT, 0:1:90);
Pc = capillary_pressure.(FT(10 * 1e-6), FT(298.15), αs);

preview_data(αs, Pc .* FT(1e-3),
             "Pc ~ α", "Contact angle (°)", "Capillary pressure (kPa)")
#------------------------------------------------------------------------------

# ### As a function of temperature (due to surface tension)
Ts = collect(FT, 0:1:50) .+ FT(273.15);
Pc = capillary_pressure.(FT(10 * 1e-6), Ts);

preview_data(Ts .- FT(273.15), Pc .* FT(1e-3),
             "Pc ~ T", "Temperature (°C)", "Capillary pressure (kPa)")
#------------------------------------------------------------------------------

# ## Diffusive coefficient
# ### Relative coefficient as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
Dr = relative_diffusive_coefficient.(Ts);

preview_data(Ts .- FT(273.15), Dr,
             "rDc ~ T", "Temperature (°C)", "Relative diffusive coefficient")
#------------------------------------------------------------------------------

# ### Relative coefficient as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
Ds = diffusive_coefficient.(Ts, [TraceGasCO₂{FT}()], [TraceGasAir{FT}()]);

preview_data(Ts .- FT(273.15), Ds,
             "Dc ~ T", "Temperature (°C)", "Diffusive coefficient (m² s⁻¹)")
#------------------------------------------------------------------------------

# ## Latent heat
# ### λ of liquid water as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
λs = latent_heat_vapor.(Ts);

preview_data(Ts .- FT(273.15), λs,
             "λ ~ T", "Temperature (°C)", "Latent heat (J kg⁻¹)")
#------------------------------------------------------------------------------

# ## Saturation vapor pressure
# ### Saturated water vapor pressure at flat surface
Ts = collect(FT, 0:1:50) .+ FT(273.15);
Ps = saturation_vapor_pressure.(Ts);

preview_data(Ts .- FT(273.15), Ps,
             "Ps ~ T", "Temperature (°C)", "Saturation vapor pressure (Pa)")
#------------------------------------------------------------------------------

# ### Saturated water vapor pressure at curved surface (Kelvin equation)
Ψs = collect(FT, -5:0.1:5) .* FT(1e6);
Ps = saturation_vapor_pressure.(FT(298.15), Ψs);

preview_data(Ψs .* FT(1e-6), Ps,
             "Ps ~ Ψs",
             "Water potential (MPa)", "Saturation vapor pressure (Pa)")
#------------------------------------------------------------------------------

# ## Surface tension
# ### Relative γ of water-air interface as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
γr = relative_surface_tension.(Ts);

preview_data(Ts .- FT(273.15), γr,
             "rγ ~ T", "Temperature (°C)", "Relative surface tension")
#------------------------------------------------------------------------------

# ### γ of water-air interface as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
γs = surface_tension.(Ts);

preview_data(Ts .- FT(273.15), γs,
             "γ ~ T", "Temperature (°C)", "Surface tension (N m⁻¹)")
#------------------------------------------------------------------------------

# ## Viscosity
# ### Relative μ of water as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
μr = relative_viscosity.(Ts);

preview_data(Ts .- FT(273.15), μr,
             "rμ ~ T", "Temperature (°C)", "Relative viscosity")
#------------------------------------------------------------------------------

# ### Relative μ of water as a function of temperature
Ts = collect(FT, 0:1:50) .+ FT(273.15);
μs = viscosity.(Ts);

preview_data(Ts .- FT(273.15), μs,
             "μ ~ T", "Temperature (°C)", "Viscosity (Pa s)")
#------------------------------------------------------------------------------


text:
type BelegunduFitnessScheme{T} <: ParetoFitnessScheme{T}
  maximise::Vector{Bool}
end

"""
Implements a pareto-based multiple-objective fitness scheme through the use
of Belegundu fitness ranking; individuals belonging to the pareto front are
assigned fitness 0, whilst those outside are assigned 1.
"""
function belegundu(s::Dict{Any, Any})
   s["of"] = eval(Base.parse(Base.get(s, "of", "Float")))
   BelegunduFitnessScheme{s["of"]}(s["maximise"])
end

scale!{F}(s::BelegunduFitnessScheme, inds::Vector{Tuple{Int, F}}) =
  for (p1_id, p1_fitness) in inds
    fitness.rank = any((p2_id, p2_fitness) -> dominates(p2_fitness.scores, p1_fitness.scores, s.maximise), inds) ? 1 : 0
  end


text:
abstract type SSMHLoop <: Loop end
struct SSMHAlg <: SamplingAlgorithm end
"Single Site Metropolis Hastings"
const SSMH = SSMHAlg()

softhard(::Type{SSMHAlg}) = IsSoft{SSMHAlg}()

isapproximate(::SSMHAlg) = true

# defΩ(::SSMH) = SimpleΩ{Vector{Int}, Float64}

"Compute a score using the change in Prior of the *single* changed site"
function proposalkernel(kernel::Function, x)
  ∇logdensity(x) = x |> transform |> jacobian |> abs |> log
  before = ∇logdensity(x)
  proposed = kernel(x)
  after = ∇logdensity(proposed)
  ratio = after - before
  proposed, ratio
end

normalkernel(rng, x, σ = 0.1) = proposalkernel(x) do x
        inv_transform(transform(x) + σ * randn(rng))
      end
# normalkernel(rng, x::Array, σ = 0.1) = normalkernel.(x, σ)

"Metropolized Independent sample"
mi(rng, x::T) where T = (rand(rng, T), 0.0)

"Changes a uniformly chosen single site with kernel"
function swapsinglesite(transitionkernel::Function, rng, ω)
  logtranstionp = 0.0
  function updater(x)
    result, logtranstionp = transitionkernel(x)
    result
  end
  update(ω, rand(1:nelem(ω)), updater), logtranstionp
end

"Changes a uniformly chosen single site with kernel"
swapsinglesite(rng, ω) = swapsinglesite(rng, ω) do x
  mi(rng, x)
end

function moveproposal(rng, ω)
  swapsinglesite(rng, ω) do x
    normalkernel(rng, x, 1.0)
  end
end

"Normal propsoal with variance `var`"
function varproposal(var)
  function k(rng, ω)
    swapsinglesite(rng, ω) do x
      normalkernel(rng, x, var)
    end
  end
end

"""
Sample from `ω::Ω` conditioned on any constraints its conditioned on.

$(SIGNATURES)

# Arguments
- `x`: Real valued random variable
- `n`: Number of samples
- `logdensity`:  τ-valued `RandVar` s.t. logerr(τ) is defined
- `propsal`: function ω::Omega -> ω::Omega
- `ωinit`: Initial omega to start chain from
"""
function Base.rand(rng::AbstractRNG,
                   ΩT::Type{OT},
                   logdensity::RandVar,
                   n::Integer,
                   alg::SSMHAlg;
                   proposal = moveproposal,
                   ωinit = ΩT(),
                   offset = 0) where {OT <: Ω}
  ω = ωinit
  plast = logdensity(ω)
  qlast = 1.0
  ωsamples = OT[]
  accepted = 0
  for i = 1:n
    ω_, logtransitionp = isempty(ω) ? (ω,0) : proposal(rng, ω)
    p_ = logdensity(ω_)
    ratio = p_ - plast + logtransitionp
    if log(rand(rng)) < ratio
      ω = ω_
      plast = p_
      accepted += 1
    end
    push!(ωsamples, deepcopy(ω))
    lens(SSMHLoop, (ω = ω, ωsamples = ωsamples, accepted = accepted, p = plast, i = i + offset))
  end
  # println("accepted, ", accepted)
  ωsamples
end


text:
export SymInteger

struct SymInteger <: Integer
    ex::Expression
end

SymInteger(name::Symbol) = SymInteger(Variable(name))
SymInteger(x::SymInteger) = x
SymInteger(x::Integer) = SymInteger(Numeric(x))
SymInteger(x::Irrational{sym}) where sym = SymInteger(Constant(sym))

DOMAIN_TYPES[:Integer] = SymInteger

term(x::SymInteger) = x.ex
Base.show(io::IO, t::SymInteger) = print(io, t.ex)

for (M, f, arity) in DiffRules.diffrules()
    op = :($M.$f)

    if arity == 1
        @eval $op(x::SymInteger) = track($op, x)
    elseif arity == 2
        @eval $op(x::SymInteger, y::SymInteger) = track($op, x, y)
    end
end

Base.promote_rule(::Type{<:SymInteger}, ::Type{<:Integer}) = SymInteger
Base.convert(::Type{<:SymInteger}, x::SymInteger) = x
Base.convert(::Type{<:SymInteger}, x::Integer) = SymInteger(x)
Base.convert(::Type{<:SymInteger}, x::Irrational{sym}) where sym = SymInteger(x)

# Int specific
Base.xor(a::SymInteger, b::SymInteger) = track(xor, a, b)


text:
export GaussianStateBHD, QuantumStateBHD

struct QuantumStateBHD{T<:AbstractMatrix}<:ContinuousMultivariateDistribution
    ρ::T
end

"""
    GaussianStateBHD

A Gaussian state distribution in intensity-to-measurement-phase quadrature coordinate
measured by balanced homodyne detection.

## Example

```julia-repl
julia> using QuantumStateBase

julia> d = GaussianStateBHD(SqueezedState(0.8, π, Matrix, dim=100));

julia> points = rand(d, 4096);
```
"""
struct GaussianStateBHD{T<:AbstractMatrix}<:ContinuousMultivariateDistribution
    ρ::T
end

Base.length(d::GaussianStateBHD) = 2

sampler(d::GaussianStateBHD) = d

Base.eltype(::GaussianStateBHD{T}) where {T} = eltype(T)

function Distributions._rand!(rng::AbstractRNG, d::GaussianStateBHD, p::AbstractVector)
    T = real(eltype(d))

    p = real(p)
    p[1] = 2π * rand(rng, T)
    p[2] = mean(d, p[1]) + std(d, p[1]) * randn(rng, T)

    return p
end

function Distributions._rand!(rng::AbstractRNG, d::GaussianStateBHD, ps::AbstractMatrix)
    T = real(eltype(d))
    n = size(ps, 2)

    ps = real(ps)
    ps[1, :] .= sort!(2π .* rand(rng, T, n))
    ps[2, :] .= mean(d, view(ps, 1, :)) .+ std(d, view(ps, 1, :)) .* randn(rng, T, n)

    return ps
end

# Distributions._logpdf(d::GaussianStateBHD, x::AbstractArray)

mean(d::GaussianStateBHD, θs) = real(mean_of_π̂ₓ(d.ρ, θs))

var(d::GaussianStateBHD, θs) = real(mean_of_π̂ₓ²(d.ρ, θs) .- mean_of_π̂ₓ(d.ρ, θs).^2)

std(d::GaussianStateBHD, θs) = real(.√(mean_of_π̂ₓ²(d.ρ, θs) .- mean_of_π̂ₓ(d.ρ, θs).^2))

# entropy(d::GaussianStateBHD)

# cov(d::GaussianStateBHD)

# ##### for Gaussian state in intensity-to-measurement-phase quadrature coordinate #####

# π̂ₓ = (â exp(-im θ) + â† exp(im θ)) / 2

function mean_of_create(ρ::AbstractMatrix)
    # ⟨ψ|â†|ψ⟩ = tr(â† ρ)
    # Creation = [
    #     0   0   0   …
    #     √1  0   0   …
    #     0   √2  0   …
    #     0   0   √3  …
    #     ⋮   ⋮   ⋮
    # ]
    diag_1_of_ρ = diag(ρ, 1)
    dim = length(diag_1_of_ρ)
    T = eltype(ρ)

    return sum(.√(T.(1:dim)) .* diag_1_of_ρ)
end

function mean_of_create²(ρ::AbstractMatrix)
    # ⟨ψ|â†â†|ψ⟩ = tr((â†)² ρ)
    # Creation² = [
    #     0       0       0     …
    #     0       0       0     …
    #     √1√2    0       0     …
    #     0       √2√3    0     …
    #     0       0       √3√4  …
    #     ⋮       ⋮       ⋮
    # ]
    diag_2_of_ρ = diag(ρ, 2)
    dim = length(diag_2_of_ρ)
    T = eltype(ρ)

    return sum(.√(T.(1:dim)) .* .√(T.(2:(dim+1))) .* diag_2_of_ρ)
end

function mean_of_annihilate(ρ::AbstractMatrix)
    # ⟨ψ|â|ψ⟩ = tr(â ρ)
    # Annihilation = [
    #     0   √1  0   0   …
    #     0   0   √2  0   …
    #     0   0   0   √3  …
    #     ⋮   ⋮   ⋮   ⋮
    # ]
    diag_n1_of_ρ = diag(ρ, -1)
    dim = length(diag_n1_of_ρ)
    T = eltype(ρ)

    return sum(.√(T.(1:dim)) .* diag_n1_of_ρ)
end

function mean_of_annihilate²(ρ::AbstractMatrix)
    # ⟨ψ|ââ|ψ⟩ = tr((â)² ρ)
    # Annihilation² = [
    #    0        0       √1√2    0       0     …
    #    0        0       0       √2√3    0     …
    #    0        0       0       0       √3√4  …
    #    ⋮        ⋮        ⋮       ⋮       ⋮
    # ]
    diag_n2_of_ρ = diag(ρ, -2)
    dim = length(diag_n2_of_ρ)
    T = eltype(ρ)

    return sum(.√(T.(1:dim)) .* .√(T.(2:(dim+1))) .* diag_n2_of_ρ)
end

function mean_of_create_annihilate(ρ::AbstractMatrix)
    # ⟨ψ|â†â|ψ⟩ = tr(â†â ρ)
    # Creation * Annihilation = [
    #     0   0   0   0   …
    #     0   1   0   0   …
    #     0   0   2   0   …
    #     0   0   0   3   …
    #     ⋮   ⋮   ⋮   ⋮
    # ]
    diag_0_of_ρ = diag(ρ)
    dim = length(diag_0_of_ρ)
    T = eltype(ρ)

    return sum(T.(0:(dim-1)) .* diag_0_of_ρ)
end

function mean_of_π̂ₓ²(ρ::AbstractMatrix, θs)
    # ⟨πₓ²⟩ = ⟨ââ exp(-2im θ) + â†â† exp(2im θ) + ââ† + â†â⟩ / 4
    # ⟨πₓ²⟩ = (exp(-2im θ)⟨â²⟩ + exp(2im θ)⟨â†²⟩ + 1 + 2⟨ââ†⟩) / 4
    # here, ⟨ââ† + â†â⟩ = 1 + 2⟨â†â⟩ due to the commutation relation
    return (
        exp.(-2im*θs) .* mean_of_annihilate²(ρ) .+
        exp.(2im*θs) .* mean_of_create²(ρ) .+
        1 .+ 2mean_of_create_annihilate(ρ)
    ) ./ 4
end

function mean_of_π̂ₓ(ρ::AbstractMatrix, θs)
    # ⟨πₓ⟩ = ⟨â exp(-im θ) + â† exp(im θ)⟩ / 2
    # ⟨πₓ⟩ = (exp(-im θ)⟨â⟩ + exp(im θ)⟨â†⟩) / 2
    return (
        exp.(-im*θs) .* mean_of_annihilate(ρ) .+
        exp.(im*θs) .* mean_of_create(ρ)
    ) ./ 2
end


text:
function weightforces!(node::Int64, params::Parameters)
  params.f[3, 1, node] = -params.nodes[node].weight * sin(params.nodes[node].dip)
  params.f[1, 1, node] = -params.nodes[node].weight * cos(params.nodes[node].dip)
end


text:
using Documenter, NFFT

makedocs(
    format = :html,
    modules = [NFFT],
    sitename = "NFFT",
    authors = "Tobias Knopp,...",
    pages = [
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Directional" => "directional.md",
        "Density" => "density.md"
    ],
    html_prettyurls = false, #!("local" in ARGS),
)

deploydocs(repo   = "github.com/tknopp/NFFT.jl.git",
           target = "build")


text:
module Inventory

using LightGraphs, MetaGraphs

struct GraphDatabase
    g::MetaGraph
    GraphDatabase() = new(MetaGraph())
end

LightGraphs.add_vertex!(gd::GraphDatabase, attrs) = add_vertex!(gd.g, Dict(attrs))
LightGraphs.add_edge!(gd::GraphDatabase, i, j, attrs) = add_edge!(gd.g, i, j, Dict(attrs))
MetaGraphs.props(gd::GraphDatabase, args...) = MetaGraphs.props(gd.g, args...)

abstract type Item end
Base.broadcastable(x::Item) = Ref(x)

function add_item!(gd::GraphDatabase, item::T; tags=String[], description="", quantity=1) where {T <: Item}
    fnames = fieldnames(T)
    description = endswith(description, ".") ? description : description * "."
    values = getproperty.(item, fnames)
    add_vertex!(gd, Dict(
        (fnames .=> values)...,
        :description => description,
        :tags => tags,
        :quantity => quantity,
        :item => item
    ))
end

function add_item!(gd::GraphDatabase, item_desc::String; tags=String[])
    add_vertex!(gd, Dict(desc_string(item_desc)..., :tags => tags))
end

abstract type Relationship end
label(::Type{T}) where {T <: Relationship} = lowercase(string(Symbol(T)))

find_vert(gd::GraphDatabase, i::Number) = i
function find_vert(gd::GraphDatabase, item::Item)
    verts = collect(filter_vertices(gd.g, :item, item))
    lv = length(verts)
    @assert lv <= 1 "More than one vertex found for $item (found $lv)"
    @assert lv > 0 "No vertex found for $item"
    first(verts)
end

add_relationship!(gd::GraphDatabase, relationship::T) where {T <: Relationship} = add_edge!(gd, find_vert(gd, relationship.src), find_vert(gd, relationship.dst), Dict(:relationship => label(T)))

find_vertices(gd::GraphDatabase, properties::Dict{Symbol, <: Any}) = collect(filter_vertices(gd.g, (g, x) -> all(props(g, x)[prop] == val for (prop, val) ∈ properties)))
find_vertices(gd::GraphDatabase, properties::Vector{Symbol}) = collect(filter_vertices(gd.g, (g, x) -> all(haskey(props(g, x), prop) for prop ∈ properties)))
find_edges(gd::GraphDatabase, relationship::Type{T}) where {T <: Relationship} = collect(filter_edges(gd.g, :relationship, label(T)))

Base.show(io::IO, gd::GraphDatabase) = print(io, "Graph database ($(length(gd.g)) items, $(length(edges(gd.g))) relationships)")

function desc_string(str)
    quantity_m = match(r"x\d+(?=(?:$|\s-))", str)
    quantity = isnothing(quantity_m) ? 1 : parse(Int, strip(quantity_m.match, 'x'))
    desc_m = match(r"\s-\s(.*)$", str)
    desc = isnothing(desc_m) ? nothing : first(desc_m.captures)
    if isnothing(quantity_m) && isnothing(desc_m)
        name_reg = r"^(.*)$"
    elseif isnothing(quantity_m)
        name_reg = r"^(.*)\s-\s"
    else
        name_reg = r"^(.*)\sx\d+"
    end
    name = match(name_reg, str).captures[1]
    Dict(:name => name, :description => desc, :quantity => quantity)
end

export GraphDatabase,
       Relationship,
       Item,
       find_vert,
       label,
       add_relationship!,
       add_item!,
       props,
       find_vertices,
       find_edges
end


text:
module MatrixFactorisation

using DataModel
using MFUtilities
using c_wrappers
using MLUtilities
using Distributions

export MatrixFactorisationModel

"""
This model implements a Bayesian matrix factorisation model for collaborative filtering following

Ahn et al. Large-Scale Distributed Bayesian Matrix Factorization using Stochastic Gradient MCMC

Ratings are modelled as Rᵢⱼ- mean_rating ∼ N(dot(Uᵢ,Vⱼ)+aᵢ+bⱼ,τ) where U,V,a,b are assumed to be normally distributed with mean 0 and precisions λᵤ (one per dim),λᵥ (one per dim),λₐ,λᵦ with Γ(α₀,β₀) hyperpriors.

Fields:
    - block: data in the format Nx3; first column user id between 1 and nusers; second column item id between 1 and nitems; third column rating.
    - nusers
    - nitems
    - d: dimension of the user and item feature vectors
    - τ: noise variance
    - α₀,β₀:  hyperprior parameters
    - n_elem_row: number of occurences per user
    - n_elem_col: number of occurences per item
    - inv_pr_pick_users, inv_pr_pick_items: needed for specialized sampler from Ahn et al.
    - batchiterator, batchstate: provides minibatches
    - ntrain: size of dataset
    - batchsize
    - mean_rating
"""
type MatrixFactorisationModel <: AbstractDataModel
    block::Array{Float64}
    nusers::Int64
    nitems::Int64
    d::Int64
    τ::Float64
    α₀::Float64
    β₀::Float64
    n_elem_row
    n_elem_col
    inv_pr_pick_users::Array{Float64}
    inv_pr_pick_items::Array{Float64}
    batchiterator::BatchIterator
    batchstate
    ntrain::Float64 # number of ratings in this block
    batchsize::Int64
    mean_rating::Float64
    function MatrixFactorisationModel(block,d,batchsize;τ= 2,α₀=1.,β₀=1., randomize=true)
        # get user and items numbers
        nusers    =  round(Int64,maximum(block[:,1]))
        nitems    = round(Int64,maximum(block[:,2]))

        # get data size
        ntrain = size(block)[1]

        # get number of user and item occurrence
        n_elem_row, n_elem_col = get_n_elem(block, nusers,nitems)

        # get auxiliary information about
        inv_pr_pick_users, inv_pr_pick_items = comp_inv_pr_pick(n_elem_row,n_elem_col,batchsize)

        # set up batch iterator
        batchiterator = BatchIterator(ntrain,batchsize,Inf,randomize=randomize)
        batchstate = start(batchiterator)

        mean_rating = mean(block[:,3])

        if randomize
            block = block[randperm(ntrain),:]
        end

        new(block,nusers,nitems,d,τ,α₀,β₀,n_elem_row,n_elem_col,inv_pr_pick_users,inv_pr_pick_items,batchiterator,batchstate,ntrain,batchsize,mean_rating)
    end
end

"""
fetchnparams returns the total number of parameters in the model.
"""
fetchnparams(dm::MatrixFactorisationModel) = (dm.nusers+dm.nitems)*(dm.d+1)

"""
getbatch returns the next minibatch
"""
function getbatch(dm::MatrixFactorisationModel)
    idx, dm.batchstate = next(dm.batchiterator,dm.batchstate)
    batch = dm.block[idx,:]
    return idx, batch
end

"""
getparams returns views of the user feature matrix U, the user biases a, the item feature matrix V and the item biases b
"""
function getparams(x::Array{Float64},dm::MatrixFactorisationModel)
    # D dx1
    # U n_u x d
    # a n_u x 1
    # V n_i x d
    # b n_i x 1

    d = dm.d
    nusers = dm.nusers
    nitems = dm.nitems

    ustart = 1
    uend = d*nusers
    astart = uend + 1
    aend = uend + nusers
    vstart = aend + 1
    vend = aend + d*nitems
    bstart = vend + 1
    bend = vend + nitems
    U = reshape(view(x,ustart:uend),(d,nusers))
    a = view(x,astart:aend)
    V = reshape(view(x,vstart:vend),(d,nitems))
    b = view(x,bstart:bend)
    return U,a,V,b
end

"""
getparams_λ returns views of the
"""
function getparams_λ(λ::Array{Float64},dm::MatrixFactorisationModel)
    d = dm.d

    @assert length(λ) == 2d+2
    λUstart = 1
    λUend = d
    λVstart = d+ 1
    λVend = 2d
    λaid = 2d+1
    λbid = 2d+2
    λU = view(λ,λUstart:λUend)
    λa = view(λ,λaid)
    λV = view(λ,λVstart:λVend)
    λb = view(λ,λbid)
    return λU,λV,λa,λb
end

"""
getgrad returns a function that gives the next stochastic gradient
"""
function DataModel.getgrad(dm::MatrixFactorisationModel;λ::Array{Float64,1}=ones(2*dm.d+2))
    function mf_gradient(x)
        err = Array(Float64,dm.batchsize);
        U,a,V,b = getparams(x,dm);
        idx, batch = getbatch(dm)
        rr = batch[:,3];
        uu = batch[:,1];
        ii = batch[:,2];

        ux, ix     = unique(uu), unique(ii)
        c_wrappers.c_comp_error!(err,rr,dm.mean_rating,U,V,a,b,uu,ii,dm.batchsize,dm.d);

        grad = zeros(fetchnparams(dm))
        gradU,grada,gradV,gradb = getparams(grad,dm)
        c_wrappers.c_comp_grad_sum!(err, uu, ii, U, V, a, b, gradU, gradV, grada, gradb, dm.batchsize, dm.d)

        return grad
    end
    return x -> dm.ntrain/dm.batchsize*mf_gradient(x) + gradprior(dm,x,λ)
end

"""
getllik returns a function that gives the next stochastic loglikehood estimate
"""
function DataModel.getllik(dm::MatrixFactorisationModel;λ::Array{Float64,1}=ones(2*dm.d+2))
    function loglikehood(x)
        err = Array(Float64,dm.batchsize);
        U,a,V,b = getparams(x,dm)
        λᵤ,λᵥ,λₐ,λᵦ = getparams_λ(λ,dm)
        idx, batch = getbatch(dm)
        rr = batch[:,3];
        uu = batch[:,1];
        ii = batch[:,2];

        c_wrappers.c_comp_error!(err,rr,dm.mean_rating,U,V,a,b,uu,ii,dm.batchsize,dm.d);
        logprior = -(fetchnparams(dm))/2*log(2*pi)+0.5*sum(log(λ))-0.5*(sum(λᵤ.*sum(U.*U,2))) - 0.5*(sum(λᵥ.*sum(V.*V,2))) - 0.5*(sum(λₐ*sum(a.*a))) - 0.5*(sum(λᵦ.*sum(b.*b)))
        llik = logprior - dm.ntrain/dm.batchsize*sum(0.5*err.^2) - dm.ntrain/2*log(2*π)

        return llik[1]
    end
    return loglikehood
end

"""
auxiliary function to calculate the prior gradient
"""
function gradprior(dm::MatrixFactorisationModel,x::Array{Float64}, λ::Array{Float64})

    λᵤ,λᵥ,λₐ,λᵦ = getparams_λ(λ,dm)

    U,a,V,b = getparams(x,dm)
    ∇ = zeros(fetchnparams(dm))
    ∇ᵤ,∇ₐ,∇ᵥ,∇ᵦ = getparams(∇,dm)
    ∇ᵤ[:] = (-λᵤ.*U)[:]
    ∇ᵥ[:] = (-λᵥ.*V)[:]
    ∇ₐ[:] = -λₐ.*(a)
    ∇ᵦ[:] = -λᵦ.*(b)
    return ∇
end

"""
evaluate returns the RMSE given parameters `x`
"""
function evaluate(dm::MatrixFactorisationModel,
  x::Vector{Float64}
  )
  prediction = predict(dm,x)
  evaluateprediction(dm,prediction)
end

"""
evaluate returns the RMSE given parameters `x` for a test model. Since some users/items might not be present in the test set, nusers/nitems has to be supplied to allow the parameters to be matched correctly.
"""
function evaluate_test(dm::MatrixFactorisationModel,
  x::Vector{Float64}, dmtrain::MatrixFactorisationModel)
  prediction = predict_test(dm,x,dmtrain)
  evaluateprediction(dm,prediction)
end

"""
evaluateprediction return the RMSE given predictions
"""
function evaluateprediction(dm::MatrixFactorisationModel,
  prediction::Vector{Float64})
  diff = dm.block[:,3] - prediction
  sse = sum(diff.*diff)
  rmse = sqrt(sse./dm.ntrain)
  Dict(:rmse => rmse)
end

"""
predict gives predicted ratings given parameters `x`
"""
function predict(dm::MatrixFactorisationModel,x::Vector{Float64})
    U,a,V,b = getparams(x,dm)
    prediction = zeros(round(Int64,dm.ntrain))
    c_wrappers.c_predict!(dm.block,U,V,a,b,prediction,1.,5.,dm.mean_rating)
    prediction
end

"""
predict gives predicted ratings given parameters `x`
"""
function predict_test(dm::MatrixFactorisationModel,x::Vector{Float64},dmtrain::MatrixFactorisationModel)
    U,a,V,b = getparams(x,dmtrain)
    prediction = zeros(round(Int64,dm.ntrain))
    c_wrappers.c_predict!(dm.block,U,V,a,b,prediction,1.,5.,dm.mean_rating)
    prediction
end

"""
lambda_sample performs a Gibbs update for the precisions λ
"""
function lambda_sample(dm::MatrixFactorisationModel,x::Array{Float64,1})
    U,a,V,b = getparams(x,dm)
    d = dm.d
    nusers  = dm.nusers
    nitems = dm.nitems
    α₀ = dm.α₀
    β₀ = dm.β₀

    αᵤ = α₀ + dm.nusers/2
    αₐ = αᵤ
    αᵥ = α₀ + dm.nitems/2
    αᵦ  = αᵥ

    βᵤ = β₀ + 0.5.*sum(U.*U,2)[:]
    βᵥ = β₀ + 0.5.*sum(V.*V,2)[:]
    βₐ = β₀ + 0.5.*sum(a.*a)[1]
    βᵦ = β₀ + 0.5.*sum(b.*b)[1]

    λᵤ = sample_gamma(ones(d)*αᵤ,βᵤ)
    λᵥ = sample_gamma(ones(d)*αᵥ,βᵥ)
    λₐ = sample_gamma(αₐ,βₐ)
    λᵦ = sample_gamma(αᵦ,βᵦ)
    return [λᵤ;λₐ;λᵥ;λᵦ]
end

"""
auxiliary function to sample gamma variables
"""
function sample_gamma(α::Real,β::Real)
    try
        @assert α > 0.0 && β > 0.0
    catch
        @show α,β
    end
    d = Gamma(α ,1/β)
    return Base.rand(d)
end

"""
auxiliary function to sample gamma variables
"""
function sample_gamma(α::Array{Float64},β::Array{Float64})
    @assert size(α) == size(β)
    return Float64[sample_gamma(a,b) for (a,b) in zip(α,β)]
end

"""
Specialised sparse sampler from Ahn et al.
"""
type SparseSGLDState
    x::Array{Float64,1}
    grad_sum::Array{Float64,1}
    λ::Array{Float64,1}
    t::Int
    stepsize::Function
    niters::Int
    function SparseSGLDState(x::Array{Float64,1},lambda::Array{Float64,1},stepsize::Function; niters::Int = 10)
        grad_sum = zeros(length(x))
        new(x,grad_sum,lambda,0,stepsize,niters)
    end
    function SparseSGLDState(x::Array{Float64,1},lambda::Array{Float64,1},stepsize::Float64; niters::Int= 10)
        grad_sum = zeros(length(x))
        f = niters::Int -> stepsize
        new(x,grad_sum,lambda,0,f,niters)
    end
end

"""
one update using the sparse SGLD sampler.
"""
function sample_sparse!(dm::MatrixFactorisationModel,
   state::SparseSGLDState)
        factor = dm.ntrain/(dm.batchsize)
        state.λ = lambda_sample(dm,state.x)
        LambdaU,LambdaV,lambdaa, lambdab = getparams_λ(state.λ,dm)
        U,a,V,b = getparams(state.x,dm);
        grad_sum_U,grad_sum_a,grad_sum_V, grad_sum_b = getparams(state.grad_sum,dm);
        uxx = []
        ixx = []
        for k = 1:state.niters

            state.t += 1
            eps = state.stepsize(state.t)
            err = Array(Float64,dm.batchsize);

            idx, batch = getbatch(dm)
            rr = batch[:,3];
            uu = batch[:,1];
            ii = batch[:,2];

            ux, ix     = unique(uu), unique(ii)
            if k == 1
                uxx = ux
                ixx = ix
            else
                uxx = vcat(uxx,ux)
                ixx = vcat(ixx,ix)
            end
            c_wrappers.c_comp_error!(err,rr,dm.mean_rating,U,V,a,b,uu,ii,dm.batchsize,dm.d);
            halfstep = 0.5 * eps
            sqrtstep = sqrt(eps)

            c_wrappers.c_comp_grad_sum!(err, uu, ii, U, V, a, b, grad_sum_U, grad_sum_V, grad_sum_a, grad_sum_b, dm.batchsize, dm.d)

            c_wrappers.c_update_para_prior!(ux, U, grad_sum_U, LambdaU, zeros(U), zeros(U), dm.inv_pr_pick_users, factor, halfstep, sqrtstep, dm.d)
            c_wrappers.c_update_para_prior!(ix, V, grad_sum_V, LambdaV, zeros(V), zeros(V), dm.inv_pr_pick_items, factor, halfstep, sqrtstep, dm.d)
            c_wrappers.c_update_para_prior!(ux, a, grad_sum_a, lambdaa, zeros(a), zeros(a), dm.inv_pr_pick_users, factor, halfstep, sqrtstep, 1)
            c_wrappers.c_update_para_prior!(ix, b, grad_sum_b, lambdab, zeros(b), zeros(b), dm.inv_pr_pick_items, factor, halfstep, sqrtstep, 1)
        end
        uxx = unique(uxx)
        ixx = unique(ixx)
        uxx, ixx
end

end


text:
using HPAT

@acc hpat function q25(d_date, file_name)
    store_sales = DataSource(DataTable{:ss_customer_sk=Int64, :ss_ticket_number=Int64, :ss_sold_date_sk=Int64, :ss_net_paid=Float64}, HDF5, file_name)
    web_sales = DataSource(DataTable{:ws_bill_customer_sk=Int64, :ws_order_number=Int64, :ws_sold_date_sk=Int64, :ws_net_paid=Float64}, HDF5, file_name)

    store_sales = store_sales[:ss_sold_date_sk > d_date]
    web_sales   = web_sales[:ws_sold_date_sk > d_date]

    store_agg = aggregate(store_sales, :ss_customer_sk,
                                           :frequency = count(union(:ss_ticket_number)),
                                           :most_recent_date = max(:ss_sold_date_sk),
                                           :amount = sum(:ss_net_paid))[:cid = :ss_customer_sk]
    web_agg = aggregate(web_sales,  :ws_bill_customer_sk,
                                           :frequency = count(union(:ws_order_number)),
                                           :most_recent_date = max(:ws_sold_date_sk),
                                           :amount = sum(:ws_net_paid))[:cid = :ws_bill_customer_sk]
    web_store_agg = [store_agg; web_agg]

    result = aggregate(web_store_agg, :cid,
                                      :recency = (37621 - max(:most_recent_date) < 60 ? 1.0 : 0.0),
                                      :frequency = sum(:frequency),
                                      :totalspend = sum(:amount))

    sort!(result, by=:cid)
end

println(q25("1-1-2015", "data.hdf5"))


text:
"""
    Base.show(io::IO, connection::CKANConnection)

Print a CKANConnection to the given IO.
"""
function Base.show(io::IO, connection::CKANConnection)
    obj_print = """

                ✔ [CKAN API connection]
                ↳ [URL: $(connection.url)] """
    print(io, obj_print)
end

"""
    Base.show(io::IO, dataset::CKANDataset)

Print a CKANDataset to the given IO.
"""
function Base.show(io::IO, dataset::CKANDataset)
    obj_print = """

                Δ [CKAN Dataset]
                ↳ [id    : $(dataset.id)]
                ↳ [title : $(dataset.metadata["title"])] """
    print(io, obj_print)
end

"""
    Base.show(io::IO, resource::CKANResource)

Print a CKANResource to the given IO.
"""
function Base.show(io::IO, resource::CKANResource)
    obj_print = """

                ρ [CKAN Resource]
                ↳ [id   : $(resource.id)]
                ↳ [type : $(resource.metadata["resource_type"])]
                ↳ [name : $(resource.metadata["name"])] """
    print(io, obj_print)
end


text:
# https://discourse.julialang.org/t/ann-urldownload-jl/37841 #

import Pkg
Pkg.instantiate()

Pkg.add("UrlDownload")
Pkg.add("CSV")
Pkg.add("DataFrames")

using UrlDownload, DataFrames, CSV

url = "https://www.stats.govt.nz/assets/Uploads/Business-price-indexes/Business-price-indexes-December-2019-quarter/Download-data/business-price-indexes-december-2019-quarter-csv.csv"
# data_frame = urldownload(url)
data_frame = urldownload(url, true)

res = urldownload(url, parser = x -> DataFrame(CSV.File(IOBuffer(x))))

res
# url = "https://covidtracking.com/api/v1/states/current.csv"
# res = urldownload(url, parser = x -> DataFrame(CSV.File(IOBuffer(x))))


text:
export RadialMap, evaluate!, evaluate, SparseRadialMap

import Base: size, show

#### Structure for the lower triangular map M

struct RadialMap
        Nx::Int64
        p::Int64
        L::LinearTransform
        C::Array{RadialMapComponent, 1}

        # Optimization parameters
        γ::Float64
        λ::Float64
        δ::Float64
        κ::Float64
end

function RadialMap(Nx::Int64, p::Int64; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        C = RadialMapComponent[]
        @inbounds for i=1:Nx
                push!(C, RadialMapComponent(i,p))
        end
        return RadialMap(Nx, p, LinearTransform(Nx), C, γ, λ, δ, κ)
end

function RadialMap(X::Array{Float64,2}, p::Int64; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        L = LinearTransform(X; diag = true)
        Nx, Ne = size(X)

        C = RadialMapComponent[]
        @inbounds for i=1:Nx
                push!(C, RadialMapComponent(i,p))
        end
        return RadialMap(Nx, p, LinearTransform(Nx), C, γ, λ, δ, κ)
end

function Base.show(io::IO, M::RadialMap)
println(io,"Radial Map of dimension Nx = $(M.Nx) and order p = $(M.p)
with parameters (γ, λ, δ, κ) = ($(M.γ), $(M.λ), $(M.δ), $(M.κ))")

end

@propagate_inbounds Base.getindex(M::RadialMap, i::Int) = getindex(M.C,i)
@propagate_inbounds Base.setindex!(M::RadialMap, C::RadialMapComponent, i::Int) = setindex!(M.C,C,i)

size(M::RadialMap) = (M.Nx, M.p)

# Evaluate the map RadialMapComponent at z = (z1,...,zNx)
function evaluate!(out, M::RadialMap, z::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1)
        Nx = M.Nx
        @assert Nx == size(z,1) "Incorrect length of the input vector"

        # We can apply the rescaling to all the components once
        if apply_rescaling == true
                transform!(M.L, z)
        end

        @inbounds  for i=start:Nx
                out[i] = M.C[i](view(z,1:i))
        end

        if apply_rescaling == true
                itransform!(M.L, z)
        end

        return out
end

evaluate(M::RadialMap, z::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1) = evaluate!(zeros(M.Nx-start+1), M, z; apply_rescaling = apply_rescaling, start = start)

(M::RadialMap)(z::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1) =  evaluate(M, z; apply_rescaling = apply_rescaling, start = start)

# Evaluate in-place the RadialMap `M`
function evaluate!(out, M::RadialMap, X::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1)
        @get M (Nx, p)
        NxX, Ne = size(X)
        @assert NxX == Nx "Wrong dimension of the ensemble matrix `X`"

        if apply_rescaling == true
                transform!(M.L, X)
        end

        @inbounds Threads.@threads for i=1:Ne
                col = view(X,:,i)
                evaluate!(view(out,:,i), M, col; apply_rescaling = false, start =start)
                # out[:,i] .= M(col, start=start)
        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end
        return out
end

evaluate(M::RadialMap, X::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1) = evaluate!(zero(X), M, X; apply_rescaling = apply_rescaling, start = start)

(M::RadialMap)(X::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1) = evaluate(M, X; apply_rescaling = apply_rescaling, start = start)

#### Sparse Structure for the lower triangular map M

struct SparseRadialMap
        Nx::Int64
        p::Array{Array{Int64,1}}
        L::LinearTransform
        C::Array{SparseRadialMapComponent, 1}

        # Optimization parameters
        γ::Float64
        λ::Float64
        δ::Float64
        κ::Float64
end

function SparseRadialMap(Nx::Int64, p::Array{Array{Int64,1}}; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        @assert Nx==size(p,1) "Wrong dimension of the SparseRadialMap"
        L = LinearTransform(Nx)
        C = SparseRadialMapComponent[]
        @inbounds  for i=1:Nx
                push!(C, SparseRadialMapComponent(i,p[i]))
        end
        return SparseRadialMap(Nx, p, L, C, γ, λ, δ, κ)
end

function SparseRadialMap(Nx::Int64, p::Array{Int64,1}; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        @assert Nx==size(p,1) "Wrong dimension of the SparseRadialMap"
        L = LinearTransform(Nx)
        C = SparseRadialMapComponent[]
        @inbounds for i=1:Nx
                push!(C, SparseRadialMapComponent(i,p[i]))
        end
        return SparseRadialMap(Nx, [fill(p[i],i) for i=1:Nx], L, C, γ, λ, δ, κ)
end

function SparseRadialMap(Nx::Int64, p::Int64; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        L = LinearTransform(Nx)
        C = SparseRadialMapComponent[]
        @inbounds for i=1:Nx
                push!(C, SparseRadialMapComponent(i,p))
        end
        return SparseRadialMap(Nx, [fill(p,i) for i=1:Nx], L, C, γ, λ, δ, κ)
end

# Methods based on an ensemble matrix X
function SparseRadialMap(X::Array{Float64,2}, p::Array{Array{Int64,1}}; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        Nx, Ne = size(X)
        @assert Nx==size(p,1) "Wrong dimension of the SparseRadialMap"
        L = LinearTransform(X)
        C = SparseRadialMapComponent[]
        @inbounds  for i=1:Nx
                push!(C, SparseRadialMapComponent(i,p[i]))
        end
        return SparseRadialMap(Nx, p, L, C, γ, λ, δ, κ)
end

function SparseRadialMap(X::Array{Float64,2}, p::Array{Int64,1}; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        Nx, Ne = size(X)
        @assert Nx==size(p,1) "Wrong dimension of the SparseRadialMap"
        L = LinearTransform(X)
        C = SparseRadialMapComponent[]
        @inbounds for i=1:Nx
                push!(C, SparseRadialMapComponent(i,p[i]))
        end
        return SparseRadialMap(Nx, [fill(p[i],i) for i=1:Nx], L, C, γ, λ, δ, κ)
end

function SparseRadialMap(X::Array{Float64,2}, p::Int64; γ::Float64=2.0, λ::Float64=0.1, δ::Float64=1e-8, κ::Float64=10.0)
        Nx, Ne = size(X)
        L = LinearTransform(X)
        C = SparseRadialMapComponent[]
        @inbounds for i=1:Nx
                push!(C, SparseRadialMapComponent(i,p))
        end
        return SparseRadialMap(Nx, [fill(p,i) for i=1:Nx], L, C, γ, λ, δ, κ)
end

function Base.show(io::IO, M::SparseRadialMap)
        println(io,"Sparse Radial Map of dimension Nx = $(M.Nx) and order p = $(M.p)
        with parameters (γ, λ, δ, κ) = ($(M.γ), $(M.λ), $(M.δ), $(M.κ))")
end

@propagate_inbounds Base.getindex(M::SparseRadialMap, i::Int) = getindex(M.C,i)
@propagate_inbounds Base.setindex!(M::SparseRadialMap, C::SparseRadialMapComponent, i::Int) = setindex!(M.C,C,i)

size(M::SparseRadialMap) = (M.Nx, M.p)

# Evaluate the map Sparse RadialMap at z = (z1,...,zNx)

function evaluate!(out, M::SparseRadialMap, z::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1)
        Nx = M.Nx
        @assert Nx==size(z,1) "Incorrect length of the input vector"

        if apply_rescaling == true
                transform!(M.L, z)
        end

        @inbounds for i=start:Nx
                if allequal(M.C[i].p,-1)
                        out[i] = z[i]
                else
                        out[i] = M.C[i](view(z,1:i))
                end
        end

        if apply_rescaling == true
                itransform!(M.L, z)
        end

        return out
end

evaluate(M::SparseRadialMap, z::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1) = evaluate!(zeros(M.Nx-start+1), M, z; apply_rescaling = apply_rescaling, start = start)

(M::SparseRadialMap)(z::AbstractVector{Float64}; apply_rescaling::Bool=true, start::Int64=1) =  evaluate(M, z; apply_rescaling = apply_rescaling, start = start)

# Evaluate the map SparseRadialMapComponent at z = (z1,...,zNx)
function evaluate!(out, M::SparseRadialMap, X::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1)
        @get M (Nx, p)
        NxX, Ne = size(X)
        @assert NxX == Nx "Wrong dimension of the ensemble matrix `X`"
        @assert size(out) == (Nx, Ne) "Wrong dimension of the output matrix `out`"

        if apply_rescaling == true
                transform!(M.L, X)
        end

        @inbounds for i=1:Ne
                col = view(X,:,i)
                evaluate!(view(out,:,i), M, col; apply_rescaling = false, start = start)
        end

        if apply_rescaling == true
                itransform!(M.L, X)
        end

        return out
end

evaluate(M::SparseRadialMap, X::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1) = evaluate!(zero(X), M, X; apply_rescaling = apply_rescaling, start = start)

(M::SparseRadialMap)(X::AbstractMatrix{Float64}; apply_rescaling::Bool=true, start::Int64=1) =  evaluate(M, X; apply_rescaling = apply_rescaling, start = start)


text:
#start script
include("findnearest.jl")
include("findneighbours.jl")
include("removeedge.jl")
include("removenode_old.jl")
include("gwr.jl")
include("plotgwr.jl")

using MAT
using PyPlot

#using CPUtime
file = matopen("../local_uniform_2d.mat")
Data = read(file, "Data")
# scatter(Data[1,:]',Data[2,:]', )

tic()
A,(),(),() = gwr(Data,200);
toc()

scatter(A[1,:]', A[2,:]')


text:
# # Modal analysis of acoustic medium in a rigid box

# Source code: [`rigid_box_tut.jl`](rigid_box_tut.jl)

# ## Description

# Example from Boundary element acoustics: Fundamentals and computer codes, TW
# Wu, page 123. Internal resonance problem. Reference frequencies: 90.7895,
# 181.579, 215.625, 233.959, 272.368, 281.895 Quadrilateral mesh.

# ## Goals

# - Construct the discrete model and solve the eigenvalue problem for the
#   natural frequencies.
# - Visualize eigenmodes.

##
# ## Definitions

# The finite element code realize on the basic functionality implemented in this
# package.
using FinEtools
# The acoustics functionality is brought in:
using FinEtoolsAcoustics

# In order to solve the eigenvalue problem we need the Arnoldi package.
import Arpack: eigs

# The input quantities are provided without units, but the units are consistent:
rho = 1.21*1e-9;# mass density
c  = 345.0*1000;# millimeters per second
bulk =  c^2*rho;
Lx = 1900.0;# length of the box, millimeters
Ly = 800.0; # length of the box, millimeters
n = 14;#
neigvs = 18;
OmegaShift = 10.0;

# The mesh covers a rectangle.
fens,fes = Q4block(Lx,Ly,n,n); # Mesh

# Construct the geometry and the pressure field.
geom = NodalField(fens.xyz)
P = NodalField(zeros(size(fens.xyz,1),1))
# There are no boundary conditions explicitly needed: The box has hard sound
# boundary conditions all around. Number the unknowns in the pressure field.
numberdofs!(P)

# The finite element machine for acoustics is constructed. The integration
# rule is appropriate to the four node quadrilaterals.
femm = FEMMAcoust(IntegDomain(fes, GaussRule(2, 2)), MatAcoustFluid(bulk,rho))

# Compute the stiffness and mass matrices:
S = acousticstiffness(femm, geom, P);
C = acousticmass(femm, geom, P);

# And now solves eigenvalue problem. Note that the acoustic mass is singular,
# and we need to use shifting to ensure convertibility of the first matrix.
evals, evecs, nconv = eigs(C+OmegaShift*S, S; nev=neigvs, which=:SM)
evals = evals .- OmegaShift;
fs = real(sqrt.(complex(evals)))./(2*pi)
println("Eigenvalues: $fs [Hz]")

# Postprocessing: export all pressure modes.
File =  "rigid_box.vtk"
scalarllist = Any[]
for n  = 1:15
    push!(scalarllist, ("Pressure_mode_$n", evecs[:,n]));
end
vtkexportmesh(File, connasarray(fes), geom.values, FinEtools.MeshExportModule.VTK.Q4;
scalars=scalarllist)
@async run(`"paraview.exe" $File`)

true


text:
#=
VT-flash algorithm
=#

struct VTFlashResult{T}
    converged::Bool
    singlephase::Bool
    RT::T
    nmol_1::Vector{T}
    V_1::T
    nmol_2::Vector{T}
    V_2::T
    iterations::Int
    fcalls::Int

    function VTFlashResult{T}(
        converged,
        singlephase,
        RT,
        nmol_1,
        V_1,
        nmol_2,
        V_2,
        iterations=-1,
        fcalls=-1,
    ) where {T}
        return new{T}(
            converged,
            singlephase,
            RT,
            copy(nmol_1),
            V_1,
            copy(nmol_2),
            V_2,
            iterations,
            fcalls,
        )
    end
end

VTFlashResult{T}(; converged, singlephase, RT, nmol_1, V_1, nmol_2, V_2, iters=-1, fcalls=-1) where {T} =
VTFlashResult{T}(converged, singlephase, RT, nmol_1, V_1, nmol_2, V_2, iters, fcalls)

"Return concentration of state with minimum energy from vt-stability tries."
function __vt_flash_init_conc_choose(
    vt_stab_tries::AbstractVector{VTStabilityResult{T}},
) where {T}
    Dmin = T(Inf)
    index_min = -1
    for (i, state) in enumerate(vt_stab_tries)
        if !state.isstable && state.energy_density < Dmin
            index_min = i
            Dmin = state.energy_density
        end
    end
    return vt_stab_tries[index_min].concentration
end

"""
Calculates pressure gradient for Brusilovsky EoS at point (N₁,..., Nₙ, V).
∇P = [∂P/∂Nᵢ..., ∂P/∂V], where i = 1,...,`ncomponents(mix)`.
"""
function __vt_flash_pressure_gradient!(
    ∇P::AbstractVector{T},
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real;
    buf::BrusilovskyThermoBuffer=thermo_buffer(mix),
) where {T}
    # I did not implement this function in src/basic_thermo.jl
    # because the gradient here does not include ∂P/∂T derivative.
    # Maybe, it should be implemented later in src/basic_thermo.jl.

    A, B, C, D, aij = eos_parameters(mix, nmol, RT; buf=buf)

    # hell arithmetics
    # does compiler smart to detect this as constants
    # if plain operation were put in ∂P/∂Nᵢ for-cycle explicitly?
    V = volume  # alias
    VmB⁻¹ = 1 / (V - B)
    ΣnmolbyVmB² = sum(nmol) * VmB⁻¹^2
    DmC = D - C
    VpC⁻¹ = 1 / (V + C)
    VpC⁻² = VpC⁻¹^2
    VpD⁻¹ = 1 / (V + D)
    VpD⁻² = VpD⁻¹^2
    AbyDmC = A / DmC
    VpC⁻¹mVpD⁻¹byDmC² = (VpC⁻¹ - VpD⁻¹) / DmC^2

    # ∂P/∂Nᵢ part
    for (i, substance) in enumerate(components(mix))
        bᵢ, cᵢ, dᵢ = substance.b, substance.c, substance.d
        ∂ᵢA = 2 * dot(nmol, @view aij[i, :])  # ∂A/∂Nᵢ

        ∇P[i] = RT * (VmB⁻¹ + bᵢ * ΣnmolbyVmB²) - (
            (∂ᵢA * DmC - A * (dᵢ - cᵢ)) * VpC⁻¹mVpD⁻¹byDmC²
            + AbyDmC * (-cᵢ * VpC⁻² + dᵢ * VpD⁻²)
        )
    end
    ∇P[end] = - RT * ΣnmolbyVmB² + AbyDmC * (VpC⁻² - VpD⁻²)
    return nothing
end

struct HessianBuffer{T<:AbstractFloat}
    thermo::BrusilovskyThermoBuffer{T}
    matrnc::Matrix{T}
    vecnc₊::Vector{T}
    vecnc₁::Vector{T}
    vecnc₂::Vector{T}
end

function HessianBuffer(mix::BrusilovskyEoSMixture{T}) where {T<:Real}
    nc = ncomponents(mix)
    return HessianBuffer{T}(
        thermo_buffer(mix),
        Matrix{T}(undef, nc, nc),
        Vector{T}(undef, nc + 1),
        Vector{T}(undef, nc),
        Vector{T}(undef, nc),
    )
end

"""
Calculates hessian for VTFlash from `state` and base `nmol`, `volume`.
The `state` must be [N₁'/N₁, ..., Nₙ'/Nₙ, V'/V] vector.
"""
function __vt_flash_hessian!(
    hess::AbstractMatrix{T},
    state::AbstractVector{T},
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,  # I-state
    volume::Real,          # I-state
    RT::Real;
    buf::HessianBuffer=HessianBuffer(mix),
) where {T}
    # TODO: make hessian symmetric for optimization

    # tip: \bbB<Tab> for 𝔹 and so on

    #     |    |   |   M   size
    #     | 𝔹  | ℂ |   --------
    # ℍ = |    |   |   𝔹   n×n
    #     |----|---|   ℂ   n×1
    #     | ℂᵀ | 𝔻 |   𝔻   1×1

    #             [        Nⱼ²          ( ∂lnΦᵢ           ∂lnΦᵢ           ) ]
    # 𝔹ᵢⱼ = RT Nᵢ [ δᵢⱼ ---------- - Nⱼ ( -----(N', V') + -----(N'', V'') ) ]
    #             [      N'ⱼ N''ⱼ       (  ∂Nⱼ             ∂Nⱼ            ) ]
    N₁ = buf.vecnc₁
    N₁ .= nmol .* @view state[1:end-1]
    V₁ = volume * state[end]

    𝔹 = @view hess[1:end-1, 1:end-1]
    ∇P = buf.vecnc₊  # (n + 1) size
    ∇P⁻ = @view ∇P[1:end-1]  # n size

    # ∇P⁻ used as buffer
    # Initialization (!) of 𝔹 with jacobian'
    # 𝔹 = jacobian'
    log_c_activity_wj!(∇P⁻, 𝔹, mix, N₁, V₁, RT; buf=buf.thermo)

    N₂ = buf.vecnc₂
    N₂ .= nmol .- N₁
    V₂ = volume - V₁
    jacobian₂ = buf.matrnc

    # ∇P⁻ used as buffer
    # 𝔹 = jacobian' + jacobian''
    log_c_activity_wj!(∇P⁻, jacobian₂, mix, N₂, V₂, RT; buf=buf.thermo)
    𝔹 .+= jacobian₂

    # 𝔹 = - Nᵢ Nⱼ * (jacobian' + jacobian'')
    # the minus missed cuz of ln Φᵢ = -ln Cₐᵢ
    𝔹 .*= nmol .* nmol'

    # 𝔹, adding diagonal term
    @inbounds for i in eachindex(nmol)
        y₁ = state[i]      # N'ᵢ / Nᵢ
        y₂ = 1 - state[i]  # N''ᵢ / Nᵢ
        𝔹[i, i] += nmol[i] ./ (y₁ * y₂)
    end
    # final 𝔹
    𝔹 .*= RT

    #            [ ∂P             ∂P             ]
    # ℂᵢ = -V Nᵢ [ --- (N', V') + --- (N'', V'') ]
    #            [ ∂Nᵢ            ∂Nᵢ            ]
    #
    #         [ ∂P            ∂P            ]
    # 𝔻 = -V² [ -- (N', V') + -- (N'', V'') ]
    #         [ ∂V            ∂V            ]
    __vt_flash_pressure_gradient!(∇P, mix, N₁, V₁, RT; buf=buf.thermo)
    ℂ = @view hess[1:end-1, end]
    ℂ .= @view ∇P[1:end-1]  # ℂ = (∂P/∂Nᵢ)'
    𝔻 = ∇P[end]  # 𝔻 = (∂P/∂V)'

    __vt_flash_pressure_gradient!(∇P, mix, N₂, V₂, RT; buf=buf.thermo)
    ℂ .+= @view ∇P[1:end-1]  # ℂ = ∇P' + ∇P''
    ℂ .*= -volume .* nmol  # final ℂ

    # hess[[end], 1:end-1] .= ℂ'  # ℂᵀ part of hessian
    hess[end, 1:end-1] .= ℂ  # seems correct and no allocs

    𝔻 += ∇P[end]  # 𝔻 = (∂P/∂V)' + (∂P/∂V)''
    𝔻 *= -volume^2  # final 𝔻
    hess[end, end] = 𝔻
    return nothing
end

function vt_flash_closures(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real,
) where {T}
    N₁ = Vector{T}(undef, ncomponents(mix))
    N₂ = Vector{T}(undef, ncomponents(mix))
    log_cₐ₁ = Vector{T}(undef, ncomponents(mix))
    log_cₐ₂ = Vector{T}(undef, ncomponents(mix))

    thermo_buf = thermo_buffer(mix)

    # calculates once
    Pbase = pressure(mix, nmol, volume, RT; buf=thermo_buf)
    log_cₐ_base = Vector{T}(undef, ncomponents(mix))
    log_c_activity!(log_cₐ_base, mix, nmol, volume, RT; buf=thermo_buf)

    "Constant vector for covolume constrain. [Nᵢbᵢ..., -V]"
    covolumes_b̃ = [(c.b for c in components(mix))..., 1]
    covolumes_b̃[1:end-1] .*= nmol
    covolumes_b̃[end] *= -volume

    # for `transform` function
    Tr_matrix = Diagonal([nmol..., volume])
    state_tr = Vector{T}(undef, size(Tr_matrix, 1))

    "Updates `N₁`, `N₂`. Returns `state_tr`, `V₁`, `V₂` from `state`."
    function transform(state::AbstractVector{T})
        mul!(state_tr, Tr_matrix, state)
        N₁ .= @view state_tr[1:end-1]
        N₂ .= nmol .- N₁
        V₁ = state_tr[end]
        V₂ = volume - V₁
        return state_tr, V₁, V₂
    end

    function constrain_step(state::AbstractVector{T}, dir::AbstractVector{T})
        αm = T(Inf)
        # # positiveness constrain (`0 < state[i] + α * dir[i] < 1`)
        @inbounds for i in eachindex(state)
            if dir[i] > 0
                α = (1 - state[i]) / dir[i]
            elseif dir[i] < 0
                α = - state[i] / dir[i]
            else
                if 0 < state[i] < 1
                    continue
                else
                    error("VTFlash: constrain_step. Zero direction $i, but state[$i] = $(state[i])")
                end
            end
            if 0 < α < αm
                αm = α
            elseif α < 0
                @warn "constrain not meet for i = $i" state[i] dir[i] α
            end
        end

        # covolume constrain
        "Covolume edge by phase 1."
        αm_cov₁ = - dot(state, covolumes_b̃) / dot(dir, covolumes_b̃)
        "Covolume edge by phase 2."
        αm_cov₂ = (sum(covolumes_b̃) - dot(state, covolumes_b̃)) / dot(dir, covolumes_b̃)

        if dot(dir, covolumes_b̃) > 0
            if 0 < αm_cov₁ < αm
                αm = αm_cov₁
            end
            if αm < αm_cov₂
                @warn "Covolume constrain of phase 2 not meet others"
            end
        else
            if αm < αm_cov₁
                @warn "Covolume constrain of phase 1 not meet others"
            end
            if 0 < αm_cov₂ < αm
                αm = αm_cov₂
            end
        end
        if αm == T(Inf)
            error("VTFlash: constrain_step. Step was not found.")
        end
        return αm
    end

    function helmholtz_diff_grad!(state::AbstractVector{T}, grad::AbstractVector{T})
        _, V₁, V₂ = transform(state)
        log_c_activity!(log_cₐ₁, mix, N₁, V₁, RT; buf=thermo_buf)
        log_c_activity!(log_cₐ₂, mix, N₂, V₂, RT; buf=thermo_buf)

        @inbounds for i in 1:length(state)-1
            Δμ = -RT * (log((N₂[i]/V₂) / (N₁[i]/V₁)) - (log_cₐ₁[i] - log_cₐ₂[i]))
            grad[i] = nmol[i] * Δμ
        end
        P₁ = pressure(mix, N₁, V₁, RT; buf=thermo_buf)
        P₂ = pressure(mix, N₂, V₂, RT; buf=thermo_buf)
        grad[end] = volume * (-P₁ + P₂)
        return grad
    end
    function helmholtz_diff!(state::AbstractVector{T}, grad::AbstractVector{T})
        _, V₁, V₂ = transform(state)

        log_c_activity!(log_cₐ₂, mix, N₂, V₂, RT; buf=thermo_buf)

        "Σ Nᵢ (μᵢ - μ₂ᵢ)"
        Ndotμ₂ = zero(T)
        @inbounds for i in 1:length(state)-1
            # μ base - μ₂
            Δμ = -RT * (log((N₂[i]/V₂)/(nmol[i]/volume)) - (log_cₐ_base[i] - log_cₐ₂[i]))
            Ndotμ₂ += nmol[i] * Δμ
        end

        P₂ = pressure(mix, N₂, V₂, RT; buf=thermo_buf)
        helmholtz_diff_grad!(state, grad)  # overwrites gradient `grad`
        ΔA = dot(grad, state) + (Pbase - P₂) * volume - Ndotμ₂
        @debug "helmholtz_diff!" state=repr(state) ΔA grad=repr(grad) norm(grad, 2)
        return ΔA, grad
    end
    return constrain_step, helmholtz_diff_grad!, helmholtz_diff!
end

"Find initial state by reducing saturation."
function __vt_flash_initial_state(
    mix::BrusilovskyEoSMixture,
    nmol::AbstractVector,
    volume::Real,
    RT::Real,
    stability_tries::AbstractVector{VTStabilityResult{T}};
    sat₁max::Real=0.25,
    steps::Integer=200,
    step_scale::Real=0.5,
    helmholtz_thresh::Real=-1e-7,
) where {T}
    constrain_step, _, helmholtz_diff! = vt_flash_closures(mix, nmol, volume, RT)

    state = Vector{T}(undef, ncomponents(mix) + 1)

    # choosing concentration with minimum of helmholtz density
    conc_test = __vt_flash_init_conc_choose(stability_tries)

    init_found = __vt_flash_initial_state!(
        state, nmol, volume, conc_test, helmholtz_diff!, constrain_step;
        sat₁max=sat₁max,
        steps=steps,
        step_scale=step_scale,
        helmholtz_thresh=helmholtz_thresh,
    )
    return init_found, state
end

function __vt_flash_initial_state!(
    state::AbstractVector{T},
    nmol::AbstractVector{T},
    volume::Real,
    conc₁::AbstractVector{T},
    helmholtz_diff!::Function,
    constrain_step::Function;
    sat₁max::Real=T(0.9),
    steps::Int=200,
    step_scale::Real=T(0.5),
    helmholtz_thresh::Real=T(-1e-7),  # must be negative value
) where {T}
    state[1:end-1] .= conc₁ * (sat₁max * volume) ./ nmol
    state[end] = sat₁max

    vec = similar(state)  # buffer vector for gradient
    scale = one(T)
    @debug "Initial state search" start_scale=scale sat₁max
    for i in 1:steps
        # upd `state`
        sat = sat₁max * scale
        state[1:end-1] .= conc₁ * (sat * volume) ./ nmol
        state[end] = sat
        # TODO: check if state feasible

        # calc helmholtz energy
        @debug "Initial state search" i state=repr(state) scale
        try
            ΔA, _ = helmholtz_diff!(state, vec)
            # check convergence
            @debug "Initial state search: ΔA calculated" i ΔA helmholtz_thresh
            if ΔA < helmholtz_thresh
                return true
            end
        catch e
            # @warn "VTFlash: initial state search" sat e
        end

        # update `scale`
        scale *= step_scale
    end
    return false
end

"""
Extracts vt-state from `optresult` (DescentMethods obj).
Sorts variables into gas and liquid.
Returns corresponding `VTFlashResult`.
"""
function __vt_flash_two_phase_result(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector{T},
    volume::Real,
    RT::Real,
    optresult,
) where {T}
    # □₁ for gas, □₂ for liquid
    state = optresult.argument
    nmol₁ = nmol .* @view state[1:end-1]
    V₁ = volume * state[end]
    nmol₂ = nmol .- nmol₁
    V₂ = volume - V₁

    P₁ = pressure(mix, nmol₁, V₁, RT)  # they should be equal
    P₂ = pressure(mix, nmol₂, V₂, RT)

    Z₁ = P₁ * V₁ / (sum(nmol₁) * RT)  # seems can be reduced to Vᵢ / sum(nmolᵢ)
    Z₂ = P₂ * V₂ / (sum(nmol₂) * RT)

    if Z₂ > Z₁  # □₂ is gas state, need exchange
        P₁, P₂ = P₂, P₁
        Z₁, Z₂ = Z₂, Z₁
        V₁, V₂ = V₂, V₁

        for i in eachindex(nmol₁, nmol₂)
            nmol₁[i], nmol₂[i] = nmol₂[i], nmol₁[i]
        end
    end

    return VTFlashResult{T}(;
            converged=optresult.converged,
            singlephase=false,
            RT=RT,
            nmol_1=nmol₁,
            V_1=V₁,
            nmol_2=nmol₂,
            V_2=V₂,
            iters=optresult.iterations,
            fcalls=optresult.calls,
    )
end

"Perform vt_flash from `unstable_state`, so only phase splitting is done."
function vt_flash(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real,
    unstable_state::AbstractVector,
) where {T}
    state = copy(unstable_state)

    # initial hessian
    hessian = Matrix{T}(undef, (size(state, 1), size(state, 1)))
    __vt_flash_hessian!(hessian, state, mix, nmol, volume, RT)

    # create closures for helmoltz energy, its gradient and constrain step
    constrain_step, _, helmholtz_diff! = vt_flash_closures(mix, nmol, volume, RT)

    # run optimizer
    optmethod = DescentMethods.CholBFGS(state)
    DescentMethods.reset!(optmethod, state, hessian)
    result = DescentMethods.optimize!(optmethod, helmholtz_diff!, state;
        gtol=1e-3,
        maxiter=100,
        constrain_step=constrain_step,
        reset=false,
    )
    return __vt_flash_two_phase_result(mix, nmol, volume, RT, result)
end

"VT-flash."
function vt_flash(
    mix::BrusilovskyEoSMixture{T},
    nmol::AbstractVector,
    volume::Real,
    RT::Real,
) where {T}
    # run vt-stability to find out whether a state single phase or not
    singlephase, vt_stab_tries = vt_stability(mix, nmol, volume, RT)

    @debug "VTFlash: VTStability result" singlephase

    if singlephase
        return VTFlashResult{T}(;
            converged=true,
            singlephase=true,
            RT=RT,
            nmol_1=nmol,
            V_1=volume,
            nmol_2=similar(nmol),
            V_2=0,
        )
    end

    # two-phase state case
    init_found, state = __vt_flash_initial_state(
        mix, nmol, volume, RT, vt_stab_tries
    )

    @debug "VTFlash: initial state search result" found=init_found state=repr(state) ΔA=helmholtz_diff!(state, similar(state))[1]

    if !init_found
        @error "VTFlash: Initial state was not found!" mixture=mix nmol=repr(nmol) volume=volume RT=RT
        error("VTFlash: Initial state was not found!")
    end

    return vt_flash(mix, nmol, volume, RT, state)
end

max_stars_repo_name: vvpisarev/CubicEoS

max_issues_repo_name: vvpisarev/CubicEoS

max_forks_repo_name: vvpisarev/CubicEoS


text:
using Jive
runtests(@__DIR__)


text:
using Test
using TOMLX
using TOMLX: Internals

macro testval(s, v)
    f = "foo = $s"
    :( @test TOMLX.parsestring($f)["foo"] == $v )
end

macro failval(s, v)
    f = "foo = $s"
    :( err = TOMLX.tryparsestring($f);
       @test err isa TOMLX.Internals.ParserError;
       @test err.type == $v;
    )
end

@testset "Numbers" begin
    @failval("00"                   , Internals.ErrParsingDateTime)
    @failval("-00"                  , Internals.ErrParsingDateTime)
    @failval("+00"                  , Internals.ErrParsingDateTime)
    @failval("00.0"                 , Internals.ErrParsingDateTime)
    @failval("-00.0"                , Internals.ErrParsingDateTime)
    @failval("+00.0"                , Internals.ErrParsingDateTime)
    @failval("9223372036854775808"  , Internals.ErrOverflowError)
    @failval("-9223372036854775809" , Internals.ErrOverflowError)

    @failval("0."        , Internals.ErrNoTrailingDigitAfterDot)
    @failval("0.e"       , Internals.ErrNoTrailingDigitAfterDot)
    @failval("0.E"       , Internals.ErrNoTrailingDigitAfterDot)
    @failval("0.0E"      , Internals.ErrGenericValueError)
    @failval("0.0e"      , Internals.ErrGenericValueError)
    @failval("0.0e-"     , Internals.ErrGenericValueError)
    @failval("0.0e+"     , Internals.ErrGenericValueError)
    # @failval("0.0e+00" , Internals.ErrGenericValueError)

    @testval("1.0"         , 1.0)
    @testval("1.0e0"       , 1.0)
    @testval("1.0e+0"      , 1.0)
    @testval("1.0e-0"      , 1.0)
    @testval("1.001e-0"    , 1.001)
    @testval("2e10"        , 2e10)
    @testval("2e+10"       , 2e10)
    @testval("2e-10"       , 2e-10)
    @testval("2_0.0"       , 20.0)
    @testval("2_0.0_0e0_0" , 20.0)
    @testval("2_0.1_0e1_0" , 20.1e10)

    @testval("1_0"    , 10)
    @testval("1_0_0"  , 100)
    @testval("1_000"  , 1000)
    @testval("+1_000" , 1000)
    @testval("-1_000" , -1000)

    @failval("0_"     , Internals.ErrLeadingZeroNotAllowedInteger)
    @failval("0__0"   , Internals.ErrLeadingZeroNotAllowedInteger)
    @failval("__0"    , Internals.ErrUnexpectedStartOfValue)
    @failval("1_0_"   , Internals.ErrTrailingUnderscoreNumber)
    @failval("1_0__0" , Internals.ErrUnderscoreNotSurroundedByDigits)
end

@testset "Booleans" begin
    @testval("true", true)
    @testval("false", false)

    @failval("true2"  , Internals.ErrExpectedNewLineKeyValue)
    @failval("false2" , Internals.ErrExpectedNewLineKeyValue)
    @failval("talse"  , Internals.ErrGenericValueError)
    @failval("frue"   , Internals.ErrGenericValueError)
    @failval("t1"     , Internals.ErrGenericValueError)
    @failval("f1"     , Internals.ErrGenericValueError)
end

@testset "Datetime" begin
    @testval("2016-09-09T09:09:09"     , DateTime(2016 , 9 , 9 , 9 , 9 , 9))
    @testval("2016-09-09T09:09:09Z"    , DateTime(2016 , 9 , 9 , 9 , 9 , 9))
    @testval("2016-09-09T09:09:09.0Z"  , DateTime(2016 , 9 , 9 , 9 , 9 , 9))
    @testval("2016-09-09T09:09:09.012" , DateTime(2016 , 9 , 9 , 9 , 9 , 9  , 12))

    @failval("2016-09-09T09:09:09.0+10:00"   , Internals.ErrOffsetDateNotSupported)
    @failval("2016-09-09T09:09:09.012-02:00" , Internals.ErrOffsetDateNotSupported)
    @failval("2016-09-09T09:09:09.0+10:00"   , Internals.ErrOffsetDateNotSupported)
    @failval("2016-09-09T09:09:09.012-02:00" , Internals.ErrOffsetDateNotSupported)

    @failval("2016-09-09T09:09:09.Z" , Internals.ErrParsingDateTime)
    @failval("2016-9-09T09:09:09Z"   , Internals.ErrParsingDateTime)
    @failval("2016-13-09T09:09:09Z"  , Internals.ErrParsingDateTime)
    @failval("2016-02-31T09:09:09Z"  , Internals.ErrParsingDateTime)
    @failval("2016-09-09T09:09:09x"  , Internals.ErrParsingDateTime)
    @failval("2016-09-09s09:09:09Z"  , Internals.ErrParsingDateTime)
    @failval("2016-09-09T09:09:09x"  , Internals.ErrParsingDateTime)
end

@testset "Time" begin
    @testval("09:09:09.99"    , Time(9 , 9 , 9 , 99))
    @testval("09:09:09.99999" , Time(9 , 9 , 9 , 999))

    @failval("09:09x09", Internals.ErrParsingDateTime)
end

# TODO: Add more dedicated value tests

@testset "String" begin
    @failval("\"foooo", Internals.ErrUnexpectedEndString)

    #=
    Found these examples of string tests somewhere
    quot0=""" """    # valid
    quot1=""" """"   # valid
    quot2=""" """""  # valid
    quot3=""" """""" # invalid
    apos0=''' '''    # valid
    apos1=''' ''''   # valid
    apos2=''' '''''  # valid
    apos3=''' '''''' # invalid

    quot4=""""""     # valid (6 in a row, empty string)
    quot5="""" """   # valid
    quot6=""""" """  # valid
    quot7="""""" """ # invalid
    apos4=''''''     # valid (6 in a row, empty string)
    apos5='''' '''   # valid
    apos6=''''' '''  # valid
    apos7='''''' ''' # invalid

    quot8="""""\""""""     # valid (5 at start, 6 at end)
    quot9="""""\"""\"""""" # valid (3 in the middle 5 at start, 6 at end)

    =#
end