text:
using WignerDMatrices
using HalfIntegers
using LinearAlgebra
using Test

import WignerDMatrices: cis_special, Piby2Radians, ZeroRadians,
TwoNPiRadians, TwoNPlusOnePiRadians
import WignerDMatrices: WignerdMatrixContainer, flatind, filledelements,
flatind_phase, sphericaldegree, SpinMatrix

@test isempty(Test.detect_ambiguities(Base, Core, WignerDMatrices))

@testset "special points" begin
	@testset "float" begin
		@testset "Piby2Radians" begin
		    @test float(Piby2Radians()) == π/2
	    	@test AbstractFloat(Piby2Radians()) == π/2
	    	@test Float64(Piby2Radians()) == π/2
	    	@test float(Piby2Radians{BigFloat}()) == big(π)/2
	    	@test Float64(Piby2Radians{BigFloat}()) == π/2
		end
		@testset "ZeroRadians" begin
			@test float(ZeroRadians()) == 0
	    	@test AbstractFloat(ZeroRadians()) == 0
	    	@test Float64(ZeroRadians()) == 0
	    	@test float(ZeroRadians{BigFloat}()) isa BigFloat
	    	@test float(ZeroRadians{BigFloat}()) == big(0.0)
		end
		@testset "TwoNPiRadians" begin
		    @test float(TwoNPiRadians(1)) == 2π
		    @test float(TwoNPiRadians{BigFloat}(1)) == big(2)*π
		    @test float(TwoNPiRadians(-1)) == -2π
		    @test float(TwoNPiRadians{BigFloat}(-1)) == big(-2)*π
		end
		@testset "TwoNPlusOnePiRadians" begin
		   	@test float(TwoNPlusOnePiRadians(1)) == 3π
		    @test float(TwoNPlusOnePiRadians{BigFloat}(1)) == big(3)*π
		   	@test float(TwoNPlusOnePiRadians(-1)) == -π
		    @test float(TwoNPlusOnePiRadians{BigFloat}(-1)) == big(-1)*π
		end
		@testset "TwoNPiRadians and ZeroRadians" begin
		    @test float(TwoNPiRadians(0)) == float(ZeroRadians())
		    @test float(TwoNPiRadians{BigFloat}(0)) == float(ZeroRadians{BigFloat}())
		end
	end
	@testset "promote" begin
		for T in [Int,Float32,Float64]
		    @test promote_rule(Piby2Radians{Float64},T) == Float64
		    @test promote_rule(ZeroRadians{Float64},T) == Float64
		end
		for T in [Int,Float32,Float64]
		    @test promote_rule(Piby2Radians{BigFloat},T) == BigFloat
		    @test promote_rule(ZeroRadians{BigFloat},T) == BigFloat
		end
		for T in [BigFloat]
			@test promote_rule(Piby2Radians{Float64},T) == BigFloat
			@test promote_rule(ZeroRadians{Float64},T) == BigFloat
		end
	end
	@testset "zero and one" begin
		function test(T, n...)
			@test zero(T(n...)) == zero(Float64)
	    	@test one(T(n...)) == one(Float64)
		    @test zero(T{BigFloat}(n...)) isa BigFloat
		    @test zero(T{BigFloat}(n...)) == zero(BigFloat)
	    	@test one(T{BigFloat}(n...)) isa BigFloat
	    	@test one(T{BigFloat}(n...)) == one(BigFloat)
		end
		@testset "Piby2Radians" begin
		    test(Piby2Radians)
		end
		@testset "ZeroRadians" begin
			test(ZeroRadians)
		end
		@testset "TwoNPiRadians" begin
		   	test(TwoNPiRadians, 1)
		end
		@testset "TwoNPlusOnePiRadians" begin
		   	test(TwoNPlusOnePiRadians, 1)
		end
	end
	@testset "negate" begin
	    @test -ZeroRadians() == ZeroRadians()
	    @test -Piby2Radians() == -π/2
	    @test -Piby2Radians{BigFloat}() == -big(π)/2
	    @test -TwoNPlusOnePiRadians{BigFloat}(0) == -big(π)
	    @test -TwoNPiRadians(1) == -2π
	    @test -TwoNPiRadians{BigFloat}(1) == -2big(π)
	    @test -TwoNPlusOnePiRadians(1) == -3π
	    @test -TwoNPlusOnePiRadians{BigFloat}(1) == -3big(π)
	end
	@testset "add" begin
		@testset "ZeroRadians" begin
			@test ZeroRadians() + ZeroRadians() == ZeroRadians()
	    	@test ZeroRadians{Float64}() + ZeroRadians{BigFloat}() === ZeroRadians{BigFloat}()
	    	@test ZeroRadians{BigFloat}() + ZeroRadians{Float64}() === ZeroRadians{BigFloat}()

	    	@test ZeroRadians() + Piby2Radians() == Piby2Radians()
	    	@test ZeroRadians() + Piby2Radians{BigFloat}() == Piby2Radians{BigFloat}()
	    	@test ZeroRadians{BigFloat}() + Piby2Radians() == Piby2Radians{BigFloat}()

	    	@test ZeroRadians() + TwoNPiRadians(1) == TwoNPiRadians(1)
	    	@test ZeroRadians() + TwoNPiRadians{BigFloat}(1) == TwoNPiRadians{BigFloat}(1)
	    	@test ZeroRadians{BigFloat}() + TwoNPiRadians(1) == TwoNPiRadians{BigFloat}(1)

	    	@test ZeroRadians() + TwoNPlusOnePiRadians(1) == TwoNPlusOnePiRadians(1)
	    	@test ZeroRadians() + TwoNPlusOnePiRadians{BigFloat}(1) == TwoNPlusOnePiRadians{BigFloat}(1)
	    	@test ZeroRadians{BigFloat}() + TwoNPlusOnePiRadians(1) == TwoNPlusOnePiRadians{BigFloat}(1)
		end
		@testset "Piby2Radians" begin
			@test Piby2Radians() + Piby2Radians() == TwoNPlusOnePiRadians(0)
	    	@test Piby2Radians{Float64}() + Piby2Radians{BigFloat}() === TwoNPlusOnePiRadians{BigFloat}(0)
	    	@test Piby2Radians{BigFloat}() + Piby2Radians{Float64}() === TwoNPlusOnePiRadians{BigFloat}(0)

	    	@test Piby2Radians() + ZeroRadians() == Piby2Radians()
	    	@test Piby2Radians{BigFloat}() + ZeroRadians() == Piby2Radians{BigFloat}()
	    	@test Piby2Radians() + ZeroRadians{BigFloat}() == Piby2Radians{BigFloat}()

	    	@test Piby2Radians() + TwoNPiRadians(1) == 5π/2
	    	@test Piby2Radians() + TwoNPiRadians{BigFloat}(1) == 5big(π)/2
	    	@test Piby2Radians{BigFloat}() + TwoNPiRadians(1) == 5big(π)/2

	    	@test Piby2Radians() + TwoNPlusOnePiRadians(1) == 7π/2
	    	@test Piby2Radians{BigFloat}() + TwoNPlusOnePiRadians(1) == 7big(π)/2
	    	@test Piby2Radians() + TwoNPlusOnePiRadians{BigFloat}(1) == 7big(π)/2
		end
		@testset "TwoNPiRadians" begin
		    @test TwoNPiRadians(1) + TwoNPiRadians(1) == TwoNPiRadians(2)
	    	@test TwoNPiRadians{Float64}(1) + TwoNPiRadians{BigFloat}(1) === TwoNPiRadians{BigFloat}(2)
	    	@test TwoNPiRadians{BigFloat}(1) + TwoNPiRadians{Float64}(1) === TwoNPiRadians{BigFloat}(2)

	    	@test TwoNPiRadians(1) + ZeroRadians() == TwoNPiRadians(1)
	    	@test TwoNPiRadians{BigFloat}(1) + ZeroRadians() == TwoNPiRadians{BigFloat}(1)
	    	@test TwoNPiRadians(1) + ZeroRadians{BigFloat}() == TwoNPiRadians{BigFloat}(1)

	    	@test TwoNPiRadians(1) + Piby2Radians() == 5π/2
	    	@test TwoNPiRadians{BigFloat}(1) + Piby2Radians() == 5big(π)/2
	    	@test TwoNPiRadians(1) + Piby2Radians{BigFloat}() == 5big(π)/2

	    	@test TwoNPiRadians(1) + TwoNPlusOnePiRadians(0) == TwoNPlusOnePiRadians(1)
	    	@test TwoNPiRadians{BigFloat}(1) + TwoNPlusOnePiRadians(0) == TwoNPlusOnePiRadians{BigFloat}(1)
	    	@test TwoNPiRadians(1) + TwoNPlusOnePiRadians{BigFloat}(0) == TwoNPlusOnePiRadians{BigFloat}(1)

	    	@test TwoNPiRadians(1) + TwoNPlusOnePiRadians(1) == TwoNPlusOnePiRadians(2)
	    	@test TwoNPiRadians{BigFloat}(1) + TwoNPlusOnePiRadians(1) == TwoNPlusOnePiRadians{BigFloat}(2)
	    	@test TwoNPiRadians(1) + TwoNPlusOnePiRadians{BigFloat}(1) == TwoNPlusOnePiRadians{BigFloat}(2)
		end
		@testset "TwoNPlusOnePiRadians" begin
		    @test TwoNPlusOnePiRadians(1) + TwoNPlusOnePiRadians(1) == TwoNPiRadians(3)
	    	@test TwoNPlusOnePiRadians{Float64}(1) + TwoNPlusOnePiRadians{BigFloat}(1) === TwoNPiRadians{BigFloat}(3)
	    	@test TwoNPlusOnePiRadians{BigFloat}(1) + TwoNPlusOnePiRadians{Float64}(1) === TwoNPiRadians{BigFloat}(3)

	    	@test TwoNPlusOnePiRadians(1) + ZeroRadians() == TwoNPlusOnePiRadians(1)
	    	@test TwoNPlusOnePiRadians{BigFloat}(1) + ZeroRadians() == TwoNPlusOnePiRadians{BigFloat}(1)
	    	@test TwoNPlusOnePiRadians(1) + ZeroRadians{BigFloat}() == TwoNPlusOnePiRadians{BigFloat}(1)

	    	@test TwoNPlusOnePiRadians(1) + Piby2Radians() == 7π/2
	    	@test TwoNPlusOnePiRadians{BigFloat}(1) + Piby2Radians() == 7big(π)/2
	    	@test TwoNPlusOnePiRadians(1) + Piby2Radians{BigFloat}() == 7big(π)/2

	    	@test TwoNPlusOnePiRadians(1) + TwoNPlusOnePiRadians(0) == TwoNPiRadians(2)
	    	@test TwoNPlusOnePiRadians{BigFloat}(1) + TwoNPlusOnePiRadians(0) == TwoNPiRadians{BigFloat}(2)
	    	@test TwoNPlusOnePiRadians(1) + TwoNPlusOnePiRadians{BigFloat}(0) == TwoNPiRadians{BigFloat}(2)

	    	@test TwoNPlusOnePiRadians(1) + TwoNPiRadians(1) == TwoNPlusOnePiRadians(2)
	    	@test TwoNPlusOnePiRadians(1) + TwoNPiRadians{BigFloat}(1) == TwoNPlusOnePiRadians{BigFloat}(2)
	    	@test TwoNPlusOnePiRadians{BigFloat}(1) + TwoNPiRadians(1) == TwoNPlusOnePiRadians{BigFloat}(2)
		end
	end
	@testset "trigonometric functions" begin
	    @testset "Piby2Radians" begin
		    for α = -10:1//2:10
		    	@test cis_special(HalfInt(α),Piby2Radians()) isa ComplexF64
		    	@test cis_special(HalfInt(α),Piby2Radians()) ≈ cis(α*π/2)
		    	@test cis_special(HalfInt(α),Piby2Radians{BigFloat}()) isa Complex{BigFloat}
		    	@test cis_special(HalfInt(α),Piby2Radians{BigFloat}()) ≈ cis(α*big(π)/2)
		    end
	    	@test cos(Piby2Radians()) == 0
	    	@test cos(Piby2Radians{BigFloat}()) isa BigFloat
	    	@test cos(Piby2Radians{BigFloat}()) == 0
	    	@test sin(Piby2Radians()) == 1
	    	@test sin(Piby2Radians{BigFloat}()) isa BigFloat
	    	@test sin(Piby2Radians{BigFloat}()) == 1
		end
		@testset "ZeroRadians" begin
			for α = -10:1//2:10
		    	@test cis_special(HalfInt(α),ZeroRadians()) isa ComplexF64
		    	@test cis_special(HalfInt(α),ZeroRadians()) ≈ cis(0)
		    	@test cis_special(HalfInt(α),ZeroRadians{BigFloat}()) isa Complex{BigFloat}
		    	@test cis_special(HalfInt(α),ZeroRadians{BigFloat}()) ≈ cis(big(0))
		    end
		    @test cos(ZeroRadians()) == 1
		    @test cos(ZeroRadians{BigFloat}()) isa BigFloat
		    @test cos(ZeroRadians{BigFloat}()) == 1
	    	@test sin(ZeroRadians()) == 0
	    	@test sin(ZeroRadians{BigFloat}()) isa BigFloat
	    	@test sin(ZeroRadians{BigFloat}()) == 0
		end
		@testset "TwoNPiRadians" begin
		    for α = -10:1//2:10
		    	@test cis_special(HalfInt(α),TwoNPiRadians(1)) isa ComplexF64
		    	@test cis_special(HalfInt(α),TwoNPiRadians(1)) ≈ cis(α*2π)
		    	@test cis_special(HalfInt(α),TwoNPiRadians{BigFloat}(1)) isa Complex{BigFloat}
		    	@test cis_special(HalfInt(α),TwoNPiRadians{BigFloat}(1)) ≈ cis(α*2big(π))
		    end
		    for n = -3:3
			    @test cos(TwoNPiRadians(n)) == 1
			    @test cos(TwoNPiRadians{BigFloat}(n)) isa BigFloat
			    @test cos(TwoNPiRadians{BigFloat}(n)) == 1
		    	@test sin(TwoNPiRadians(n)) == 0
		    	@test sin(TwoNPiRadians{BigFloat}(n)) isa BigFloat
		    	@test sin(TwoNPiRadians{BigFloat}(n)) == 0
		    end
		end
		@testset "TwoNPlusOnePiRadians" begin
		    for α = -10:1//2:10
		    	@test cis_special(HalfInt(α),TwoNPlusOnePiRadians(1)) isa ComplexF64
		    	@test cis_special(HalfInt(α),TwoNPlusOnePiRadians(1)) ≈ cis(α*3π)
		    	@test cis_special(HalfInt(α),TwoNPlusOnePiRadians{BigFloat}(1)) isa Complex{BigFloat}
		    	@test cis_special(HalfInt(α),TwoNPlusOnePiRadians{BigFloat}(1)) ≈ cis(α*3big(π))
		    end
		    for n = -3:3
			    @test cos(TwoNPlusOnePiRadians(n)) == -1
			    @test cos(TwoNPlusOnePiRadians{BigFloat}(n)) isa BigFloat
			    @test cos(TwoNPlusOnePiRadians{BigFloat}(n)) == -1
		    	@test sin(TwoNPlusOnePiRadians(n)) == 0
		    	@test sin(TwoNPlusOnePiRadians{BigFloat}(n)) isa BigFloat
		    	@test sin(TwoNPlusOnePiRadians{BigFloat}(n)) == 0
		    end
		end
	end
	@testset "wignerdmatrixelement" begin
		function testapprox(m,n,dj_m_n,dj_m_n2)
	    	@test begin
	    		res = isapprox(dj_m_n,dj_m_n2,atol=1e-13,rtol=sqrt(eps(Float64)))
	    		if !res
	    			@show m n dj_m_n dj_m_n2
	    		end
	    		res
	    	end
	    end
	    function tests(j,v)
	    	@testset "Piby2Radians" begin
		        for m in -j:j, n in -j:j
		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,π/2 ,v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,Piby2Radians(),v)

		        	testapprox(m,n,dj_m_n,dj_m_n2)
		        end
		    end
		    @testset "ZeroRadians" begin
		        for m in -j:j, n in -j:j
		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,0,v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,ZeroRadians(),v)

		        	testapprox(m,n,dj_m_n,dj_m_n2)
		        end
		    end
		    @testset "TwoNPiRadians" begin
		        for m in -j:j, n in -j:j
		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,2π,v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,TwoNPiRadians(1),v)

		        	testapprox(m,n,dj_m_n,dj_m_n2)

		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,0,v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,TwoNPiRadians(0),v)

		        	testapprox(m,n,dj_m_n,dj_m_n2)

		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,ZeroRadians(),v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,TwoNPiRadians(0),v)

		        	testapprox(m,n,dj_m_n,dj_m_n2)
		        end
		    end
		    @testset "TwoNPlusOnePiRadians" begin
		        for m in -j:j, n in -j:j
		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,3π,v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,WignerDMatrices.ThreePi,v)

		        	dj_m_n = WignerDMatrices.wignerdmatrixelement(j,m,n,π,v)
		        	dj_m_n2 = WignerDMatrices.wignerdmatrixelement(j,m,n,WignerDMatrices.Pi,v)

		        	testapprox(m,n,dj_m_n,dj_m_n2)
		        end
		    end
		end

		@testset "integer j" begin
		    for j = 0:10
				A = zeros(ComplexF64,2j+1,2j+1)
			    v = WignerDMatrices.eigenvecsJy!(j,A)

			    tests(j,v)
			end
		end
		@testset "half-integer j" begin
			for j = 1//2:19//2
				A = zeros(ComplexF64,Int(2j+1),Int(2j+1))
			    v = WignerDMatrices.eigenvecsJy!(HalfInt(j),A)

		    	tests(j,v)
		    end
		end
	end
end
@testset "WignerdMatrixContainer" begin
    @testset "flatind" begin
    	@testset "j = 1/2" begin
    	    j = HalfInt(1/2)
    	    @test filledelements(j) == 2
    	    @test flatind(j, 1, 1) == 1
            @test flatind(j, 2, 1) == 2
    	end
    	@testset "j = 1" begin
        	j = 1
        	@test filledelements(j) == 4
            @test flatind(j, 1, 1) == 1
            @test flatind(j, 2, 1) == 2
            @test flatind(j, 3, 1) == 3
            @test flatind(j, 2, 2) == 4
        end
    	@testset "j = 3/2" begin
        	j = HalfInt(3/2)
        	@test filledelements(j) == 6
            @test flatind(j, 1, 1) == 1
            @test flatind(j, 2, 1) == 2
            @test flatind(j, 3, 1) == 3
            @test flatind(j, 4, 1) == 4
            @test flatind(j, 2, 2) == 5
            @test flatind(j, 3, 2) == 6
        end
        @testset "j = 2" begin
        	j = 2
        	@test filledelements(j) == 9
            @test flatind(j, 1, 1) == 1
            @test flatind(j, 2, 1) == 2
            @test flatind(j, 3, 1) == 3
            @test flatind(j, 4, 1) == 4
            @test flatind(j, 5, 1) == 5
            @test flatind(j, 2, 2) == 6
            @test flatind(j, 3, 2) == 7
            @test flatind(j, 4, 2) == 8
            @test flatind(j, 3, 3) == 9
        end
    end
    @testset "flating phase" begin
        @testset "j = 1/2" begin
			j = HalfInt(1/2)
			@testset "saved" begin
				@test flatind_phase(j, 1, 1) == (1, 1)
				@test flatind_phase(j, 2, 1) == (2, 1)
			end
			@testset "right" begin
				@test flatind_phase(j, 1, 2) == (2, -1)
				@test flatind_phase(j, 2, 2) == (1, 1)
			end
    	end
    	@testset "j = 1" begin
        	j = 1
        	@testset "saved" begin
	            @test flatind_phase(j, 1, 1) == (1, 1)
	            @test flatind_phase(j, 2, 1) == (2, 1)
	            @test flatind_phase(j, 3, 1) == (3, 1)
            	@test flatind_phase(j, 2, 2) == (4, 1)
        	end
        	@testset "top" begin
            	@test flatind_phase(j, 1, 2) == (2, -1)
        	end
        	@testset "bottom" begin
            	@test flatind_phase(j, 3, 2) == (2, 1)
        	end
        	@testset "right" begin
            	@test flatind_phase(j, 1, 3) == (3, 1)
	            @test flatind_phase(j, 2, 3) == (2, -1)
	            @test flatind_phase(j, 3, 3) == (1, 1)
        	end
        end
    	@testset "j = 3/2" begin
        	j = HalfInt(3/2)
        	@testset "saved" begin
	        	@test flatind_phase(j, 1, 1) == (1, 1)
	            @test flatind_phase(j, 2, 1) == (2, 1)
	            @test flatind_phase(j, 3, 1) == (3, 1)
	            @test flatind_phase(j, 4, 1) == (4, 1)
	            @test flatind_phase(j, 2, 2) == (5, 1)
	            @test flatind_phase(j, 3, 2) == (6, 1)
        	end
        	@testset "top" begin
            	@test flatind_phase(j, 1, 2) == (2, -1)
            	@test flatind_phase(j, 1, 3) == (3, 1)
        	end
        	@testset "bottom" begin
            	@test flatind_phase(j, 4, 2) == (3, 1)
            	@test flatind_phase(j, 4, 3) == (2, 1)
        	end
        	@testset "right" begin
	            @test flatind_phase(j, 2, 3) == (6, -1)
	            @test flatind_phase(j, 3, 3) == (5, 1)
	            @test flatind_phase(j, 1, 4) == (4, -1)
	            @test flatind_phase(j, 2, 4) == (3, 1)
	            @test flatind_phase(j, 3, 4) == (2, -1)
	            @test flatind_phase(j, 4, 4) == (1, 1)
        	end
        end
        @testset "j = 2" begin
        	j = 2
        	@testset "saved" begin
	        	@test flatind_phase(j, 1, 1) == (1, 1)
	            @test flatind_phase(j, 2, 1) == (2, 1)
	            @test flatind_phase(j, 3, 1) == (3, 1)
	            @test flatind_phase(j, 4, 1) == (4, 1)
	            @test flatind_phase(j, 5, 1) == (5, 1)
	            @test flatind_phase(j, 2, 2) == (6, 1)
	            @test flatind_phase(j, 3, 2) == (7, 1)
	            @test flatind_phase(j, 4, 2) == (8, 1)
            	@test flatind_phase(j, 3, 3) == (9, 1)
        	end
        	@testset "top" begin
            	@test flatind_phase(j, 1, 2) == (2, -1)
            	@test flatind_phase(j, 1, 3) == (3, 1)
            	@test flatind_phase(j, 1, 4) == (4, -1)
            	@test flatind_phase(j, 2, 3) == (7, -1)
        	end
        	@testset "bottom" begin
            	@test flatind_phase(j, 5, 2) == (4, 1)
            	@test flatind_phase(j, 5, 3) == (3, 1)
            	@test flatind_phase(j, 5, 4) == (2, 1)
            	@test flatind_phase(j, 4, 3) == (7, 1)
        	end
        	@testset "right" begin
	            @test flatind_phase(j, 2, 4) == (8, 1)
	            @test flatind_phase(j, 3, 4) == (7, -1)
	            @test flatind_phase(j, 4, 4) == (6, 1)
	            @test flatind_phase(j, 1, 5) == (5, 1)
	            @test flatind_phase(j, 2, 5) == (4, -1)
	            @test flatind_phase(j, 3, 5) == (3, 1)
	            @test flatind_phase(j, 4, 5) == (2, -1)
	            @test flatind_phase(j, 5, 5) == (1, 1)
        	end
        end
    end
end
@testset "Jy" begin
	jmax = 2
    A = zeros(ComplexF64,2jmax+1,2jmax+1)

    @testset "Jy_matrix_zbasis!" begin
    	j = 0
	    h = WignerDMatrices.Jy_matrix_zbasis!(j,A)
	    @test size(h) == (2j+1,2j+1)
	    @test ishermitian(h)
	    @test iszero(h[1])

	    j = 1
	    h = WignerDMatrices.Jy_matrix_zbasis!(j,A)
	    @test size(h) == (2j+1,2j+1)
	    @test ishermitian(h)
	    @test all(iszero,diag(h))
	    @test h[1,2] == im*WignerDMatrices.X(j,-j+1)/2
	    @test h[2,3] == im*WignerDMatrices.X(j,-j+2)/2

        j = 2
        h = WignerDMatrices.Jy_matrix_zbasis!(j,A)
	    @test size(h) == (2j+1,2j+1)
	    @test ishermitian(h)
	    @test all(iszero,diag(h))

        # check that this runs
        v = WignerDMatrices.Jy_eigen_nomatch!(1, A)
    end
end
@testset "WignerdMatrix" begin
	@testset "Constructor" begin
		β = π/3
		@testset "integer j" begin
			j = 2
			jh = HalfInt(j)
			d = WignerdMatrix(j,β)
			@test WignerDMatrices.sphericaldegree(d) == j
			@test d.beta == β
			@test WignerDMatrices.eulerangles(d) == (zero(β),β,zero(β))
			@test axes(d) == (-jh:jh,-jh:jh)
			@test axes(d,1) == axes(d,2) == -jh:jh
			@test size(d) == (Integer(2j+1),Integer(2j+1))
			@test size(d,1) == size(d,2) == Integer(2j+1)

			@test d == d

			for i in eachindex(d.dj)
				@test d[i] == d.dj[i]
			end

			@test SpinMatrix(d) == SpinMatrix(collect(d))
		end
		@testset "half-integer j" begin
			j = 3//2
			jh = HalfInt(j)
			d = WignerdMatrix(j,β)
			@test WignerDMatrices.sphericaldegree(d) == j
			@test d.beta == β
			@test WignerDMatrices.eulerangles(d) == (zero(β),β,zero(β))
			@test axes(d) == (-jh:jh,-jh:jh)
			@test axes(d,1) == axes(d,2) == -jh:jh
			@test size(d) == (Integer(2j+1),Integer(2j+1))
			@test size(d,1) == size(d,2) == Integer(2j+1)

			for i in eachindex(d.dj)
				@test d[i] == d.dj[i]
			end

			@test SpinMatrix(d) == SpinMatrix(collect(d))
		end
	end
	function checkdjvalues(test,j,n)
		for β in LinRange(0, 4π, n)
			d = WignerdMatrix(j,β)
			test(d,β)
		end

		@testset "Piby2Radians" begin
			β = Piby2Radians()
			d = WignerdMatrix(j,β)
			test(d,β)
		end
		@testset "ZeroRadians" begin
			β = ZeroRadians()
			d = WignerdMatrix(j,β)
			test(d,β)
		end
		@testset "TwoNPiRadians" begin
		    for n in -2:2
		   		β = TwoNPiRadians(n)
				d = WignerdMatrix(j,β)
				test(d,β)
			end
		end
		@testset "TwoNPlusOnePiRadians" begin
		    for n in -2:2
		   		β = TwoNPlusOnePiRadians(n)
				d = WignerdMatrix(j,β)
				test(d,β)
			end
		end
	end
	@testset "d1//2_mn(β)" begin
		function test(d,β)
			@test isapprox(d[-1/2,-1/2],cos(β/2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[1/2,-1/2],-sin(β/2),atol=1e-14,rtol=1e-8)

			@test isapprox(d[-1/2,1/2],sin(β/2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[1/2,1/2],cos(β/2),atol=1e-14,rtol=1e-8)
		end

		n = 100
		j = 1//2
		checkdjvalues(test,j,n)
	end
    @testset "d1_mn(β)" begin
		function test(d,β)
			@test isapprox(d[1,1],(1+cos(β))/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[1,0],-sin(β)/√2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[1,-1],(1-cos(β))/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[0,1],sin(β)/√2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[0,0],cos(β),atol=1e-14,rtol=1e-8)
			@test isapprox(d[0,-1],-sin(β)/√2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[-1,1],(1-cos(β))/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1,0],sin(β)/√2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1,-1],(1+cos(β))/2,atol=1e-14,rtol=1e-8)
		end

		n = 100
		j = 1
		checkdjvalues(test,j,n)
	end
	@testset "d3//2_mn(β)" begin
	    function test(d,β)
			@test isapprox(d[-3/2,-3/2],cos(β/2)^3,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1/2,-3/2],-√3*sin(β/2)*cos(β/2)^2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[1/2,-3/2],√3*sin(β/2)^2*cos(β/2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[3/2,-3/2],-sin(β/2)^3,atol=1e-14,rtol=1e-8)

			@test isapprox(d[-3/2,-1/2],√3*sin(β/2)*cos(β/2)^2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1/2,-1/2],cos(β/2)*(3cos(β/2)^2-2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[1/2,-1/2],sin(β/2)*(3sin(β/2)^2-2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[3/2,-1/2],√3*sin(β/2)^2*cos(β/2),atol=1e-14,rtol=1e-8)

			@test isapprox(d[-3/2,1/2],√3*sin(β/2)^2*cos(β/2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1/2,1/2],-sin(β/2)*(3sin(β/2)^2-2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[1/2,1/2],cos(β/2)*(3cos(β/2)^2-2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[3/2,1/2],-√3*sin(β/2)*cos(β/2)^2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[-3/2,3/2],sin(β/2)^3,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1/2,3/2],√3*sin(β/2)^2*cos(β/2),atol=1e-14,rtol=1e-8)
			@test isapprox(d[1/2,3/2],√3*sin(β/2)*cos(β/2)^2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[3/2,3/2],cos(β/2)^3,atol=1e-14,rtol=1e-8)
		end

		n = 100
		j = 3//2
		checkdjvalues(test,j,n)
	end
	@testset "d2_mn(β)" begin
		function test(d::WignerdMatrix,β)
			@test isapprox(d[2,1],-sin(β)*(1+cos(β))/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[2,0],1/2*√(3/2)*sin(β)^2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[2,-1],-sin(β)*(1-cos(β))/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[1,1],(2cos(β)^2+cos(β)-1)/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[1,0],-√(3/2)*sin(β)*cos(β),atol=1e-14,rtol=1e-8)
			@test isapprox(d[1,-1],-(2cos(β)^2-cos(β)-1)/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[0,1],√(3/2)*sin(β)*cos(β),atol=1e-14,rtol=1e-8)
			@test isapprox(d[0,0],1/2*(3cos(β)^2-1),atol=1e-14,rtol=1e-8)
			@test isapprox(d[0,-1],-√(3/2)*sin(β)*cos(β),atol=1e-14,rtol=1e-8)

			@test isapprox(d[-1,1],-(2cos(β)^2-cos(β)-1)/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1,0],√(3/2)*sin(β)*cos(β),atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1,-1],(2cos(β)^2+cos(β)-1)/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[-2,1],sin(β)*(1-cos(β))/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-2,0],1/2*√(3/2)*sin(β)^2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-2,-1],sin(β)*(1+cos(β))/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[2,2],(1+cos(β))^2/4,atol=1e-14,rtol=1e-8)
			@test isapprox(d[2,-2],(1-cos(β))^2/4,atol=1e-14,rtol=1e-8)

			@test isapprox(d[1,2],sin(β)*(1+cos(β))/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[1,-2],-sin(β)*(1-cos(β))/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[0,2],1/2*√(3/2)*sin(β)^2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[0,-2],1/2*√(3/2)*sin(β)^2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[-1,2],sin(β)*(1-cos(β))/2,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-1,-2],-sin(β)*(1+cos(β))/2,atol=1e-14,rtol=1e-8)

			@test isapprox(d[-2,2],(1-cos(β))^2/4,atol=1e-14,rtol=1e-8)
			@test isapprox(d[-2,-2],(1+cos(β))^2/4,atol=1e-14,rtol=1e-8)
		end

	    n = 100
		j = 2
		checkdjvalues(test,j,n)
	end
	@testset "linear indexing" begin
	    d1by2 = WignerdMatrix(0.5, rand())
	    d1 = WignerdMatrix(1, rand())

	    for d in [d1by2, d1]
		    for (ind,I) in enumerate(eachindex(d))
		    	@test d[ind] == d[I]
		    end
		end
	end
	@testset "shift angles" begin
		@testset "β + 4π*n" begin
			for j in 0:1//2:2, β in -4π:π/4:4π
		    	d1 = WignerdMatrix(j, β)
		    	d2 = WignerdMatrix(j, β + 4π)
		    	d3 = WignerdMatrix(j, β - 4π)
		    	@test collect(d1) ≈ collect(d2)
		    	@test collect(d1) ≈ collect(d3)
		    end
		end
		@testset "β + (2n+1)*2π" begin
		    for j in 0:1//2:2, β in -4π:π/2:4π
		    	d1 = WignerdMatrix(j, β)
		    	d2 = WignerdMatrix(j, β + 2π)
		    	d3 = WignerdMatrix(j, β - 2π)
		    	d4 = WignerdMatrix(j, β + 6π)
		    	d5 = WignerdMatrix(j, β - 6π)
		    	@test collect(d1) ≈ (-1)^2j .* collect(d2)
		    	@test collect(d1) ≈ (-1)^2j .* collect(d3)
		    	@test collect(d1) ≈ (-1)^2j .* collect(d4)
		    	@test collect(d1) ≈ (-1)^2j .* collect(d5)
		    end
		end
		@testset "negative β" begin
		    for j in 0:1//2:2, β in -4π:π/4:4π
		    	d1 = WignerdMatrix(j, β)
		    	d2 = WignerdMatrix(j, -β)
		    	for m in axes(d1,1), n in axes(d1,2)
		    		@test isapprox(d1[m,n],d2[n,m],atol=1e14,rtol=1e-8)
		    	end
		    end
		end
		@testset "π - β" begin
		    for j in 0:1//2:2, β in -4π:π/4:4π
		    	d1 = WignerdMatrix(j, β)
		    	d2 = WignerdMatrix(j, π - β)
		    	for m in axes(d1,1), n in axes(d1,2)
		    		@test isapprox(d2[m,n],(-1)^(j-n) * d1[-m,n],atol=1e-14,rtol=1e-8)
		    	end
		    end
		end
		@testset "β + (2n+1)π" begin
		    for j in 0:1//2:2, β in -4π:π/4:-4π, k = -1:1
		    	d1 = WignerdMatrix(j, β)
		    	d2 = WignerdMatrix(j, β + (2k+1)π)
		    	d3 = WignerdMatrix(j, β - (2k+1)π)
		    	for m in axes(d1,1), n in axes(d1,2)
		    		@test isapprox(d2[m,n],(-1)^((2k+1)j-n) * d1[m,-n],atol=1e-14,rtol=1e-8)
		    		@test isapprox(d3[m,n],(-1)^(-(2k+1)j-n) * d1[m,-n],atol=1e-14,rtol=1e-8)
		    	end
		    end
		end
	end
	@testset "similar" begin
	    d = WignerdMatrix{Float32}(1, π/3)
	    d′ = similar(d)
	    @test d′ isa WignerdMatrix{Float32}
	    @test sphericaldegree(d′) == sphericaldegree(d)
	    @test d′.beta == d.beta
	    @test typeof(d′.dj) == typeof(d.dj)
	    @test size(d′.dj) == size(d.dj)

	    d′ = similar(d, Float64)
	    @test d′ isa WignerdMatrix{Float64}
	    @test sphericaldegree(d′) == sphericaldegree(d)
	    @test d′.beta == d.beta
	    @test eltype(d′) == Float64
	    @test size(d′.dj) == size(d.dj)
	end
	@testset "LinearAlgebra" begin
		@testset "product" begin
			β = π/3
			for j in [1/2,1,3/2,2]
				d = WignerdMatrix(j,β)
				@test collect(d*d) ≈ collect(d)*collect(d)
			end
		end
		@testset "symmetry" begin
			for j = [1/2, 1]
				d = WignerdMatrix(j ,π/3)
			    @test !issymmetric(d)
				@test !ishermitian(d)
			end
		end
		@testset "det" begin
			for j = [1/2, 1]
				d = WignerdMatrix(j ,π/3)
		    	@test det(d) == 1
		    end
		end
		@testset "trace" begin
			for β in LinRange(-4π, 4π, 100)
			    d = WignerdMatrix(1/2,β)
			    @test isapprox(tr(d), 2cos(β/2), atol=1e-14, rtol=1e-8)
			    @test isapprox(tr(d),sum(d[m,m] for m=-1//2:1//2),atol=1e-14,rtol=1e-8)

			    d = WignerdMatrix(1,β)
			    @test isapprox(tr(d) , 1 + 2cos(β), atol=1e-14, rtol=1e-8)
			    @test isapprox(tr(d),sum(d[m,m] for m=-1:1),atol=1e-14,rtol=1e-8)
			    @test isapprox(tr(d),tr(collect(d)),atol=1e-14,rtol=1e-8)

			    d = WignerdMatrix(3/2,β)
			    @test isapprox(tr(d) , 2(cos(β/2) + cos(3β/2)), atol=1e-14, rtol=1e-8)
			    @test isapprox(tr(d),sum(d[m,m] for m=-3//2:3//2),atol=1e-14,rtol=1e-8)
			    @test isapprox(tr(d),tr(collect(d)),atol=1e-14,rtol=1e-8)

			    d = WignerdMatrix(2,β)
			    @test isapprox(tr(d),1 + 2cos(β) + 2cos(2β),atol=1e-14,rtol=1e-8)
			    @test isapprox(tr(d),sum(d[m,m] for m=-2:2),atol=1e-14,rtol=1e-8)
			    @test isapprox(tr(d), tr(collect(d)),atol=1e-14,rtol=1e-8)
			end

			for n = 0:10
			    d1 = WignerdMatrix(n/2, ZeroRadians())
			    d2 = WignerdMatrix(n/2, TwoNPiRadians(0))
			    d3 = WignerdMatrix(n/2, TwoNPiRadians(2))
			    @test tr(d1) == tr(d2) == tr(d3) == Float64(n + 1)

				d1 = WignerdMatrix(n/2, π)
				d2 = WignerdMatrix(n/2, TwoNPlusOnePiRadians(0))
				rexp = isinteger(n/2) ? (-1)^Int(n/2) : 0
			    @test tr(d1) == tr(d2) == rexp

			    d = WignerdMatrix(n/2, TwoNPiRadians(1))
			    rexp = isinteger(n/2) ? n+1 : -(n+1)
			    @test tr(d) == rexp
			end
		end
		@testset "inv" begin
			rotangles = -4π:π/4:4π
			js = 0:1//2:4
			for β = rotangles, j in js
				d = WignerdMatrix(j,β)
				dinv = inv(d)
				@test dinv isa WignerdMatrix
				@test collect(dinv) ≈ inv(collect(d))
				@test collect(dinv * d) ≈ Diagonal(ones(size(d,1)))
			end

			for j in js
				d = WignerdMatrix(j, ZeroRadians())
				@test inv(d) == d
				@test inv(d).beta === ZeroRadians()

				d = WignerdMatrix(j, TwoNPlusOnePiRadians(0))
				dinv = WignerdMatrix(j, -TwoNPlusOnePiRadians(0))
				@test inv(d) == dinv
				@test inv(d).beta === -TwoNPlusOnePiRadians(0)

				d = WignerdMatrix(j, TwoNPiRadians(1))
				dinv = WignerdMatrix(j, -TwoNPiRadians(1))
				@test inv(d) == dinv
				@test inv(d).beta === -TwoNPiRadians(1)

				d = WignerdMatrix(j, TwoNPlusOnePiRadians(1))
				dinv = WignerdMatrix(j, -TwoNPlusOnePiRadians(1))
				@test inv(d) == dinv
				@test inv(d).beta === -TwoNPlusOnePiRadians(1)
			end
		end
		@testset "diag" begin
			function testdiag(d)
				j = sphericaldegree(d)
				@test begin
					res = diag(d) == [d[m,m] for m=-j:j] == diag(collect(d))
					if !res
						@show j, d.beta
					end
					res
				end
			end
			for j = 0:1//2:10
		    	for β in LinRange(-4π, 4π, 100)
			    	d = WignerdMatrix(j, β)
			    	testdiag(d)
			    end

			    for T in [ZeroRadians, Piby2Radians]
			    	d = WignerdMatrix(j, T())
			    	testdiag(d)
			    end

			    d = WignerdMatrix(j, π)
			    testdiag(d)

			    for T in [TwoNPiRadians, TwoNPlusOnePiRadians], n=-3:3
			    	d = WignerdMatrix(j, T(n))
			    	testdiag(d)
			    end
		    end
		end
	end
end
@testset "WignerDMatrix" begin
    @testset "Constructor" begin
    	β = π/3
		α,γ = π/6, π/2
		j = 2
		D = WignerDMatrix(j,α,β,γ)

        @test D.alpha == α
        @test D.gamma == γ
        @test D.dj.beta == β
        @test WignerDMatrices.eulerangles(D) == (α,β,γ)
        @test D.dj isa WignerdMatrix
        @test D.dj == WignerdMatrix(j,β)
        @test WignerDMatrices.sphericaldegree(D) == j

        @test SpinMatrix(D) == SpinMatrix(collect(D))

        D2 = WignerDMatrix(j,α,β,γ)
        @test D == D2

        D3 = WignerDMatrix(j,0,β,0)
        d = WignerdMatrix(j,β)
        @test D3 == d
        @test d == D3

        D4 = WignerDMatrix(α, d, γ)
        @test D4.dj === d

        @testset "special angles" begin
        	β = π/3
			α,γ = π/6, π/2
			j = 2
            D = WignerDMatrix(j,TwoNPlusOnePiRadians(0),β,γ)
            @test D[1,1] == -D.dj[1,1] * cis(-γ)
            D = WignerDMatrix(j,TwoNPlusOnePiRadians(0),β,ZeroRadians())
            @test D[1,1] == -D.dj[1,1]
            D = WignerDMatrix(j,α,β,TwoNPlusOnePiRadians(0))
            @test D[1,1] == -D.dj[1,1] * cis(-α)
        end
    end
    @testset "indexing" begin
    	β = π/3
		α,γ = π/6, π/2
		j = 2
		D = WignerDMatrix(j,α,β,γ)
        @test D[1,1] ≈ D.dj[1,1]*cis(-(D.alpha + D.gamma))

        # Cartesian indexing
        for n in axes(D,2), m in axes(D,1)
        	Dmn = WignerDMatrices.wignerDmatrixelement(j, m, n, (α,β,γ))
        	@test isapprox(D[m,n], Dmn, atol=1e-13, rtol=1e-8)
        end

        for (i,I) in enumerate(eachindex(D))
        	@test D[i] == D[I]
        end
    end
    @testset "similar" begin
	    D = WignerDMatrix{ComplexF64}(1, 0, π/3, π/3)
	    D′ = similar(D)
	    @test D′ isa WignerDMatrix{ComplexF64}
	    @test sphericaldegree(D′.dj) == sphericaldegree(D.dj)
	    @test D′.dj.beta == D.dj.beta
	    @test typeof(D′.dj) == typeof(D.dj)
	    @test size(D′.dj) == size(D.dj)

	    D′ = similar(D, ComplexF32)
	    @test D′ isa WignerDMatrix{ComplexF32}
	    @test sphericaldegree(D′.dj) == sphericaldegree(D.dj)
	    @test typeof(D′.dj) == typeof(D.dj)
	    @test size(D′.dj) == size(D.dj)
	end
	@testset "LinearAlgebra" begin
		rotangles = -4π:π/3:4π
		@testset "product" begin
			β = π/3
			α,γ = π/6, π/2
			for j in [1/2,1,3/2,2]
				D = WignerDMatrix(j,α,β,γ)
				@test collect(D*D) == collect(D)*collect(D)
				@test collect(D*D.dj) ≈ collect(D)*collect(D.dj)
				@test collect(D.dj*D) ≈ collect(D.dj)*collect(D)
			end
		end
		@testset "symmetry" begin
			for j = [1/2, 1]
				D = WignerDMatrix(j, π/2, π/3, π/4)
			    @test !issymmetric(D)
				@test !ishermitian(D)
			end
		end
		@testset "det" begin
			for j in 0:1//2:4, α in rotangles,
				β in rotangles, γ in rotangles

				D = WignerDMatrix(j, α, β, γ)
		    	@test det(D) == 1
		    end
		end
		@testset "trace" begin
			for j in 0:1//2:10, α in rotangles,
				β in rotangles, γ in rotangles
			    D = WignerDMatrix(j,α,β,γ)
			    @test begin
			    	res = isapprox(tr(D),sum(D[m,m] for m=-j:j),atol=1e-12,rtol=1e-8)
			    	if !res
			    		@show j, α, β, γ
			    	end
			    	res
			    end
			    @test isapprox(tr(D),tr(collect(D)),atol=1e-12,rtol=1e-8)
			end

			for j in 0:1//2:10, α in rotangles, γ in rotangles
				D1 = WignerDMatrix(j,α,ZeroRadians(),γ)
				D2 = WignerDMatrix(j,α,0,γ)
				D3 = WignerDMatrix(j,α,TwoNPiRadians(0),γ)
			    @test isapprox(tr(D1),sum(D1[m,m] for m=-j:j),atol=1e-12,rtol=1e-8)
			    @test isapprox(tr(D1),tr(collect(D1)),atol=1e-12,rtol=1e-8)
			    @test tr(D1) == tr(D2) == tr(D3)

			    D1 = WignerDMatrix(j,α,TwoNPlusOnePiRadians(0),γ)
			    D2 = WignerDMatrix(j,α,π,γ)
			    D3 = WignerDMatrix(j,α,TwoNPlusOnePiRadians(0),γ)
			    @test isapprox(tr(D1),sum(D1[m,m] for m=-j:j),atol=1e-12,rtol=1e-8)
			    @test isapprox(tr(D1),tr(collect(D1)),atol=1e-12,rtol=1e-8)
			    @test tr(D1) == tr(D2) == tr(D3)
			end
		end
		@testset "inv" begin
			for j in 0:1//2:4, α in rotangles,
				β in rotangles, γ in rotangles

				D = WignerDMatrix(j,α,β,γ)
				Dinv = inv(D)
				@test Dinv isa WignerDMatrix
				@test collect(Dinv) ≈ inv(collect(D))
			end
		end
		@testset "diag" begin
		    function testdiag(d)
				j = sphericaldegree(d)
				@test begin
					res = diag(d) == [d[m,m] for m=-j:j] == diag(collect(d))
					if !res
						@show j, d.beta
					end
					res
				end
			end
			for j = 0:1//2:10
		    	for α in rotangles, γ in rotangles

		    		for β in rotangles
				    	D = WignerDMatrix(j, α, β, γ)
				    	testdiag(D)
				    end
				    for T in [ZeroRadians, Piby2Radians]
				    	D = WignerDMatrix(j, α, T(), γ)
				    	testdiag(D)
				    end

				    D = WignerDMatrix(j, α, π, γ)
				    testdiag(D)

				    for T in [TwoNPiRadians, TwoNPlusOnePiRadians], n=-3:3
				    	D = WignerDMatrix(j, α, T(n), γ)
				    	testdiag(D)
				    end
			    end
		    end
		end
	end
end
@testset "subarray" begin
	j = half(1); α, β, γ = 2π*rand(), 2π*rand(), 2π*rand()

	function testviewindexing(d)
		dv = @view d[0.5,0.5]
        @test dv[] == d[0.5,0.5]

        dv = @view d[:,0.5]
        for i in axes(d,1)
			@test dv[i] == d[i,0.5]
		end

		dv = @view d[-0.5:0.5, 0.5]
        for (ind,i) in enumerate(axes(d,1))
			@test dv[ind] == d[i,0.5]
		end

		dv = @view d[:,:]
		for i in eachindex(d)
			@test dv[i] == d[i]
		end
	end
    @testset "WignerdMatrix" begin
        d = WignerdMatrix(j, β)
        testviewindexing(d)
    end
    @testset "WignerDMatrix" begin
        d = WignerDMatrix(j, α, β, γ)
        testviewindexing(d)
    end
end
@testset "broadcasting" begin
    @testset "wignerdmatrix" begin
        d1by2 = WignerdMatrix(1/2, rand())
        d1 = WignerdMatrix(1, rand())

        for d in [d1by2, d1]
	        d2 = d .+ d
	        for I in eachindex(d)
	        	@test d2[I] == 2d[I]
	        end
	        d2 = d .+ 0
	        for I in eachindex(d)
	        	@test d2[I] == d[I]
	        end
	    end
    end
    @testset "wignerDmatrix" begin
        D1by2 = WignerDMatrix(1/2, rand(3)...)
        D1 = WignerDMatrix(1, rand(3)...)

        for d in [D1by2, D1]
	        d2 = d .+ d
	        for I in eachindex(d)
	        	@test d2[I] == 2d[I]
	        end
	        d2 = d .+ 0
	        for I in eachindex(d)
	        	@test d2[I] == d[I]
	        end
	    end
    end
end
@testset "show" begin
    io = IOBuffer()

    function testshow(io, d)
    	summary(io,d)
	    show(io,MIME"text/plain"(), d)
	    show(io, d)
    end

    for β in [rand(), ZeroRadians(), Piby2Radians(),
    	TwoNPiRadians(1), TwoNPlusOnePiRadians(1)]

	    d = WignerdMatrix(1, β)
	    testshow(io, d)

	    D = WignerDMatrix(0, d, 0)
    	testshow(io, D)

	    d = WignerdMatrix(0.5, β)
	    testshow(io, d)

	    D = WignerDMatrix(0, d, 0)
    	testshow(io, D)
	end

    D = WignerDMatrix(1/2, 0, 0, 0)
	testshow(io, D)

    Dv = @view D[:,:]
    testshow(io, Dv)
end


text:
get_ν(β) = β/(1 - β)
get_β(ν) = ν/(1 + ν)
get_k(p, β) = (β - p*β + p)/(2β - p*β + p - 1)
get_d(setup::Setup) = convert(Integer, (size(setup.m, 1) - 1)/size(setup.m, 2))

ext: jl


text:
struct BitsTypeBool
    a::Bool
end

struct BitsTypeChar
    a::Char
end

struct BitsTypeUInt8
    a::UInt8
end

struct BitsTypeUInt16
    a::UInt16
end

struct BitsTypeUInt32
    a::UInt32
end

struct BitsTypeUInt64
    a::UInt64
end

struct BitsTypeUInt
    a::UInt
end

struct BitsTypeInt8
    a::Int8
end

struct BitsTypeInt16
    a::Int16
end

struct BitsTypeInt32
    a::Int32
end

struct BitsTypeInt64
    a::Int64
end

struct BitsTypeInt
    a::Int
end

struct BitsTypeFloat32
    a::Float32
end

struct BitsTypeFloat64
    a::Float64
end

@testset "Single-field bits types" begin
    @test begin
        b = JlrsReflect.reflect([BitsTypeBool])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeBool] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeBool")]
        pub struct BitsTypeBool {
            pub a: ::jlrs::wrappers::inline::bool::Bool,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeChar])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeChar] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeChar")]
        pub struct BitsTypeChar {
            pub a: ::jlrs::wrappers::inline::char::Char,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeUInt8])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeUInt8] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeUInt8")]
        pub struct BitsTypeUInt8 {
            pub a: u8,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeUInt16])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeUInt16] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeUInt16")]
        pub struct BitsTypeUInt16 {
            pub a: u16,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeUInt32])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeUInt32] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeUInt32")]
        pub struct BitsTypeUInt32 {
            pub a: u32,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeUInt64])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeUInt64] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeUInt64")]
        pub struct BitsTypeUInt64 {
            pub a: u64,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeUInt])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeUInt] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeUInt")]
        pub struct BitsTypeUInt {
            pub a: u64,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeInt8])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeInt8] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeInt8")]
        pub struct BitsTypeInt8 {
            pub a: i8,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeInt16])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeInt16] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeInt16")]
        pub struct BitsTypeInt16 {
            pub a: i16,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeInt32])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeInt32] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeInt32")]
        pub struct BitsTypeInt32 {
            pub a: i32,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeInt64])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeInt64] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeInt64")]
        pub struct BitsTypeInt64 {
            pub a: i64,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeInt])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeInt] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeInt")]
        pub struct BitsTypeInt {
            pub a: i64,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeFloat32])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeFloat32] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeFloat32")]
        pub struct BitsTypeFloat32 {
            pub a: f32,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([BitsTypeFloat64])
        sb = JlrsReflect.StringWrappers(b)

        sb[BitsTypeFloat64] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck, IntoJulia)]
        #[jlrs(julia_type = "Main.BitsTypeFloat64")]
        pub struct BitsTypeFloat64 {
            pub a: f64,
        }"""
    end
end


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Implementation of Pigments and BRDFs

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

"""
    Pigment

This abstract type represents a pigment, i.e., a function that associates a color with
each point on a parametric surface ``(u,v)``.

Each subtype of this type must be a callable like `(p::Pigment)(uv::Vec2D)` and must return
the color of the surface as a `RGB{Float32}` in a given [`Vec2D`](@ref) point.

See also: [`UniformPigment`](@ref), [`CheckeredPigment`](@ref), [`ImagePigment`](@ref)
"""
abstract type Pigment end

"""
    (p::Pigment)(uv::Vec2D)

Return the color of the surface in the given point [`Vec2D`](@ref).
"""
function (p::Pigment)(uv::Vec2D)
    p(uv...)
end

#################
# UniformPigment

"""
    UniformPigment <: Pigment

A uniform [`Pigment`](@ref) over the whole surface.
"""
Base.@kwdef struct UniformPigment <: Pigment
    color::RGB{Float32} = WHITE
end

@doc """
    UniformPigment(color::RGB{Float32})

Constructor for a [`UniformPigment`](@ref) instance.
""" UniformPigment(::RGB{Float32})

@doc """
    UniformPigment(; color::RGB{Float32} = WHITE)

Constructor for a [`UniformPigment`](@ref) instance.
""" UniformPigment(; ::RGB{Float32})

"""
    (up::UniformPigment)(u::Float32, v::Float32)

Return the color of the surface in the given point ``(u,v)``.
"""
(up::UniformPigment)(::Float32, ::Float32) = up.color

###################
# CheckeredPigment

"""
    CheckeredPigment{N} <: Pigment

A checkered [`Pigment`](@ref). The number of rows/columns in the checkered pattern is tunable with `N`,
but you cannot have a different number of repetitions along the u/v directions.
"""
struct CheckeredPigment{N} <: Pigment
    color_on::RGB{Float32}
    color_off::RGB{Float32}

    function CheckeredPigment{N}(color_on::RGB{Float32}, color_off::RGB{Float32}) where {N}
        @assert isa(N, Integer)
        new{N}(color_on, color_off)
    end
end

@doc """
    CheckeredPigment{N}(color_on::RGB{Float32}, color_off::RGB{Float32}) where {N}

Constructor for a [`CheckeredPigment`](@ref) instance.
""" CheckeredPigment{N}(::RGB{Float32}, ::RGB{Float32}) where {N}

"""
    CheckeredPigment{N}(; color_on::RGB{Float32} = WHITE,
                      color_off::RGB{Float32} = BLACK) where {N}

Constructor for a [`CheckeredPigment`](@ref) instance.
"""
function CheckeredPigment{N}(; color_on::RGB{Float32} = WHITE, color_off::RGB{Float32} = BLACK) where {N}
    CheckeredPigment{N}(color_on, color_off)
end

"""
    CheckeredPigment(; N::Int = 2, color_on::RGB{Float32} = WHITE,
                      color_off::RGB{Float32} = BLACK) where {N}

Constructor for a [`CheckeredPigment`](@ref) instance.
"""
function CheckeredPigment(; N::Int = 2, color_on::RGB{Float32} = WHITE, color_off::RGB{Float32} = BLACK)
    CheckeredPigment{N}(color_on, color_off)
end

"""
    (cp::CheckeredPigment{N})(u::Float32, v::Float32) where {N}

Return the color of the surface in the given point ``(u,v)``.
"""
function (cp::CheckeredPigment{N})(u::Float32, v::Float32) where {N}
    ceil(u*N)%2 == ceil(v*N)%2 ? cp.color_on : cp.color_off
end

###############
# ImagePigment

"""
    ImagePigment <: Pigment

A textured [`Pigment`](@ref). The texture is given through a PFM image.
"""
struct ImagePigment <: Pigment
    image::HdrImage
end

@doc """
    ImagePigment(image::HdrImage)

Constructor for a [`ImagePigment`](@ref) instance.
""" ImagePigment(::HdrImage)

"""
    ImagePigment(; image::HdrImage = HdrImage(1, 1))

Constructor for a [`ImagePigment`](@ref) instance.
"""
function ImagePigment(; image::HdrImage = HdrImage(1, 1))
    ImagePigment(image)
end

"""
    (ip::ImagePigment)(u::Float32, v::Float32)

Return the color of the surface in the given point ``(u,v)``.
"""
function (ip::ImagePigment)(u::Float32, v::Float32)
    size_row, size_col = size(ip.image)
    f(x) = iszero(x) ? x + 1 : x
    row = u * size_row |> ceil |> Int |> f
    col = v * size_col |> ceil |> Int |> f
    ip.image[row, col]
end

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

"""
    BRDF

An abstract type representing a Bidirectional Reflectance Distribution Function.

Each subtype of this type must include a field `pigment::`[`Pigment`](@ref) storing the pigment on which the BRDF operates.
Each subtype of this type must implement an `at(::NewBRDF, ::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D)` function, where `NewBRDF` should be swubstituted with your new type name. This function evaluates the BRDF of a point with given normal, input and output directions and uv coordinates (which are used to evaluate)

See also: [`DiffuseBRDF`](@ref), [`SpecularBRDF`](@ref),
"""
abstract type BRDF end

##############
# DiffuseBRDF

"""
    DiffuseBRDF <: BRDF

A class representing an ideal diffuse [`BRDF`](@ref) (also called "Lambertian").
"""
Base.@kwdef struct DiffuseBRDF <: BRDF
    pigment::Pigment = UniformPigment()
end

@doc """
    DiffuseBRDF(pigment::Pigment, reflectance::Float32)

Constructor for a [`DiffuseBRDF`](@ref) instance.
""" DiffuseBRDF(::Pigment, ::Float32)

@doc """
    DiffuseBRDF(; pigment::Pigment = UniformPigment(),
                  reflectance::Float32 = 1f0)

Constructor for a [`DiffuseBRDF`](@ref) instance.
""" DiffuseBRDF(; ::Pigment, ::Float32)

"""
    at(brdf::DiffuseBRDF, normal::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D)

Get the radiance, given a point `uv` ([`Vec2D`](@ref)) on the surface with a [`DiffuseBRDF`](@ref)., an incoming
direction `in_dir` and outcoming direction ([`Vec`](@ref)), a `normal` of the surface point ([`Normal`](@ref)).
"""
function at(brdf::DiffuseBRDF, #=normal=#::Normal, #=in_dir=#::Vec, #=out_dir=#::Vec, uv::Vec2D)
    brdf.pigment(uv) * (1f0 / π)
end

"""
    scatter_ray(::DiffuseBRDF, pcg::PCG, incoming_dir::Vec, interaction_point::Point, normal::Normal, depth::Int)

Scatter a ray on the surface.
"""
function scatter_ray(::DiffuseBRDF, pcg::PCG, incoming_dir::Vec,
                     interaction_point::Point, normal::Normal, depth::Int)
    e1, e2, e3 = create_onb_from_z(normal)
    cos_θ_sq = rand(pcg, Float32)
    cos_θ = sqrt(cos_θ_sq)
    φ = 2f0 * π * rand(pcg, Float32)

    Ray(
        interaction_point,
        e1 * cos(φ) * cos_θ + e2 * sin(φ) * cos_θ + e3 * sqrt(1f0 - cos_θ_sq),
        1.0f-3,
        Inf32,
        depth
    )
end

###############
# SpecularBRDF

"""
    SpecularBRDF <: BRDF

A class representing an ideal mirror [`BRDF`](@ref).
"""
Base.@kwdef struct SpecularBRDF <: BRDF
    pigment::Pigment = UniformPigment()
    threshold_angle_rad::Float32 = π / 1800f0
end

@doc """
    SpecularBRDF(pigment::Pigment, threshold_angle_rad::Float32)

Constructor for a [`SpecularBRDF`](@ref) instance.
""" SpecularBRDF(::Pigment, ::Float32)

@doc """
    SpecularBRDF(; pigment::Pigment = UniformPigment(),
                   threshold_angle_rad::Float32 = π / 1800f0)

Constructor for a [`SpecularBRDF`](@ref) instance.
""" SpecularBRDF(; ::Pigment, ::Float32)

"""
    at(brdf::SpecularBRDF, normal::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D)

Get the radiance, given a point `uv` ([`Vec2D`](@ref)) on the surface with a [`SpecularBRDF`](@ref)., an incoming
direction `in_dir` and outcoming direction ([`Vec`](@ref)), a `normal` of the surface point ([`Normal`](@ref)).
"""
function at(brdf::SpecularBRDF, normal::Normal, in_dir::Vec, out_dir::Vec, uv::Vec2D)
    θ_in = normalized_dot(normal, in_dir) |> acos
    θ_out = normalized_dot(normal, out_dir) |> acos
    abs(θ_in - θ_out) < brdf.threshold_angle_rad < brdf.threshold_angle_rad ? brdf.pigment(uv) : BLACK
end

"""
    scatter_ray(::SpecularBRDF, pcg::PCG, incoming_dir::Vec, interaction_point::Point, normal::Normal, depth::Int)

Scatter a ray on the surface.
"""
function scatter_ray(::SpecularBRDF, pcg::PCG, incoming_dir::Vec,
                     interaction_point::Point, normal::Normal, depth::Int)
    ray_dir = normalize(incoming_dir)
    normal = normalize(normal)

    Ray(interaction_point,
        ray_dir - normal * 2f0 * (normal ⋅ ray_dir),
        1f-3,
        Inf32,
        depth)
end

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

"""
    Material

A material with a `brdf::BRDF` ([`BRDF`](@ref)) and and `emitted_radiance::Pigment` ([`Pigment`](@ref)).
"""
Base.@kwdef struct Material
    brdf::BRDF = DiffuseBRDF()
    emitted_radiance::Pigment = UniformPigment(BLACK)
end

@doc """
    Material(brdf::BRDF, emitted_radiance::Pigment)

Constructor for a [`Material`](@ref) instance.
""" Material(::BRDF, ::Pigment)

@doc """
    Material(; brdf::BRDF = DiffuseBRDF(), emitted_radiance::Pigment = UniformPigment(BLACK))

Constructor for a [`Material`](@ref) instance.
""" Material(; ::BRDF, ::Pigment)


text:
"""
    DeferrableLoad

Deferrable Load
"""
struct DeferrableLoad <: Resource
    index::Int  # Unique index of device
    num_timesteps::Int  # Number of time-steps in horizon

    ncycles::Int  # Number of cycles
    cycle_begin::Vector{Int}  # Beginning of each cycle
    cycle_end::Vector{Int}    # End of each cycle
    cycle_energy_min::Vector{Float64}  # Minimal energy consumption for each cycle
    cycle_energy_max::Vector{Float64}  # Maximal energy consumption for each cycle

    pwr_min::Vector{Float64}  # Minimum power consumption at each time
    pwr_max::Vector{Float64}  # Maximum power consumption at each time
    binflag::Bool  # Indicates whether device is hard or soft on-off

    function DeferrableLoad(;
        index::Integer=0,
        T::Integer=0,
        ncycles::Integer=0,
        cycle_begin::Vector{Int}=Int[],
        cycle_end::Vector{Int}=Int[],
        cycle_energy_min::Vector{Float64}=Float64[],
        cycle_energy_max::Vector{Float64}=Float64[],
        pwr_min::Vector{Float64}=Float64[],
        pwr_max::Vector{Float64}=Float64[],
        binflag::Bool=true
    )
        # Dimension checks
        T == length(pwr_min) || throw(DimensionMismatch(
            "T=$T but pwr_min has size $(length(pwr_min))"
        ))
        T == length(pwr_max) || throw(DimensionMismatch(
            "T=$T but pwr_max has size $(length(pwr_max))"
        ))
        ncycles == length(cycle_begin) || throw(DimensionMismatch(
            "N=$ncycles but cycle_begin has size $(length(cycle_begin))"
        ))
        ncycles == length(cycle_end) || throw(DimensionMismatch(
            "N=$ncycles but cycle_end has size $(length(cycle_end))"
        ))
        ncycles == length(cycle_energy_min) || throw(DimensionMismatch(
            "N=$ncycles but cycle_energy_max has size $(length(cycle_energy_max))"
        ))
        ncycles == length(cycle_energy_max) || throw(DimensionMismatch(
            "N=$ncycles but cycle_energy_max has size $(length(cycle_energy_max))"
        ))

        return new(
            index, T,
            ncycles, cycle_begin, cycle_end,
            cycle_energy_min, cycle_energy_max,
            pwr_min, pwr_max,
            binflag
        )
    end
end

"""
    addmodel!(l,
        var2idx, obj, varlb, varub, vartypes,
        row2idx, rowlb, rowub, constrI, constrJ, constrV,
        constr_
    )

Add self to existing model.

# Arguments
- `l`
- `var2idx`
- `obj`
- `varlb`
- `vartypes`
- `row2idx`
- `rowlb`
- `rowub`
- `constrI`
- `constrJ`
- `constrV`
- `constr_`
"""
function add_resource_to_model!(
    m::MOI.ModelLike,
    h::House,
    l::DeferrableLoad,
    var2idx, con2idx
)
    T = l.num_timesteps

    # ==========================================
    #    I. Add local variables
    # ==========================================
    pnet = MOI.add_variables(m, T)  # Net load
    defer = MOI.add_variables(m, T)  # On-Off indicator
    # Update list of variable indices
    for t in 1:T
        var2idx[(:curt, l.index, :pnet, t)] = pnet[t]
        var2idx[(:curt, l.index, :uind, t)] = defer[t]
    end

    # Variable bounds
    for t in 1:T
        MOI.add_constraint(m,
            MOI.SingleVariable(pnet[t]),
            MOI.Interval(l.pwr_min[t], l.pwr_max[t])
        )  # min-max load

        if l.binflag
            # Add binary constraint
            MOI.add_constraint(m,
                MOI.SingleVariable(defer[t]),
                MOI.ZeroOne()
            )
        else
            # Just add bounds
            MOI.add_constraint(m,
                MOI.SingleVariable(defer[t]),
                MOI.Interval(0.0, 1.0)
            )
        end
    end

    # ==========================================
    #    II. Update linking constraints
    # ==========================================
    for t in 1:T
        cidx = con2idx[(:house, h.index, :link, t)]
        MOI.modify(m, cidx, MOI.ScalarCoefficientChange(pnet[t], -1.0))
    end

    # ==========================================
    #    III. Add local constraints
    # ==========================================
    # Energy consumption for each cycle
    # E_min[k] <= p[t_1] + ... + p[t_k] <= E_max[k]
    # where t_1, t_k are the first and last period of cycle k
    for k in 1:l.ncycles
        # create constraint for cycle k
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction{Float64}([
                MOI.ScalarAffineTerm(1.0, pnet[t])
                for t in (l.cycle_begin[k]):(l.cycle_end[k])
            ], 0.0),
            MOI.Interval(l.cycle_energy_min[k], l.cycle_energy_max[k])
        )
    end

    # On-off
    for t in 1:T
        # Pmin[t] * u[t] <= p[t], i.e., p[t] - Pmin[t] * u[t] >= 0
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction{Float64}([
                MOI.ScalarAffineTerm(1.0, pnet[t]),
                MOI.ScalarAffineTerm(-l.pwr_min[t], defer[t])
            ], 0.0),
            MOI.GreaterThan(0.0)
        )
        # p[t] <= Pmax * u[t], i.e., p[t] - Pmax[t] * u[t] <= 0
        MOI.add_constraint(m,
            MOI.ScalarAffineFunction{Float64}([
                MOI.ScalarAffineTerm(1.0, pnet[t]),
                MOI.ScalarAffineTerm(-l.pwr_max[t], defer[t])
            ], 0.0),
            MOI.LessThan(0.0)
        )
    end

    # ==========================================
    #    IV. Add sub-resources to current model
    # ==========================================
    # (None here)

    return nothing
end


text:
export regularize

"derivative of regularizar"
function regularize(W::AFArray{Float32})
  reg=W # L1:sign(W)
end


text:
n = 500
B = rand(n, n)
x = rand(n)

A = x*x'
y = B*x

println(A[1,1])

A = eye(n)
y = A\x

println(sum(abs.(x-y)))


text:
module Grids

using DocStringExtensions
export Grid, Cube

include("directional_funcs.jl")
include("stretch_funcs.jl")
include("stretch_params.jl")
include("stretch_param_match.jl")
include("coordinates.jl")

"""
    Grid

An N-dimensional Cartesian, structured, computational grid.
"""
struct Grid{FT, T, N, CE, CO}
  c::NTuple{N,T}
  cent::CE
  corn::CO
end
function Grid(c::NTuple{N,T}) where {N,T}
    cent = CenterSpace(c)
    corn = VertexSpace(c)
    args = (eltype(c[1].c.h), T, N, typeof.((cent, corn))...)
    Grid{args...}(c, cent, corn)
end

abstract type Dimension end
struct XDim <: Dimension end
struct YDim <: Dimension end
struct ZDim <: Dimension end

struct Space{X,Y,Z}
    x::X
    y::Y
    z::Z
end
CenterSpace(c::Tuple) = Space(
        Space1D{XDim}(c[1].c.h),
        Space1D{YDim}(c[2].c.h),
        Space1D{ZDim}(c[3].c.h)
    )
VertexSpace(c::Tuple) = Space(
        Space1D{XDim}(c[1].n.h),
        Space1D{YDim}(c[2].n.h),
        Space1D{ZDim}(c[3].n.h)
    )

struct Space1D{dim,H}
    h::H
end
Space1D{dim}(h::H) where {H,dim} = Space1D{dim,H}(h)

# Allow for field-like indexing:
Base.@propagate_inbounds Base.getindex(c::Space1D{XDim}, i,j,k) = Base.getindex(c.h, i)
Base.@propagate_inbounds Base.getindex(c::Space1D{YDim}, i,j,k) = Base.getindex(c.h, j)
Base.@propagate_inbounds Base.getindex(c::Space1D{ZDim}, i,j,k) = Base.getindex(c.h, k)

function Cube(a::FT, b::FT, n::Int; kwargs...) where {FT}
    c = Coordinates(a,b,n; kwargs...)
    return Grid((c,c,c))
end

function Base.size(grid::Grid, t::Tuple)
    @assert length(t) == 3
    return (length(grid.c[k], dl) for (k,dl) in enumerate(t))
end

function Base.show(io::IO, g::Grid{FT,T,N}) where {FT,T,N}
  println(io, "---------------------- Grid")
  for i in 1:N
    println(io, g.c[i])
  end
  println(io, "---------------------- end grid")
end

Base.length(grid::Grid) = grid.c[1].n.s*grid.c[2].n.s*grid.c[3].n.s

export all_coords
function all_coords(grid)
    return (;xn = coords(grid.c[1], Vertex1D()).h,
             yn = coords(grid.c[2], Vertex1D()).h,
             zn = coords(grid.c[3], Vertex1D()).h,
             xc = coords(grid.c[1], Center1D()).h,
             yc = coords(grid.c[2], Center1D()).h,
             zc = coords(grid.c[3], Center1D()).h)
end

end


text:
using Revise
using Plots
using DynamicalSystems
using DifferentialEquations
using Basins

# http://dx.doi.org/10.1155/2013/256092
# Multiple Coexisting Attractors and Hysteresis in the Generalized Ueda Oscillator
@inline @inbounds function ueda(u, p, t)
    b = p[1]; F = p[2]; ω = p[3]; α = p[4]
	x=u[1]; y=u[2]; z=u[3];
    dx = y
    dy = -b*y-sign(x)*abs(x)^α+F*sin(z)
    dz = ω
    return SVector{3}(dx, dy, dz)
end

ω=1
α=4
p= [0.05, 7.5, ω, α]
ds = ContinuousDynamicalSystem(ueda, rand(3), p)
integ  = integrator(ds; alg=Tsit5(),  reltol=1e-8, save_everystep=false)

xres=100
yres=100

xg = range(-2.,2.,length=xres)
yg = range(-7.,7.,length=yres)

@time basin = Basins.basins_map2D(xg, yg, integ; T=2π/ω)

plot(xg,yg,basin', seriestype=:heatmap)


text:
"""
Loads dependencies and bootstraps a Genie app. Exposes core Genie functionality.
"""
module Genie

import Revise

push!(LOAD_PATH, @__DIR__)

import Inflector

include("Configuration.jl")
const config = Configuration.Settings(app_env = ENV["GENIE_ENV"])

include("constants.jl")

import Sockets
import Logging

using Reexport

include(joinpath(@__DIR__, "genie_types.jl"))

include("HTTPUtils.jl")
include("Exceptions.jl")
include("App.jl")
include("genie_module.jl")
include("Util.jl")
include("FileTemplates.jl")
include("Toolbox.jl")
include("Generator.jl")
include("Encryption.jl")
include("Cookies.jl")
include("Input.jl")
include("Router.jl")
include("Renderer.jl")
include("WebChannels.jl")
include("WebThreads.jl")
include("Headers.jl")
include("Assets.jl")
include("AppServer.jl")
include("Commands.jl")
include("Responses.jl")
include("Requests.jl")
include("Flash.jl")
include("Plugins.jl")
include("Deploy.jl")

# === #
# EXTRAS #

include("Cache.jl")
config.cache_storage == :File && include("cache_adapters/FileCache.jl")

include("Sessions.jl")

export serve, up, down, loadapp, genie
@reexport using .Router

"""
    serve(path::String = Genie.config.server_document_root, params...; kwparams...)

Serves a folder of static files located at `path`. Allows Genie to be used as a static files web server.
The `params` and `kwparams` arguments are forwarded to `Genie.startup()`.

# Arguments
- `path::String`: the folder of static files to be served by the server
- `params`: additional arguments which are passed to `Genie.startup` to control the web server
- `kwparams`: additional keyword arguments which are passed to `Genie.startup` to control the web server

# Examples
```julia-repl
julia> Genie.serve("public", 8888, async = false, verbose = true)
[ Info: Ready!
2019-08-06 16:39:20:DEBUG:Main: Web Server starting at http://127.0.0.1:8888
[ Info: Listening on: 127.0.0.1:8888
[ Info: Accept (1):  🔗    0↑     0↓    1s 127.0.0.1:8888:8888 ≣16
```
"""
function serve(path::String = Genie.config.server_document_root, params...; kwparams...)
  cd(path)
  path = "."

  Router.route("/") do
    Router.serve_static_file(path, root = path)
  end
  Router.route(".*") do
    Router.serve_static_file(Router.params(:REQUEST).target, root = path)
  end

  up(params...; kwparams...)
end

### NOT EXPORTED ###

"""
    newapp(path::String = "."; autostart::Bool = true, fullstack::Bool = false, dbsupport::Bool = false, mvcsupport::Bool = false) :: Nothing

Scaffolds a new Genie app, setting up the file structure indicated by the various arguments.

# Arguments
- `path::String`: the name of the app and the path where to bootstrap it
- `autostart::Bool`: automatically start the app once the file structure is created
- `fullstack::Bool`: the type of app to be bootstrapped. The fullstack app includes MVC structure, DB connection code, and asset pipeline files.
- `dbsupport::Bool`: bootstrap the files needed for DB connection setup via the SearchLight ORM
- `mvcsupport::Bool`: adds the files used for HTML+Julia view templates rendering and working with resources

# Examples
```julia-repl
julia> Genie.newapp("MyGenieApp")
2019-08-06 16:54:15:INFO:Main: Done! New app created at MyGenieApp
2019-08-06 16:54:15:DEBUG:Main: Changing active directory to MyGenieApp
2019-08-06 16:54:15:DEBUG:Main: Installing app dependencies
 Resolving package versions...
  Updating `~/Dropbox/Projects/GenieTests/MyGenieApp/Project.toml`
  [c43c736e] + Genie v0.10.1
  Updating `~/Dropbox/Projects/GenieTests/MyGenieApp/Manifest.toml`

2019-08-06 16:54:27:INFO:Main: Starting your brand new Genie app - hang tight!
 _____         _
|   __|___ ___|_|___
|  |  | -_|   | | -_|
|_____|___|_|_|_|___|

┌ Info:
│ Starting Genie in >> DEV << mode
└
[ Info: Logging to file at MyGenieApp/log/dev.log
[ Info: Ready!
2019-08-06 16:54:32:DEBUG:Main: Web Server starting at http://127.0.0.1:8000
2019-08-06 16:54:32:DEBUG:Main: Web Server running at http://127.0.0.1:8000
```
"""
const newapp = Generator.newapp

const newapp_webservice = Generator.newapp_webservice
const newapp_mvc = Generator.newapp_mvc
const newapp_fullstack = Generator.newapp_fullstack
const new = Generator.new

"""
    loadapp(path::String = "."; autostart::Bool = false) :: Nothing

Loads an existing Genie app from the file system, within the current Julia REPL session.

# Arguments
- `path::String`: the path to the Genie app on the file system.
- `autostart::Bool`: automatically start the app upon loading it.

# Examples
```julia-repl
shell> tree -L 1
.
├── Manifest.toml
├── Project.toml
├── bin
├── bootstrap.jl
├── config
├── env.jl
├── genie.jl
├── log
├── public
├── routes.jl
└── src

5 directories, 6 files

julia> using Genie

julia> Genie.loadapp(".")
 _____         _
|   __|___ ___|_|___
|  |  | -_|   | | -_|
|_____|___|_|_|_|___|

┌ Info:
│ Starting Genie in >> DEV << mode
└
[ Info: Logging to file at MyGenieApp/log/dev.log
```
"""
function loadapp(path::String = "."; autostart::Bool = false) :: Nothing
  Core.eval(Main, quote
      include(joinpath($path, $(Genie.BOOTSTRAP_FILE_NAME)))
  end)

  autostart && (Core.eval(Main.UserApp, :(up())))

  nothing
end

"""
    startup(port::Int = Genie.config.server_port, host::String = Genie.config.server_host;
        ws_port::Int = Genie.config.websockets_port, async::Bool = ! Genie.config.run_as_server) :: Nothing

Starts the web server. Alias for `AppServer.startup`

# Arguments
- `port::Int`: the port used by the web server
- `host::String`: the host used by the web server
- `ws_port::Int`: the port used by the Web Sockets server
- `async::Bool`: run the web server task asynchronously

# Examples
```julia-repl
julia> startup(8000, "127.0.0.1", async = false)
[ Info: Ready!
Web Server starting at http://127.0.0.1:8000
```
"""
const startup = AppServer.startup
const up = startup

const down = AppServer.down

### PRIVATE ###

"""
    run() :: Nothing

Runs the Genie app by parsing the command line args and invoking the corresponding actions.
Used internally to parse command line arguments.
"""
function run(; server::Union{Sockets.TCPServer,Nothing} = nothing) :: Nothing
  Commands.execute(Genie.config, server = server)

  nothing
end

"""
    genie() :: Nothing

Formerly, this has been `genie.jl` in an app's base directory.
"""
function genie(; context = @__MODULE__) :: Nothing
  haskey(ENV, "GENIE_ENV") || (ENV["GENIE_ENV"] = "dev")

  if !haskey(ENV, "HOST")
    ENV["HOST"] = (ENV["GENIE_ENV"] == "dev") ? "127.0.0.1" : "0.0.0.0"
  end

  ### EARLY BIND TO PORT FOR HOSTS WITH TIMEOUT ###
  EARLYBINDING = if haskey(ENV, "EARLYBIND") && lowercase(ENV["EARLYBIND"]) == "true" && haskey(ENV, "PORT")
    printstyled("\nEarly binding to host $(ENV["HOST"]) and port $(ENV["PORT"]) \n", color = :light_blue, bold = true)
    try
      Sockets.listen(parse(Sockets.IPAddr, ENV["HOST"]), parse(Int, ENV["PORT"]))
    catch ex
      @show ex

      printstyled("\nFailed early binding!\n", color = :red, bold = true)
      nothing
    end
  else
    nothing
  end

  ### OFF WE GO! ###
  ROOT_PATH = pwd()
  push!(LOAD_PATH, ROOT_PATH, "src")

  load(context = context)
  run(server = EARLYBINDING)

  nothing
end

end


text:
using ExprSearch
using DecisionTrees

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
setrotations!(cc::AbstractNsolt, (θ, μ)) = setrotations(cc, θ, μ)

setrotations(cc::AbstractNsolt, args...) = setrotations!(similar(cc), args...)

function getrotations(nsolt::AbstractNsolt{T,D}) where {T,D}
    initpms = getrotations_init(nsolt)
    proppms = map(enumerate(nstages(nsolt))) do (d, nstg)
        map(1:nstg) do k
            getrotations_prop(nsolt, d, k)
        end
    end
    pms = vcat(initpms, proppms...)

    return (vcat(map(t->t[1], pms)...), vcat(map(t->t[2], pms)...))
end

getrotations_init(nsolt::RnsoltTypeI) = [ mat2rotations(nsolt.W0), mat2rotations(nsolt.U0) ]
getrotations_init(nsolt::RnsoltTypeII) = [ mat2rotations(nsolt.W0), mat2rotations(nsolt.U0) ]
getrotations_init(nsolt::CnsoltTypeI) = [ mat2rotations(nsolt.V0) ]
getrotations_init(nsolt::CnsoltTypeII) = [ mat2rotations(nsolt.V0) ]

function getrotations_prop(nsolt::RnsoltTypeI, d, k)
    mat2rotations(nsolt.Udks[d][k])
end
function getrotations_prop(nsolt::RnsoltTypeII, d, k)
    pms = mat2rotations.([ nsolt.Wdks[d][k], nsolt.Udks[d][k] ])
    return (vcat(map(t->t[1], pms)...), vcat(map(t->t[2], pms)...))
end
function getrotations_prop(nsolt::CnsoltTypeI, d, k)
    pms = mat2rotations.([ nsolt.Wdks[d][k], nsolt.Udks[d][k] ])
    return (vcat(map(t->t[1], pms)..., nsolt.θdks[d][k]), vcat(map(t->t[2], pms)...))
end
function getrotations_prop(nsolt::CnsoltTypeII, d, k)
    pms = mat2rotations.([ nsolt.Wdks[d][k], nsolt.Udks[d][k], nsolt.Ŵdks[d][k], nsolt.Ûdks[d][k] ])
    return (vcat(map(t->t[1], pms)..., nsolt.θ1dks[d][k], nsolt.θ2dks[d][k]), vcat(map(t->t[2], pms)...))
end

function setrotations!(nsolt::AbstractNsolt{T,D}, θ::AbstractArray, μ::AbstractArray) where {T,D}
    npiθ, npiμ = nparamsinit(nsolt)
    nppsθ, nppsμ = nparamsperstage(nsolt)
    dlmθ = intervals([ sum(npiθ), (sum(nppsθ) .* nstages(nsolt))... ])
    dlmμ = intervals([ sum(npiμ), (sum(nppsμ) .* nstages(nsolt))... ])
    θinit, μinit = @view(θ[dlmθ[1]]), @view(μ[dlmμ[1]])
    θprop, μprop = map(t->@view(θ[t]), dlmθ[2:end]), map(t->@view(μ[t]), dlmμ[2:end])

    setrotations_init!(nsolt, θinit, μinit)

    foreach(1:D, nstages(nsolt), θprop, μprop) do d, nstg, θd, μd
        dlmθk = intervals(fill(sum(nppsθ), nstg))
        dlmμk = intervals(fill(sum(nppsμ), nstg))
        θdk = map(k->@view(θd[k]), dlmθk)
        μdk = map(k->@view(μd[k]), dlmμk)
        foreach(1:nstg, θdk, μdk) do k, θ, μ
            setrotations_prop!(nsolt, d, k, θ, μ)
        end
    end

    return nsolt
end

function nparamsinit(nsolt::RnsoltTypeI)
    ([ ngivensangles.(nsolt.nChannels)... ], [ nsolt.nChannels... ])
end
function nparamsinit(nsolt::RnsoltTypeII)
    ([ ngivensangles.(nsolt.nChannels)... ], [ nsolt.nChannels... ])
end
function nparamsinit(nsolt::CnsoltTypeI)
    ([ ngivensangles(nsolt.nChannels) ], [ nsolt.nChannels ])
end
function nparamsinit(nsolt::CnsoltTypeII)
    ([ ngivensangles(nsolt.nChannels) ], [ nsolt.nChannels ])
end

function nparamsperstage(nsolt::RnsoltTypeI)
    ([ ngivensangles(nsolt.nChannels[2]) ], [ nsolt.nChannels[2] ])
end
function nparamsperstage(nsolt::RnsoltTypeII)
    ([ ngivensangles.(nsolt.nChannels)... ], [ nsolt.nChannels... ])
end
function nparamsperstage(nsolt::CnsoltTypeI)
    hP = fld(nsolt.nChannels, 2)
    (
        [ fill(ngivensangles(hP), 2)..., fld(nsolt.nChannels, 4) ],
        [ fill(hP, 2)... ],)
end
function nparamsperstage(nsolt::CnsoltTypeII)
    fP, cP = fld(nsolt.nChannels,2), cld(nsolt.nChannels, 2)
    (
        [ fill(ngivensangles(fP),2)..., fill(ngivensangles(cP),2)..., fill(fld(nsolt.nChannels,4),2)... ],
        [ fP, fP, cP, cP] ,)
end

function setrotations_init!(nsolt::RnsoltTypeI, θ::AbstractArray, μ::AbstractArray)
    dlsθ, dlsμ = intervals.(nparamsinit(nsolt))
    nsolt.W0 .= rotations2mat(θ[dlsθ[1]], μ[dlsμ[1]])
    nsolt.U0 .= rotations2mat(θ[dlsθ[2]], μ[dlsμ[2]])
    return nsolt
end

function setrotations_init!(nsolt::RnsoltTypeII, θ::AbstractArray, μ::AbstractArray)
    dlsθ, dlsμ = intervals.(nparamsinit(nsolt))
    nsolt.W0 .= rotations2mat(θ[dlsθ[1]], μ[dlsμ[1]])
    nsolt.U0 .= rotations2mat(θ[dlsθ[2]], μ[dlsμ[2]])
    return nsolt
end

function setrotations_init!(nsolt::CnsoltTypeI, θ::AbstractArray, μ::AbstractArray)
    nsolt.V0 .= rotations2mat(θ, μ)
    return nsolt
end

function setrotations_init!(nsolt::CnsoltTypeII, θ::AbstractArray, μ::AbstractArray)
    nsolt.V0 .= rotations2mat(θ, μ)
    return nsolt
end

function setrotations_prop!(nsolt::RnsoltTypeI, d::Integer, k::Integer, θ::AbstractArray, μ::AbstractArray)
    nsolt.Udks[d][k] .= rotations2mat(θ, μ)
    return nsolt
end

function setrotations_prop!(nsolt::RnsoltTypeII, d::Integer, k::Integer, θ::AbstractArray, μ::AbstractArray)
    dlsθ, dlsμ = intervals.(nparamsperstage(nsolt))
    nsolt.Wdks[d][k] .= rotations2mat(θ[dlsθ[1]], μ[dlsμ[1]])
    nsolt.Udks[d][k] .= rotations2mat(θ[dlsθ[2]], μ[dlsμ[2]])
    return nsolt
end

function setrotations_prop!(nsolt::CnsoltTypeI, d::Integer, k::Integer, θ::AbstractArray, μ::AbstractArray)
    dlsθ, dlsμ = intervals.(nparamsperstage(nsolt))
    nsolt.Wdks[d][k] .= rotations2mat(θ[dlsθ[1]], μ[dlsμ[1]])
    nsolt.Udks[d][k] .= rotations2mat(θ[dlsθ[2]], μ[dlsμ[2]])
    nsolt.θdks[d][k] .= θ[dlsθ[3]]
    return nsolt
end

function setrotations_prop!(nsolt::CnsoltTypeII, d::Integer, k::Integer, θ::AbstractArray, μ::AbstractArray)
    dlsθ, dlsμ = intervals.(nparamsperstage(nsolt))
    nsolt.Wdks[d][k] .= rotations2mat(θ[dlsθ[1]], μ[dlsμ[1]])
    nsolt.Udks[d][k] .= rotations2mat(θ[dlsθ[2]], μ[dlsμ[2]])
    nsolt.Ŵdks[d][k] .= rotations2mat(θ[dlsθ[3]], μ[dlsμ[3]])
    nsolt.Ûdks[d][k] .= rotations2mat(θ[dlsθ[4]], μ[dlsμ[4]])
    nsolt.θ1dks[d][k] .= θ[dlsθ[5]]
    nsolt.θ2dks[d][k] .= θ[dlsθ[6]]
    return nsolt
end

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
using MechanicalUnits
using Test
# Note that the printing functionality has been moved from here to a clone/ fork of Unitful.
# This is now mostly redundant.
shortp(x) = repr(x, context = :color=>true)
longp(x) = repr(:"text/plain", x, context = :color=>true)

global const sInt = typeof(Int(1)) == Int64 ? "Int64" : "Int32"

@testset "Most basic" begin
    testunits = [m ,  s ,  kg ]
    teststrings = ["m" ,  "s" ,  "kg" ]
    for (u, su) in zip(testunits, teststrings)
        st = "2\e[36m$su\e[39m"
        @test shortp(2u) == st
        @test longp(2u) == st
    end
    # Products of most basic quantities
    q1, q2, res = m , s , "6\e[36mm\e[39m∙\e[36ms\e[39m"
    @test shortp(q1 * q2) == longp(q1 * q2)
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = m , kg , "6\e[36mkg\e[39m∙\e[36mm\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = s , m , "6\e[36mm\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = s , kg , "6\e[36mkg\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = kg , m , "6\e[36mkg\e[39m∙\e[36mm\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = kg , s , "6\e[36mkg\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
end

@testset "Special temperature units, Affine units" begin
    testunits = [°C, °F]
    teststrings = ["°C", "°F"]
    for (u, su) in zip(testunits, teststrings)
        st = "3\e[36m$su\e[39m"
        @test shortp(3u) == st
        @test longp(3u) == st
    end
    # These quantities can't be multiplied, so that's it.
end

@testset "Fahrenheit intervals, Kelvin, directly derived units" begin
    testunits = [Ra, K, h, μm, minute]
    teststrings = ["Ra", "K", "h", "μm", "minute"]
    for (u, su) in zip(testunits, teststrings)
        st = "2\e[36m$su\e[39m"
        @test shortp(2u) == st
        @test longp(2u) == st
    end
    q1, q2, res = m , s , "6\e[36mm\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = m , kg , "6\e[36mkg\e[39m∙\e[36mm\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = s , m , "6\e[36mm\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = s , kg , "6\e[36mkg\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = kg , m , "6\e[36mkg\e[39m∙\e[36mm\e[39m"
    @test shortp(2q1 * 3q2) == res
    q1, q2, res = kg , s , "6\e[36mkg\e[39m∙\e[36ms\e[39m"
    @test shortp(2q1 * 3q2) == res
end

@testset "Rational quantities" begin
    shortp((1//100)kg) == "(1//100)\e[36mkg\e[39m"
    longp((1//100)kg) == "(1//100)\e[36mkg\e[39m"
end
@testset "Constructors from type signatures" begin
    # Unitfu now prints shorter type signatures, which
    # on reconstruction makes FreeUnits instead of the original.
    # The shorter output is more readable.
end

@testset "Type signatures, exponents -4 to 4" begin
    dimdi = Dict([m => " ᴸ", s => " ᵀ", kg => " ᴹ",
            Ra => " ᶿ", K => " ᶿ", h => " ᵀ",
            μm => " ᴸ", minute => " ᵀ"])
    for bu in ["m", "s", "kg", "Ra", "K", "h", "μm", "minute"],
        di in ["⁻⁴" , "⁻³", "⁻²" , "⁻¹", "²", "³", "⁴"]
        res = shortp(typeof(eval(Symbol(bu*di))))
        expec1 = "FreeUnits{(\e[36m" * bu * di* "\e[39m,),"
        expec2a = ""
        expec2b = " "
        expec3 = dimdi[eval(Meta.parse(bu)) ]
        expec4 =  di * ", nothing}"
        @test res == expec1*expec2a*expec3*expec4 ||
              res == expec1*expec2b*expec3*expec4
    end
end

@testset "Arrays with units" begin
    a1 = [1 2]m
    st ="[2 4]\e[36mm\e[39m"
    @test shortp(2a1) == st

    st = "1×2 Matrix{Quantity{Int64,  ᴸ, FreeUnits{(\e[36mm\e[39m,),  ᴸ, nothing}}}:\n 2  4"
    @test longp(2a1) == st

    a2 = [1 2]m*s^-1
    st = "[2 4]\e[36mm\e[39m∙\e[36ms⁻¹\e[39m"
    @test shortp(2a2) == st

    a3 = [10, 6, 2, -2]m
    st = "[20, 12, 4, -4]\e[36mm\e[39m"
    @test shortp(2a3) == st

    st = "4-element Vector{Quantity{Int64,  ᴸ, FreeUnits{(\e[36mm\e[39m,),  ᴸ, nothing}}}:\n 20\n 12\n  4\n -4"
    @test longp(2a3) == st
end

@testset "Tuples with units" begin
    a1 = (1, 2)m
    st ="(2, 4)\e[36mm\e[39m"
    @test shortp(2 .*a1) == st
    @test longp(2 .*a1) == st
    a2 = (1, 2)m*s^-1
    st = "(2, 4)\e[36mm\e[39m∙\e[36ms⁻¹\e[39m"
    @test shortp(2 .*a2) == st
    a3 = (10, 6, 2, -2)m
    st = "(10, 6, 2, -2)\e[36mm\e[39m"
    @test shortp(a3) == st
end

@testset "Ranges" begin
    a1 = 1m:2m:5m
    st = "(2:4:10)\e[36mm\e[39m"
    @test shortp(2 .*a1) == st
    @test longp(2 .*a1) == st
    a2 = (5:-2:-1)m*s^-1
    st = "[10, 6, 2, -2]\e[36mm\e[39m∙\e[36ms⁻¹\e[39m"
    @test shortp(collect(2 .*a2)) == st
end

#
@testset "Dimensions" begin
    u  = s*m*kg*K
    @test shortp(u) == "\e[36mkg\e[39m∙\e[36mK\e[39m∙\e[36mm\e[39m∙\e[36ms\e[39m"
    @test shortp(dimension(u)) == " ᴸ∙ ᴹ∙ ᶿ∙ ᵀ"
    st = "FreeUnits{(\e[36mkg\e[39m, \e[36mK\e[39m, \e[36mm\e[39m, \e[36ms\e[39m),  ᴸ∙ ᴹ∙ ᶿ∙ ᵀ, nothing}"
    @test shortp(typeof(u)) == st

    @test shortp(typeof(dimension(u))) == "Dimensions{(Dimension{:Length}(1//1), Dimension{:Mass}(1//1), Dimension{:Temperature}(1//1), Dimension{:Time}(1//1))}"
    @test shortp(dimension(u^2)) == " ᴸ²∙ ᴹ²∙ ᶿ²∙ ᵀ²"
    @import_expand A mol
    v  = A * mol
    @test shortp(v) == "\e[36mA\e[39m∙\e[36mmol\e[39m"
    @test shortp(dimension(v)) == " ᴺ∙ ᴵ"
    st = "FreeUnits{(\e[36mA\e[39m, \e[36mmol\e[39m),  ᴺ∙ ᴵ, nothing}"
    @test shortp(typeof(v)) == st
    @test shortp(typeof(dimension(v))) == "Dimensions{(Dimension{:Amount}(1//1), Dimension{:Current}(1//1))}"
    @test shortp(dimension(v^2)) == " ᴺ²∙ ᴵ²"
end

max_stars_repo_name: hustf/MechanicalUnits

max_issues_repo_name: hustf/MechanicalUnits


text:
#-----------------------------------------------------------------------
#	SAMPLE GENERATION
#-----------------------------------------------------------------------

using LinearAlgebra, Serialization, BenchmarkTools

#include("../src/squirrel.jl")
include("../src/metric.jl")

g  	= metric.g
gk 	= metric.ge

Nsamp   = 1000

nb	= 24
tol	= 1e-14
ξ1	= 1e-18
ξ2	= 1e1
sufx 	= ""
dir 	= "../res/"
Nfs	= string(Nsamp)

#-----------------------------------------------------------------------
#	SAMPLES IN ANALOGUE GEOMETRY WITH ATMOSPHERIC & IONOSPHERIC EFFECTS
#-----------------------------------------------------------------------

tc	= squirrel.seval.gen(Nsamp,g,6)

pfx	= "tct"

tctloc	= dir*pfx*"-"*Nfs*sufx

Serialization.serialize(tctloc,squirrel.seval.tc2tup(tc))

#-----------------------------------------------------------------------
#	SAMPLES IN KERR GEOMETRY
#-----------------------------------------------------------------------

tck 	= squirrel.seval.gen(Nsamp,gk,6)

pfx	= "tck"

tckloc	= dir*pfx*"-"*Nfs*sufx

Serialization.serialize(tckloc,squirrel.seval.tc2tup(tck))

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


text:
module Dykstra
using LinearAlgebra
using SparseArrays
#################### LinearConstraints
struct LinearConstraint{T}
    A::SparseMatrixCSC{T,Int64}
    b::SparseVector{T, Int64}
    P::SparseMatrixCSC{T, Int64}
    q::SparseVector{T, Int64}
    rez::SparseVector{T, Int64}
end

function get_AiAA(A)
    X = A*A'
    if X == Diagonal(X)
        @assert sum(diag(X) .== 0)==0 "A*A' must not have a zero determinant"
        # @info "X is diagonal with non-zero diagonal values ! huge win."
        return A'spdiagm(1 ./ diag(X))
    else
        return A'sparse(inv(Matrix(X)))
    end
end

function LinearConstraint(A,b)
    AiAA = get_AiAA(A)
    P = -AiAA*A
    q = AiAA*b
    return LinearConstraint(sparse(A),sparse(b),sparse(P),sparse(q),spzeros(eltype(q),length(q)))
end
@inline function proj!(x,X::LinearConstraint)
    # X.rez .= x + X.P*x + X.q
    return x + X.P*x + X.q
end
@inline check(x,X::LinearConstraint) = all(abs.(X.A*x .- X.b) .<= 100*sqrt(eps(eltype(x))))

#################### DecomposedLinearConstraints
struct DecomposedLinearConstraint{T}
    A::SparseMatrixCSC{T,Int64}
    b::SparseVector{T, Int64}
    p::SparseVector{T, Int64}
    M::SparseMatrixCSC{Bool, Int64}
    q::SparseVector{T, Int64}
    rez::SparseVector{T, Int64}
end
function decompose(P;check=true)
    # check predicate and decompose:
    if check
        @assert all([sum(unique(P[i,:]).!=0)<=1 for i in 1:size(P,1)]) "This constraint does not seem to be decomposable. Maybe use LinearConstraint instead."
    end
    b = sparse([sum(unique(P[i,:])) for i in 1:size(P,1)])
    M = sparse(BitMatrix(P .!= 0))
    return b,M
end
function DecomposedLinearConstraint(A,b;check=true)
    AiAA = get_AiAA(A)
    P = sparse(-AiAA*A)
    p,M = decompose(P,check=check)
    q = sparse(+AiAA*b)
    return DecomposedLinearConstraint(A,sparse(b),p,M,q,spzeros(eltype(q),length(q)))
end
@inline function proj!(x,X::DecomposedLinearConstraint)
    #X.rez .= x + X.q + X.p.*(X.M*x)
    return x + X.q + X.p.*(X.M*x)
end
@inline check(x,X::DecomposedLinearConstraint) = all(abs.(X.A*x .- X.b) .<= 100*sqrt(eps(eltype(x))))

#################### PositivityConstraints
struct PositivityConstraint{T}
    rez::SparseVector{T, Int64}
    tol::T
end
function PositivityConstraint(x)
    return PositivityConstraint(spzeros(eltype(x),length(x)),sqrt(eps(eltype(x))))
end
@inline function proj!(x,X::PositivityConstraint)
    # X.rez .= x
    # SparseArrays.fkeep!(X.rez, (i,x) -> x>X.tol)
    SparseArrays.fkeep!(x, (i,x) -> x>X.tol)
    return x
end
@inline check(x,X::PositivityConstraint) = all(x.>=0)

function project!(x,ctrs;err=nothing)
    old_x = zero(x)
    err = isnothing(err) ? 10*sqrt(eps(eltype(x))) : err # warninng : this 10 should probably not be there;
    sn = [zero(x) for _ in ctrs]
    rez = zero(x)
    ctnue=!all([check(x,c) for c in ctrs])
    while ctnue
        old_x .= x
        for (i,ctr) in enumerate(ctrs)
            rez .= proj!(x-sn[i],ctr)
            sn[i] += rez - x
            x .= rez
        end
        ctnue = !all(abs.(old_x .-x) .< err)
        # display("$ctnue, $(sum(abs.(old_x .-x) .< err)) over $(length(old_x))")
    end
    return x
end

function linearize_thorin_transportation_problem(a::Vector{T},b::Vector{T}) where T
    n1 = length(a)
    n2 = length(b)
    p = (n1+1)*(n2+1)
    eq_A = zeros(T,(n1+n2,p-1))
    eq_b = [a...,b...]
    for i in 0:n2
        for j in 1:n1
            eq_A[j,j+i*(1+n1)] = 1
        end
    end
    for i in 1:n2
        eq_A[n1+i,(1:(n1+1)).+n1.+(i-1)*(n1+1)] .= 1
    end
    return eq_A,eq_b
end

end # module Dykstra


text:
## Various solver API's used in the past.  These functions are due to be standardized, and obsolete code / functions removed.

"""
    $SIGNATURES

Perform inference over the Bayes tree according to `opt::SolverParams`.

Notes
- Variety of options, including fixed-lag solving -- see `getSolverParams(fg)` for details.

Example
```julia
# without [or with] compute recycling
tree, smt, hist = solveTree!(fg [,tree])
```

Related

solveCliq!, wipeBuildNewTree!
"""
function solveTree!(dfgl::G,
                    oldtree::BayesTree=emptyBayesTree();
                    delaycliqs::Vector{Symbol}=Symbol[],
                    recordcliqs::Vector{Symbol}=Symbol[],
                    skipcliqids::Vector{Symbol}=Symbol[],
                    maxparallel::Int=50  ) where G <: DFG.AbstractDFG
  #
  @info "Solving over the Bayes (Junction) tree."
  smtasks=Vector{Task}()
  hist = Dict{Int, Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}}()
  opt = DFG.getSolverParams(dfgl)

  if opt.isfixedlag
      @info "Quasi fixed-lag is enabled (a feature currently in testing)!"
      fifoFreeze!(dfgl)
  end

  # current incremental solver builds a new tree and matches against old tree for recycling.
  tree = wipeBuildNewTree!(dfgl, drawpdf=opt.drawtree, show=opt.showtree, maxparallel=maxparallel, filepath=joinpath(getSolverParams(dfgl).logpath,"bt.pdf"))
  # setAllSolveFlags!(tree, false)

  @info "Do tree based init-inference on tree"
  if opt.async
    smtasks = asyncTreeInferUp!(dfgl, tree, oldtree=oldtree, N=opt.N, drawtree=opt.drawtree, recordcliqs=recordcliqs, limititers=opt.limititers, downsolve=opt.downsolve, incremental=opt.incremental, skipcliqids=skipcliqids, delaycliqs=delaycliqs )
  else
    smtasks, hist = initInferTreeUp!(dfgl, tree, oldtree=oldtree, N=opt.N, drawtree=opt.drawtree, recordcliqs=recordcliqs, limititers=opt.limititers, downsolve=opt.downsolve, incremental=opt.incremental, skipcliqids=skipcliqids, delaycliqs=delaycliqs )
  end
  @info "Finished tree based init-inference"

  # transfer new tree to outside parameter
  oldtree.bt = tree.bt
  oldtree.btid = tree.btid
  oldtree.cliques = tree.cliques
  oldtree.frontals = tree.frontals

  return oldtree, smtasks, hist
end

"""
    $SIGNATURES

Perform inference over one clique in the Bayes tree according to `opt::SolverParams`.

Example
```julia
tree = wipeBuildNewTree!(fg)
smt, hist = solveCliq!(fg, tree, :x1 [,cliqHistories=hist] )
```

Related

solveTree!, wipeBuildNewTree!
"""
function solveCliq!(dfgl::G,
                    tree::BayesTree,
                    cliqid::Symbol;
                    recordcliq::Bool=false,
                    # cliqHistories = Dict{Int,Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}}(),
                    maxparallel::Int=50  ) where G <: DFG.AbstractDFG
  #
  # hist = Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}()
  opt = DFG.getSolverParams(dfgl)

  if opt.isfixedlag
      @info "Quasi fixed-lag is enabled (a feature currently in testing)!"
      fifoFreeze!(dfgl)
  end

  # if !isTreeSolved(treel, skipinitialized=true)
  cliq = whichCliq(tree, cliqid)
  cliqtask = @async tryCliqStateMachineSolve!(dfgl, tree, cliq.index, drawtree=opt.drawtree, limititers=opt.limititers, downsolve=opt.downsolve,recordcliqs=(recordcliq ? [cliqid] : Symbol[]), incremental=opt.incremental) # N=N
  # end # if

  # post-hoc store possible state machine history in clique (without recursively saving earlier history inside state history)
  # assignTreeHistory!(tree, cliqHistories)

  # cliqHistories
  return cliqtask
end

"""
    $SIGNATURES

Perform up and down message passing (multi-process) algorithm for full sum-product solution of all continuous marginal beliefs.

Notes
- For legacy versions of tree traversal, see `inferOverTreeIterative!` instead.
"""
function inferOverTree!(dfg::G,
                        bt::BayesTree;
                        oldtree::BayesTree=emptyBayesTree(),
                        N::Int=100,
                        upsolve::Bool=true,
                        downsolve::Bool=true,
                        dbg::Bool=false,
                        drawpdf::Bool=false,
                        treeinit::Bool=false,
                        incremental::Bool=false,
                        limititers::Int=1000,
                        skipcliqids::Vector{Symbol}=Symbol[],
                        recordcliqs::Vector{Symbol}=Symbol[]  ) where G <: AbstractDFG
  #

  @info "Solving over the Bayes (Junction) tree."
  smtasks=Vector{Task}()
  ch = Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}()
  setAllSolveFlags!(bt, false)

  @info "Do tree based init-inference on tree"
  if dbg
    smtasks = asyncTreeInferUp!(dfg, bt, oldtree=oldtree, N=N, drawtree=drawpdf, recordcliqs=recordcliqs, limititers=limititers, downsolve=downsolve, incremental=incremental, skipcliqids=skipcliqids )
  else
    smtasks, ch = initInferTreeUp!(dfg, bt, oldtree=oldtree, N=N, drawtree=drawpdf, recordcliqs=recordcliqs, limititers=limititers, downsolve=downsolve, incremental=incremental, skipcliqids=skipcliqids )
  end
  @info "Finished tree based init-inference"

  return smtasks, ch
end

"""
    $SIGNATURES

Perform up and down message passing (multi-process) algorithm for full sum-product solution of all continuous marginal beliefs.

Notes
- Legacy support function, use `inferOverTree!` instead that is based on the state machine method.
- Previous versions of the code used iterative loops to traverse the Bayes (Junction) tree.
- Even older code is available as `inferOverTreeR!`
"""
function inferOverTreeIterative!(dfg::G,
                                 bt::BayesTree;
                                 N::Int=100,
                                 dbg::Bool=false,
                                 drawpdf::Bool=false  ) where G <: AbstractDFG
  #
  # @info "Batch rather than incremental solving over the Bayes (Junction) tree."
  # setAllSolveFlags!(bt, false)
  @info "Ensure all nodes are initialized"
  ensureAllInitialized!(dfg)
  @info "Do multi-process upward pass of inference on tree"
  upMsgPassingIterative!(ExploreTreeType(dfg, bt, bt.cliques[1], nothing, NBPMessage[]),N=N, dbg=dbg, drawpdf=drawpdf);
  @info "Do multi-process downward pass of inference on tree"
  downMsgPassingIterative!(ExploreTreeType(dfg, bt, bt.cliques[1], nothing, NBPMessage[]),N=N, dbg=dbg, drawpdf=drawpdf);
  return smtasks, ch
end

"""
    $SIGNATURES

Perform up and down message passing (single process, recursive) algorithm for full sum-product solution of all continuous marginal beliefs.
"""
function inferOverTreeR!(fgl::G,
                         bt::BayesTree;
                         N::Int=100,
                         dbg::Bool=false,
                         drawpdf::Bool=false,
                         treeinit::Bool=false  )::Tuple{Vector{Task},Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}} where G <: AbstractDFG
  #
  @info "Batch rather than incremental solving over the Bayes (Junction) tree."
  setAllSolveFlags!(bt, false)
  @info "Ensure all nodes are initialized"
  smtasks = Vector{Task}()
  ch = Vector{Tuple{DateTime, Int, Function, CliqStateMachineContainer}}()
  if treeinit
    smtasks, ch = initInferTreeUp!(fgl, bt, N=N, drawtree=drawpdf)
  else
    @info "Do conventional recursive up inference over tree"
    ensureAllInitialized!(fgl)
    upMsgPassingRecursive(ExploreTreeType(fgl, bt, bt.cliques[1], nothing, NBPMessage[]), N=N, dbg=dbg, drawpdf=drawpdf);
  end
  @info "Do recursive down inference over tree"
  downMsgPassingRecursive(ExploreTreeType(fgl, bt, bt.cliques[1], nothing, NBPMessage[]), N=N, dbg=dbg, drawpdf=drawpdf);
  return smtasks, ch
end

"""
    $(SIGNATURES)

Perform multimodal incremental smoothing and mapping (mm-iSAM) computations over given factor graph `fgl::FactorGraph` on the local computer.  A pdf of the Bayes (Junction) tree will be generated in the working folder with `drawpdf=true`
"""
function batchSolve!(dfg::G,
                     oldtree::BayesTree=emptyBayesTree();
                     upsolve::Bool=true,
                     downsolve::Bool=true,
                     drawpdf::Bool=false,
                     show::Bool=false,
                     N::Int=100,
                     recursive::Bool=false,
                     dbg::Bool=false,
                     treeinit::Bool=false,
                     incremental::Bool=false,
                     limititers::Int=1000,
                     skipcliqids::Vector{Symbol}=Symbol[],
                     recordcliqs::Vector{Symbol}=Symbol[],
                     returntasks::Bool=false,
                     maxparallel::Int=50  ) where G <: AbstractDFG
  #
  @warn "deprecated batchSolve! in favor of new solveTree! interface with the same and more functionality."
  if DFG.getSolverParams(dfg).isfixedlag
      @info "Quasi fixed-lag is enabled (a feature currently in testing)!"
      fifoFreeze!(dfg)
  end
  tree = wipeBuildNewTree!(dfg, drawpdf=false, maxparallel=maxparallel)
  drawpdf ? drawTree(tree, show=show) : nothing
  # show ? showTree() : nothing

  smtasks = Vector{Task}()
  if recursive
    # recursive is a single core method that is slower but occasionally helpful for better stack traces during debugging
    smtasks, ch = inferOverTreeR!(dfg, tree, N=N, drawpdf=drawpdf, dbg=dbg, treeinit=treeinit)
  else
    smtasks, ch = inferOverTree!(dfg, tree, oldtree=oldtree, N=N, drawpdf=drawpdf, dbg=dbg, treeinit=treeinit,
                                 limititers=limititers, recordcliqs=recordcliqs, upsolve=upsolve,
                                 downsolve=downsolve, incremental=incremental, skipcliqids=skipcliqids  )
  end

  # later development allows tasks for each cliq state machine to be returned also
  return tree, smtasks
end


text:
@testset "Differentials (on Sphere(2))" begin
# The Adjoint Differentials test using the same variables as the differentials
# test
x = SnPoint([1.,0.,0.])
y = SnPoint([0.,1.,0.])
M = Sphere(2)
ξ = log(M,x,y)
# Text differentials (1) Dx of Log_xy
@test AdjDxLog(M,x,x,ξ) == -ξ
@test AdjDyLog(M,x,x,ξ) == ξ
@test AdjDxExp(M,x,zeroTVector(M,x),ξ) == ξ
@test AdjDξExp(M,x,zeroTVector(M,x),ξ) == ξ
for t in [0,0.15,0.33,0.66,0.9]
    @test AdjDxGeo(M,x,x,t,ξ) == (1-t)*ξ
    @test norm(M,x,AdjDyGeo(M,x,x,t,ξ) - t*ξ)  ≈ 0 atol=10.0^(-16)
end
Mp = Power(M,2)
xP = PowPoint([x,y,x])
yP = PowPoint([x,x,y])
ξP = PowTVector([ξ, zeroTVector(M,x), -ξ])
@test norm(Mp,xP, AdjDforwardLogs(Mp,xP,ξP)
     - PowTVector([-ξ, zeroTVector(M,x),zeroTVector(M,x)])  )  ≈ 0 atol=3*10.0^(-16)
end


text:
import PipeLayout: serialmerge, pipesplit, isapproxin

@testset "split pipe utilities" begin
    diams = map(Diameter, [(1.0, 1.0), (1.2, 1.2), (1.4, 1.5)])

    @test serialmerge(diams, [1.0, 0.0, 0.0]) ≈ 1.0
    @test serialmerge(diams, [0.0, 1.0, 0.0]) ≈ 1.2
    @test serialmerge(diams, [0.0, 0.0, 1.0]) ≈ 1.4
    d = serialmerge(diams, [0.5, 0.5, 0.0])
    @test isapproxin(d, 1.0, 1.2)

    @test pipesplit(diams, 1.0) ≈ [1.0, 0.0, 0.0]
    @test pipesplit(diams, 1.2) ≈ [0.0, 1.0, 0.0]
    @test pipesplit(diams, 1.4) ≈ [0.0, 0.0, 1.0]
    @test pipesplit(diams, d)   ≈ [0.5, 0.5, 0.0]
end


text:
using Flux
using Flux: onehot, chunk, batchseq, throttle, logitcrossentropy
using StatsBase: wsample
using Base.Iterators: partition
using Dates: today
import Downloads
import BSON

# hyperparameters
Base.@kwdef struct Args
    η::Float32 = 1e-2           # learning rate
    seqlen::Integer = 50        # length of batch sequences
    nbatch::Integer = 50        # number of batches
    throttle::Integer = 30      # throttle timeout
    savepath = "runs/"          # runs saved to this directory
end

function get_data(args)
    # Download the data in this directory as "input.txt"
    url = "https://cs.stanford.edu/people/karpathy/char-rnn/shakespeare_input.txt"
    isfile("input.txt") || Downloads.download(url, "input.txt")

    # Collect all the characters in the text
    characters = collect(String(read("input.txt")))

    # Getting an alphabet of all unique characters + the underscore
    alphabet = [unique(characters)..., '_']

    text = [onehot(ch, alphabet) for ch in characters]
    stop = onehot('_', alphabet)

    N = length(alphabet)

    # partitioning the data as a sequence of batches,
    # which are collected into an array of batches.
    Xs = collect(partition(batchseq(chunk(text, args.nbatch), stop), args.seqlen))
    Ys = collect(partition(batchseq(chunk(text[2:end], args.nbatch), stop), args.seqlen))

    Xs, Ys, N, alphabet
end

function model(N::Integer)
    Chain(
        LSTM(N, 128),
        LSTM(128, 128),
        Dense(128, N)
    )
end

function train(; kws...)
    # initialize hyperparameters
    args = Args(; kws...)

    # get data
    Xs, Ys, N, alphabet = get_data(args)

    # construct model
    m = model(N)

    # loss and optimization functions
    loss(xs, ys) = logitcrossentropy.([m(x) for x ∈ xs], ys) |> sum
    opt = ADAM(args.η)

    evalcb = () -> @show loss(Xs[5], Ys[5])

    # Training
    @info "Training begun."
    Flux.train!(loss, params(m), zip(Xs, Ys), opt, cb=throttle(evalcb, args.throttle))

    # Saving model
    !ispath(args.savepath) && mkpath(args.savepath)
    modelpath = joinpath(args.savepath, "model$(today()).bson")
    let model = cpu(m)
        BSON.@save modelpath model
    end

    m, alphabet
end

function sample(m, alphabet, len; seed="")
    m = cpu(m)
    Flux.reset!(m)
    buf = IOBuffer()
    seed == "" && (seed = string(rand(alphabet)))
    write(buf, seed)
    c = wsample(alphabet, softmax([m(onehot(c, alphabet)) for c ∈ collect(seed)][end]))
    for i ∈ 1:len
        write(buf, c)
        c = wsample(alphabet, softmax(m(onehot(c, alphabet))))
    end
    String(take!(buf))
end

if abspath(PROGRAM_FILE) == @__FILE__
    cd(@__DIR__)
    m, alphabet = @time train()
    sample(m, alphabet, 1000) |> println
end

max_stars_repo_name: benharmonics/ModelZoo

max_issues_repo_name: benharmonics/ModelZoo

max_forks_repo_name: benharmonics/ModelZoo