text:
using MultivariateOrthogonalPolynomials, StaticArrays, BlockArrays, BlockBandedMatrices, ArrayLayouts,
        QuasiArrays, Test, ClassicalOrthogonalPolynomials, BandedMatrices, FastTransforms, LinearAlgebra
import MultivariateOrthogonalPolynomials: tri_forwardrecurrence, grid, TriangleRecurrenceA, TriangleRecurrenceB, TriangleRecurrenceC, xy_muladd!

@testset "Triangle" begin
    @testset "basics" begin
        P = JacobiTriangle()
        @test copy(P) ≡ P

        xy = axes(P,1)
        x,y = first.(xy),last.(xy)
        @test xy[SVector(0.1,0.2)] == SVector(0.1,0.2)
        @test x[SVector(0.1,0.2)] == 0.1
        @test y[SVector(0.1,0.2)] == 0.2
    end
    p = (n,k,a,b,c,x,y) -> jacobip(n-k,2k+b+c+1,a,2x-1) * (1-x)^k * jacobip(k,c,b,2y/(1-x)-1)

    @testset "evaluation" begin
        @testset "versus explicit" begin
            x,y = xy = SVector(0.1,0.2)
            for (a,b,c) in ((0,0,0), (1,0,0), (0,1,0), (0,0,1), (0.1,0.2,0.3))
                P = JacobiTriangle(a,b,c)

                for n = 0:5, k=0:n
                    @test P[xy,Block(n+1)[k+1]] ≈ p(n,k,a,b,c,x,y) atol=1E-13
                end
            end
        end

        @testset "forwardrecurrnce" begin
            P = JacobiTriangle()
            xy = SVector(0.1,0.2)
            P_N = P[xy, Block.(Base.OneTo(10))]
            @test P_N == P[xy,Block.(1:10)]
            for N = 1:10
                @test P[xy, Block(N)] ≈ P_N[Block(N)] ≈ p.(N-1, 0:N-1, 0,0,0, xy...)
                for j=1:N
                    P[xy,Block(N)[j]] ≈ p(N-1,j-1,0,0,0,xy...)
                end
            end
            @test P[xy,1] == 1
            @test P[xy,2] ≈ p(1,0,0,0,0,xy...)
            @test P[xy,3] ≈ p(1,1,0,0,0,xy...)
            @test P[xy,4] ≈ p(2,0,0,0,0,xy...)

            @test P[[SVector(0.1,0.2),SVector(0.2,0.3)],1] ≈ [1,1]
            @test P[[SVector(0.1,0.2),SVector(0.2,0.3)],Block(2)] ≈ [P[SVector(0.1,0.2),2] P[SVector(0.1,0.2),3];
                                                                    P[SVector(0.2,0.3),2] P[SVector(0.2,0.3),3]]
        end
        @testset "function" begin
            P = JacobiTriangle()
            xy = SVector(0.1,0.2)
            c = PseudoBlockVector([1; Zeros(∞)], (axes(P,2),))
            f = P*c
            @test @inferred(f[xy]) == 1.0
            c = PseudoBlockVector([1:3; Zeros(∞)], (axes(P,2),))
            f = P*c
            @test f[xy] ≈ P[xy,1:3]'*(1:3)
            c = PseudoBlockVector([1:6; Zeros(∞)], (axes(P,2),))
            f = P*c
            @test f[xy] ≈ P[xy,1:6]'*(1:6)

            c = PseudoBlockVector([randn(5050); Zeros(∞)], (axes(P,2),))
            f = P*c
            @test f[xy] ≈ P[xy,1:5050]'*c[1:5050]

            c = PseudoBlockVector([1:10; zeros(∞)], (axes(P,2),))
            f = P*c
            xy = SVector(0.1,0.2)
            @test f[xy] ≈ dot(P[xy,1:10],1:10)
            @test f[[xy, xy.+0.1]] ≈ [f[xy], f[xy.+0.1]]
            @test f[permutedims([xy, xy.+0.1])] ≈ [f[xy] f[xy.+0.1]]

            @testset "block structure missing" begin
                f = P * [1:5; zeros(∞)]
                @test f.args[2][Block(2)] == 2:3
                @test f[xy] ≈ P[xy,1:5]'*(1:5)

                f = P * [1:5; Zeros(∞)]
                @test f.args[2][Block(2)] == 2:3
                @test f[xy] ≈ P[xy,1:5]'*(1:5)
            end
        end
    end

    @testset "transform" begin
        @testset "grid" begin
            P = JacobiTriangle()
            N = M = 20
            P_N = P[:,Block.(Base.OneTo(N))]
            x = [sinpi((2N-2n-1)/(4N))^2 for n in 0:N-1]
            w = [sinpi((2M-2m-1)/(4M))^2 for m in 0:M-1]
            g = [SVector(x[n+1], x[N-n]*w[m+1]) for n in 0:N-1, m in 0:M-1]
            @test grid(P_N) ≈ g
        end

        @testset "relation with transform" begin
            P = JacobiTriangle()
            c = PseudoBlockVector([1:10; zeros(∞)], (axes(P,2),))
            f = P*c
            N = 5
            P_N = P[:,Block.(Base.OneTo(N))]
            g = grid(P_N)
            F = f[g]
            Pl = plan_tri2cheb(F, 0, 0, 0)
            PA = plan_tri_analysis(F)
            @time U = Pl\(PA*F)
            @test MultivariateOrthogonalPolynomials.tridenormalize!(U,0,0,0) ≈ [1 3 6 10 0; 2 5 9 0 0; 4 8 0 0 0; 7 0 0 0 0; 0 0 0 0 0]
        end

        @testset "expansions" begin
            P = JacobiTriangle()
            xy = axes(P,1)
            x,y = first.(xy),last.(xy)
            N = 20
            P_N = P[:,Block.(Base.OneTo(N))]
            u = P_N * (P_N \ (exp.(x) .* cos.(y)))
            @test u[SVector(0.1,0.2)] ≈ exp(0.1)*cos(0.2)

            P_n = P[:,1:200]
            u = P_n * (P_n \ (exp.(x) .* cos.(y)))
            @test u[SVector(0.1,0.2)] ≈ exp(0.1)*cos(0.2)

            @time u = P * (P \ (exp.(x) .* cos.(y)))
            @test u[SVector(0.1,0.2)] ≈ exp(0.1)*cos(0.2)
        end
    end

    @testset "operators" begin
        P = JacobiTriangle()
        xy = axes(P,1)

        ∂ˣ = PartialDerivative{1}(xy)
        ∂ʸ = PartialDerivative{2}(xy)

        @test eltype(∂ˣ) == eltype(∂ʸ) == Float64

        Dˣ = JacobiTriangle(1,0,1) \ (∂ˣ * P)
        Dʸ = JacobiTriangle(0,1,1) \ (∂ʸ * P)

        M = P'P
        @test blockbandwidths(M) == subblockbandwidths(M) == (0,0)
        @test M[1,1] ≈ 1/2
        @test M[2,2] ≈ 1/4
        @test M[3,3] ≈ 1/12
        @test M[4,4] ≈ 1/6
        @test M[5,5] ≈ 1/18

        Rx = JacobiTriangle(1,0,0) \ P
        Ry = JacobiTriangle(0,1,0) \ P
        Rz = JacobiTriangle(0,0,1) \ P

        x,y = 0.1,0.2

        for n=1:5, k=0:n-1
            @test p(n,k,0,0,0,x,y) ≈ p(n-1,k,1,0,0,x,y) *  Rx[Block(n)[k+1], Block(n+1)[k+1]] + p(n,k,1,0,0,x,y) *  Rx[Block(n+1)[k+1], Block(n+1)[k+1]]
        end
        for n=0:5
            k = n
            @test p(n,k,0,0,0,x,y) ≈ p(n,k,1,0,0,x,y) *  Rx[Block(n+1)[k+1], Block(n+1)[k+1]]
        end

        @test P[SVector(x,y),Block.(1:5)]' ≈ JacobiTriangle(1,0,0)[SVector(x,y),Block.(1:5)]' * Rx[Block.(1:5),Block.(1:5)]
        @test P[SVector(x,y),Block.(1:5)]' ≈ JacobiTriangle(0,1,0)[SVector(x,y),Block.(1:5)]' * Ry[Block.(1:5),Block.(1:5)]
        @test P[SVector(x,y),Block.(1:5)]' ≈ JacobiTriangle(0,0,1)[SVector(x,y),Block.(1:5)]' * Rz[Block.(1:5),Block.(1:5)]

        c = [randn(100); zeros(∞)]
        @test (P*c)[SVector(x,y)] ≈ (JacobiTriangle(1,0,0)*(Rx*c))[SVector(x,y)]
        @test (P*c)[SVector(x,y)] ≈ (JacobiTriangle(0,1,0)*(Ry*c))[SVector(x,y)]
        @test (P*c)[SVector(x,y)] ≈ (JacobiTriangle(0,0,1)*(Rz*c))[SVector(x,y)]

        Lx = P \ WeightedTriangle(1,0,0)
        Ly = P \ WeightedTriangle(0,1,0)
        Lz = P \ WeightedTriangle(0,0,1)

        for n=0:5, k=0:n
            @test x*p(n,k,1,0,0,x,y) ≈ p(n,k,0,0,0,x,y) *  Lx[Block(n+1)[k+1], Block(n+1)[k+1]] + p(n+1,k,0,0,0,x,y) *  Lx[Block(n+2)[k+1], Block(n+1)[k+1]]
        end

        @test WeightedTriangle(1,0,0)[SVector(x,y),Block.(1:5)]' ≈ P[SVector(x,y),Block.(1:6)]' * Lx[Block.(1:6),Block.(1:5)]
        @test WeightedTriangle(0,1,0)[SVector(x,y),Block.(1:5)]' ≈ P[SVector(x,y),Block.(1:6)]' * Ly[Block.(1:6),Block.(1:5)]
        @test WeightedTriangle(0,0,1)[SVector(x,y),Block.(1:5)]' ≈ P[SVector(x,y),Block.(1:6)]' * Lz[Block.(1:6),Block.(1:5)]

        @test (WeightedTriangle(1,0,0)*c)[SVector(x,y)] ≈ (P*(Lx*c))[SVector(x,y)]
        @test (WeightedTriangle(0,1,0)*c)[SVector(x,y)] ≈ (P*(Ly*c))[SVector(x,y)]
        @test (WeightedTriangle(0,0,1)*c)[SVector(x,y)] ≈ (P*(Lz*c))[SVector(x,y)]

        ∂ˣ² = (∂ˣ)^2
        ∂ʸ² = (∂ʸ)^2
        @test ∂ˣ² isa ApplyQuasiMatrix{<:Any,typeof(^)}
        @test ∂ʸ² isa ApplyQuasiMatrix{<:Any,typeof(^)}

        Dˣ² = JacobiTriangle(2,0,2) \ (∂ˣ² * P)
        Dˣ² = JacobiTriangle(2,0,2) \ (∂ˣ * (∂ˣ * P))
        Dʸ² = JacobiTriangle(0,2,2) \ (∂ʸ * (∂ʸ * P))

        @testset "jacobi" begin
            P = JacobiTriangle()
            xy = axes(P,1)
            x,y = first.(xy),last.(xy)
            X = P \ (x .* P)
            Y = P \ (y .* P)

            N = 100
            @test X.args[1][Block.(Base.OneTo(N)), Block.(Base.OneTo(N))] isa BandedBlockBandedMatrix

            @test X[Block(2,2)] isa BandedMatrix

            @test X[Block.(1:5),Block.(1:5)] ≈ (Lx * Rx)[Block.(1:5),Block.(1:5)]
            @test Y[Block.(1:5),Block.(1:5)] ≈ (Ly * Ry)[Block.(1:5),Block.(1:5)]

            @test blockbandwidths(X) == blockbandwidths(Y) == (1,1)
            @test subblockbandwidths(X) == (0,0)
            @test subblockbandwidths(Y) == (1,1)

            @testset "Test with exact" begin
                x,y = 0.1,0.2
                N = 0
                @test [p.(N,0:N,0,0,0,x,y); p.(N+1,0:N+1,0,0,0,x,y)]'* X[Block.(N+1:N+2),1] ≈ x
                @test [p.(N,0:N,0,0,0,x,y); p.(N+1,0:N+1,0,0,0,x,y)]'* Y[Block.(N+1:N+2),1] ≈ y
                for N = 1:5
                    @test [p.(N-1,0:N-1,0,0,0,x,y); p.(N,0:N,0,0,0,x,y); p.(N+1,0:N+1,0,0,0,x,y)]'* X[Block.(N:N+2),Block(N+1)] ≈ x * p.(N,0:N,0,0,0,x,y)'
                    @test [p.(N-1,0:N-1,0,0,0,x,y); p.(N,0:N,0,0,0,x,y); p.(N+1,0:N+1,0,0,0,x,y)]'* Y[Block.(N:N+2),Block(N+1)] ≈ y * p.(N,0:N,0,0,0,x,y)'
                end
            end

            @testset "Recurrence" begin
                @testset "B⁺" begin
                    B⁺ = function(N)
                        Bˣ = X[Block(N+1,N)]'
                        Bʸ = Y[Block(N+1,N)]'
                        B¹ = Bˣ[:,1:N]
                        b₂ = Bʸ[N,N+1]
                        b₁ = if N == 1
                            -b₂ \ [Bʸ[N,N]/Bˣ[N,N]]
                        else
                            -b₂ \ [zeros(N-2); Bʸ[N,N-1]/Bˣ[N-1,N-1]; Bʸ[N,N]/Bˣ[N,N]]
                        end
                        [inv(B¹) zeros(N,N-1) zeros(N);
                            b₁' zeros(1,N-1) inv(b₂)]
                    end

                    # B⁺! = function (w, N, v)
                    #     w[1:N] .= X[Block(N+1,N)][band(0)] .\ view(v,1:N)
                    # end

                    for N = 1:5
                        Bˣ = X[Block(N+1,N)]'; Bʸ = Y[Block(N+1,N)]'; B = [Bˣ; Bʸ]
                        Aˣ = X[Block(N,N)]'; Aʸ = Y[Block(N,N)]'; A = [Aˣ; Aʸ]
                        if N > 1
                            Cˣ = X[Block(N-1,N)]'; Cʸ = Y[Block(N-1,N)]'; C = [Cˣ; Cʸ]
                        end
                        @test B⁺(N) * B ≈ I
                        @test TriangleRecurrenceA(N, X, Y) ≈ B⁺(N)
                        @test TriangleRecurrenceB(N, X, Y) ≈ B⁺(N)*[-Aˣ; -Aʸ]
                        if N > 1
                            @test TriangleRecurrenceC(N, X, Y) ≈ B⁺(N)*[Cˣ; Cʸ]
                        end

                        x,y = 0.1,0.2
                        v = randn(N)
                        w = randn(N+1)
                        @test xy_muladd!((x,y), TriangleRecurrenceA(N,X,Y),  v, 2.0, copy(w)) ≈
                            B⁺(N) * [x*Eye(N); y*Eye(N)]*v + 2w
                        @test xy_muladd!((x,y), TriangleRecurrenceA(N,X,Y)', w, 2.0, copy(v)) ≈
                            (B⁺(N) * [x*Eye(N); y*Eye(N)])'*w + 2v

                        @test mul!(w, TriangleRecurrenceB(N, X, Y), v) ≈ B⁺(N)*[-Aˣ; -Aʸ]*v
                        @test mul!(v, TriangleRecurrenceB(N, X, Y)', w) ≈ (B⁺(N)*[-Aˣ; -Aʸ])'*w
                        @test muladd!(3.0, TriangleRecurrenceB(N, X, Y), v, 2.0, copy(w)) ≈ 3B⁺(N)*[-Aˣ; -Aʸ]*v + 2w
                        @test muladd!(3.0, TriangleRecurrenceB(N, X, Y)', w, 2.0, copy(v)) ≈ 3(B⁺(N)*[-Aˣ; -Aʸ])'*w + 2v

                        if N > 1
                            u = randn(N-1)
                            @test muladd!(3.0, TriangleRecurrenceC(N, X, Y), u, 2.0, copy(w)) ≈
                                3B⁺(N)*[Cˣ; Cʸ]*u + 2w
                            @test muladd!(3.0, TriangleRecurrenceC(N, X, Y)', w, 2.0, copy(u)) ≈
                                3(B⁺(N)*[Cˣ; Cʸ])'*w + 2u

                            # need in-place to minimise buffers in Clenshaw
                            w̃ = copy(w)
                            @test lmul!(TriangleRecurrenceC(N, X, Y)', w̃) === w̃
                            @test w̃ ≈ (B⁺(N)*[Cˣ; Cʸ])'*w
                        end
                    end

                    @testset "comparison with exact" begin
                        x,y = 0.1,0.2

                        N = 1
                        Aˣ = X[Block(N,N)]'
                        Aʸ = Y[Block(N,N)]'
                        P_0 = [1]
                        P_1 = B⁺(1)*[x*I-Aˣ; y*I-Aʸ]*P_0
                        @test P_1 ≈ p.(N,0:N,0,0,0,x,y)

                        N = 2
                        Aˣ = X[Block(N,N)]'; Aʸ = Y[Block(N,N)]'; A = [Aˣ-x*I; Aʸ-y*I]
                        Bˣ = X[Block(N+1,N)]'; Bʸ = Y[Block(N+1,N)]'; B = [Bˣ; Bʸ]
                        Cˣ = X[Block(N-1,N)]'; Cʸ = Y[Block(N-1,N)]'; C = [Cˣ; Cʸ]
                        @test norm(C * p.(N-2,0:N-2,0,0,0,x,y) + A * p.(N-1,0:N-1,0,0,0,x,y) + B * p.(N,0:N,0,0,0,x,y)) ≤ 10eps()

                        P_2 = B⁺(N)*([x*I-Aˣ; y*I-Aʸ]*P_1 - [Cˣ; Cʸ]*P_0 )
                        @test p.(N,0:N,0,0,0,x,y) ≈ P_2

                        N = 3
                        Aˣ = X[Block(N,N)]'; Aʸ = Y[Block(N,N)]'; A = [Aˣ-x*I; Aʸ-y*I]
                        Cˣ = X[Block(N-1,N)]'; Cʸ = Y[Block(N-1,N)]'; C = [Cˣ; Cʸ]
                        P_3 = B⁺(N)*([x*I-Aˣ; y*I-Aʸ]*P_2 - [Cˣ; Cʸ]*P_1 )
                        @test p.(N,0:N,0,0,0,x,y) ≈ P_3

                        @testset "simplify" begin
                            Bˣ = X[Block(N+1,N)]'; Bʸ = Y[Block(N+1,N)]'; B = [Bˣ; Bʸ]

                            b₂ = Bʸ[N,N+1]
                            Bc = -Bʸ[N,N-1]/(b₂*Bˣ[N-1,N-1])*Cˣ[N-1,N-1] + Cʸ[N,N-1]/b₂
                            @test B⁺(N) * [Cˣ; Cʸ] ≈ [Diagonal(Bˣ[band(0)][1:end-1] .\ Cˣ[band(0)]);
                                                        zeros(1, N-1);
                                                      [zeros(1,N-2) Bc]]

                            Ba_1 = -Bʸ[N,N-1]/(b₂*Bˣ[N-1,N-1]) * (x - Aˣ[N-1,N-1]) + (-Aʸ[N,N-1])/b₂
                            Ba_2 = -Bʸ[N,N]/(b₂*Bˣ[N,N]) * (x - Aˣ[N,N]) + (y-Aʸ[N,N])/b₂

                            @test B⁺(N) * [x*I-Aˣ; y*I-Aʸ] ≈ [Diagonal(Bˣ[band(0)] .\ (x .- Aˣ[band(0)]));
                                                             [zeros(1,N-2) Ba_1 Ba_2]]

                            w = Vector{Float64}(undef,N+1)
                            w[1:N-1] .= (Bˣ[band(0)][1:end-1] .\ Cˣ[band(0)] .* P_1)
                            w[N] = 0
                            w[N+1] = Bc * P_1[end]
                            @test w ≈ B⁺(N) * [Cˣ; Cʸ] * P_1
                            w[1:N] .= Bˣ[band(0)] .\ (x .- Aˣ[band(0)]) .* P_2
                            w[N+1] = Ba_1 * P_2[end-1] + Ba_2 * P_2[end]
                            @test w ≈ B⁺(N) * [x*I-Aˣ; y*I-Aʸ] * P_2
                            w[1:N-1] .-= (Bˣ[band(0)][1:end-1] .\ Cˣ[band(0)] .* P_1)
                            w[N+1] -= Bc * P_1[end]
                            @test w ≈ P_3
                        end
                    end
                end

                @testset "forward Recurrence" begin
                    @test tri_forwardrecurrence(1, X, Y, 0.1, 0.2) ≈ [1.0]
                    @test tri_forwardrecurrence(2, X, Y, 0.1, 0.2) ≈ [1.0; p.(1,0:1,0,0,0,0.1,0.2)]
                    @test tri_forwardrecurrence(4, X, Y, 0.1, 0.2)[Block(4)] ≈ p.(3,0:3,0,0,0,0.1,0.2)
                end
            end

            @testset "truncations" begin
                N = 100;
                KR,JR = Block.(1:N),Block.(1:N)

                @test Dˣ[KR,JR] isa BandedBlockBandedMatrix
                @test Lx[KR,JR] isa BandedBlockBandedMatrix
                @test Rx[KR,JR] isa BandedBlockBandedMatrix
                @test Ly[KR,JR] isa BandedBlockBandedMatrix
                @test Ry[KR,JR] isa BandedBlockBandedMatrix
                @test X[KR,JR] isa BandedBlockBandedMatrix
                @test Y[KR,JR] isa BandedBlockBandedMatrix
            end

            @testset "other parameters" begin
                P = JacobiTriangle(1,0,0)
                xy = axes(P,1)
                x,y = first.(xy),last.(xy)
                X = P \ (x .* P)
                Y = P \ (y .* P)

                P_ex = PseudoBlockVector{Float64}(undef, 1:5)
                for n = 0:4, k=0:n
                    P_ex[Block(n+1)[k+1]] = p(n,k,1,0,0,0.1,0.2)
                end
                @test P_ex'*X[Block.(1:5),Block.(1:4)] ≈ 0.1 * P_ex[Block.(1:4)]'
                @test P_ex'*Y[Block.(1:5),Block.(1:4)] ≈ 0.2 * P_ex[Block.(1:4)]'
                @test P[SVector(0.1,0.2),Block.(1:5)]'*X[Block.(1:5),Block.(1:4)] ≈ 0.1 * P[SVector(0.1,0.2),Block.(1:4)]'
            end
        end

        @testset "higher order conversion" begin
            P = JacobiTriangle()
            Q = JacobiTriangle(1,1,1)
            R = Q \ P
            L = P \ Weighted(Q)

            xy = SVector(0.1,0.2)
            @test P[xy,1:10]' ≈ Q[xy,1:10]' * R[1:10,1:10]
            @test Weighted(Q)[xy,1:10]' ≈ P[xy,1:50]'*L[1:50,1:10]

            Q = JacobiTriangle(0,0,2)
            R = Q \ P
            L = P \ Weighted(Q)
            xy = SVector(0.1,0.2)
            @test P[xy,1:10]' ≈ Q[xy,1:10]' * R[1:10,1:10]
            @test Weighted(Q)[xy,1:10]' ≈ P[xy,1:50]'*L[1:50,1:10]
        end

        @testset "WeightedBasis" begin
            P = JacobiTriangle()
            Q = JacobiTriangle(1,1,1)
            w_0 = TriangleWeight(0,0,0)
            @test all(isone,w_0)
            @test !all(isone,TriangleWeight(0,0,1))
            @test w_0 == w_0
            @test w_0 .* P == w_0 .* P
            @test w_0 .* P == P
            @test P == w_0 .* P
            @test w_0 .* P == Weighted(P)
            @test Weighted(P)== w_0 .* P
            @test P == Weighted(P)
            @test Weighted(P)== P

            @test Weighted(P) \ Weighted(P) == P \ Weighted(P) == Weighted(P) \ P ==
                        P \ P == (w_0 .* P) \ P == P \ (w_0 .* P) == (w_0 .* P) \ (w_0 .* P) ==
                        (w_0 .* P) \ Weighted(P) == Weighted(P) \ (w_0 .* P)

            @test ((w_0 .* Q) \ P)[1:10,1:10] == ((w_0 .* Q) \ (w_0 .* P))[1:10,1:10] == (Q \ (w_0 .* P))[1:10,1:10] == (Q \ P)[1:10,1:10]
        end

        @testset "general (broken)" begin
            P = JacobiTriangle()
            Q = JacobiTriangle(0.1,0.2,0.3)
            w = TriangleWeight(0.2,0.3,0.4)
            @test_throws ErrorException Q\P
            @test_throws ErrorException Weighted(Q)\Weighted(P)
            @test_throws ErrorException (w.*Q)\(w.*P)
            @test_throws ErrorException Weighted(Q)\P
            @test_throws ErrorException Q\Weighted(P)
            @test_throws ErrorException (w .* Q)\P
            @test_throws ErrorException Q\( w.* P)
            @test_throws ErrorException Weighted(Q)\(w .* P)
            @test_throws ErrorException (w .* Q)\Weighted(P)
        end
    end

    @testset "AngularMomentum" begin
        P = JacobiTriangle()
        P¹ = JacobiTriangle(1,1,1)
        xy = axes(P,1)
        x,y = first.(xy),last.(xy)
        ∂ˣ = PartialDerivative{1}(xy)
        ∂ʸ = PartialDerivative{2}(xy)
        L1 = x .* ∂ʸ
        L2 = y .* ∂ˣ
        L = x .* ∂ʸ - y .* ∂ˣ
        A = P¹ \ (L1 * P)
        B = P¹ \ (L2 * P)
        C = P¹ \ (L * P)
        @test C[1:10,1:10] ≈ A[1:10,1:10] - B[1:10,1:10]
    end
end


text:
function get_variable_info(model::JuMP.Model)
    varinds = OrderedDict()
    integer = OrderedDict()
    lb = OrderedDict()
    ub = OrderedDict()
    inits = OrderedDict()
    for (k, vs) in model.obj_dict
        if vs isa VariableRef
            v = vs
            ind = v.index.value
            varinds[k] = ind
            integer[k] = is_binary(v) || is_integer(v)
            lb[k] = if has_lower_bound(v)
                lower_bound(v)
            elseif is_binary(v)
                0.0
            else
                -Inf
            end
            ub[k] = if has_upper_bound(v)
                upper_bound(v)
            elseif is_binary(v)
                1.0
            else
                Inf
            end
            inits[k] = if start_value(v) !== nothing
                start_value(v)
            elseif ub[k] == Inf && lb[k] == -Inf
                0.0
            elseif ub[k] == Inf
                lb[k] + 1.0
            elseif lb[k] == -Inf
                lb[k] - 1.0
            else
                (ub[k] + lb[k]) / 2
            end
        elseif vs isa AbstractArray{<:VariableRef}
            inds = map(v -> v.index.value, vs)
            if inds isa JuMP.Containers.DenseAxisArray
                varinds[k] = OrderedDict(inds.data .=> keys(inds))
            else
                varinds[k] = OrderedDict(vec(inds) .=> 1:length(inds))
            end
            integer[k] = map(vs) do v
                is_binary(v) || is_integer(v)
            end
            lb[k] = map(vs) do v
                if has_lower_bound(v)
                    lower_bound(v)
                elseif is_binary(v)
                    0.0
                else
                    -Inf
                end
            end
            ub[k] = map(vs) do v
                if has_upper_bound(v)
                    upper_bound(v)
                elseif is_binary(v)
                    1.0
                else
                    Inf
                end
            end
            inits[k] = map(vs) do v
                ind = varinds[k][v.index.value]
                if start_value(v) !== nothing
                    start_value(v)
                elseif ub[k][ind] == Inf && lb[k][ind] == -Inf
                    0.0
                elseif ub[k][ind] == Inf
                    lb[k][ind] + 1.0
                elseif lb[k][ind] == -Inf
                    lb[k][ind] - 1.0
                else
                    (ub[k][ind] + lb[k][ind]) / 2
                end
            end
        else
            continue
        end
    end
    sort!(integer)
    sort!(lb)
    sort!(ub)
    sort!(inits)
    _lb = OrderedDict(k => v for (k, v) in lb)
    _ub = OrderedDict(k => v for (k, v) in ub)
    _inits = OrderedDict(k => v for (k, v) in inits)
    _integer = OrderedDict(k => v for (k, v) in integer)
    return _lb, _ub, _inits, _integer
end

function get_constraint_info(model::JuMP.Model, nvars)
    Iineq = Int[]
    Jineq = Int[]
    Vineq = Float64[]
    bineq = Float64[]

    Ieq = Int[]
    Jeq = Int[]
    Veq = Float64[]
    beq = Float64[]

    ineqcounter = 0
    eqcounter = 0
    inv_geq_constraints = OrderedDict()
    inv_leq_constraints = OrderedDict()
    inv_eq_constraints = OrderedDict()
    for (_, c) in model.obj_dict
        if c isa ConstraintRef
            cs = [c]
        elseif c isa AbstractArray{<:ConstraintRef}
            cs = c
        else
            cs = nothing
            continue
        end
        foreach(cs) do c
            set = constraint_object(c).set
            if set isa MOI.GreaterThan
                inv_geq_constraints[JuMP.index(c).value] = c
            elseif set isa MOI.LessThan
                inv_leq_constraints[JuMP.index(c).value] = c
            elseif set isa MOI.EqualTo
                inv_eq_constraints[JuMP.index(c).value] = c
            else
                throw("Set type not supported.")
            end
        end
    end
    sort!(inv_geq_constraints)
    sort!(inv_leq_constraints)
    sort!(inv_eq_constraints)
    geq_constraints = values(inv_geq_constraints)
    leq_constraints = values(inv_leq_constraints)
    eq_constraints = values(inv_eq_constraints)

    for c in geq_constraints
        func = constraint_object(c).func
        @assert func isa AffExpr
        set = constraint_object(c).set
        ineqcounter += 1
        push!(bineq, -set.lower)
        for (var, val) in func.terms
            push!(Iineq, ineqcounter)
            push!(Jineq, var.index.value)
            push!(Vineq, -val)
        end
    end
    for c in leq_constraints
        func = constraint_object(c).func
        @assert func isa AffExpr
        set = constraint_object(c).set
        ineqcounter += 1
        push!(bineq, set.upper)
        for (var, val) in func.terms
            push!(Iineq, ineqcounter)
            push!(Jineq, var.index.value)
            push!(Vineq, val)
        end
    end
    for c in eq_constraints
        func = constraint_object(c).func
        @assert func isa AffExpr
        set = constraint_object(c).set
        eqcounter += 1
        push!(beq, set.value)
        for (var, val) in func.terms
            push!(Ieq, eqcounter)
            push!(Jeq, var.index.value)
            push!(Veq, val)
        end
    end
    Aineq = sparse(Iineq, Jineq, Vineq, ineqcounter, nvars)
    Aeq = sparse(Ieq, Jeq, Veq, eqcounter, nvars)

    return Aineq, bineq, Aeq, beq
end

function get_objective_info(model::JuMP.Model, nvars)
    obj = objective_function(model)
    @assert obj isa AffExpr
    sense = objective_sense(model)
    I = Int[]
    V = Float64[]
    for (var, val) in obj.terms
        push!(I, var.index.value)
        if sense == MathOptInterface.MIN_SENSE
            push!(V, val)
        elseif sense == MathOptInterface.MAX_SENSE
            push!(V, -val)
        end
    end
    c = sparsevec(I, V, nvars)
    return c
end

function DictModel(model::JuMP.Model)
    lb, ub, inits, integer = get_variable_info(model)
    nvars = length(flatten(lb)[1])
    Aineq, bineq, Aeq, beq = get_constraint_info(model, nvars)
    c = get_objective_info(model, nvars)
    ks = collect(keys(inits))
    dict_model = DictModel()
    for k in ks
        addvar!(dict_model, k, lb[k], ub[k], init = inits[k], integer = integer[k])
    end
    if length(c) > 0
        set_objective!(dict_model, x -> dot(c, flatten(x, ks)[1]))
    end
    if length(bineq) > 0
        add_ineq_constraint!(dict_model, x -> Aineq * flatten(x, ks)[1] - bineq)
    end
    if length(beq) > 0
        add_eq_constraint!(dict_model, x -> Aeq * flatten(x, ks)[1] - beq)
    end
    return dict_model
end

ext: jl


text:
module Sloth

using Knet
using MAT
using Images, FileIO
import Knet.Ops20: BNMoments
import Knet: atype

const SlothArray = Union{Array{T}, KnetArray{T}} where T <: AbstractFloat
const SlothParam = Param{T} where T <: SlothArray
const SlothWeight = Union{SlothArray, SlothParam}
const SlothBias = Union{SlothArray, SlothParam, AbstractFloat}
const IntHyperparam = Union{Int, Tuple{Vararg{Int}}}
F(x::T) where T <: AbstractFloat = eltype(atype())(x)
dir() = abspath(@__DIR__)
dir(args...) = abspath(joinpath(dir(), args...))

include("layers.jl")
include("rnn.jl")
include("beautify.jl")
include("data.jl")
include("vgg.jl")

end # module


text:
export FixedMarginalConstraint

"""
FixedMarginalConstraint

One of the form constraint objects.
Provides a constraint on the marginal distribution such that it remains fixed during inference.
Can be viewed as blocking of updates of a specific edge associated with the marginal.

See also: [`constrain_form`](@ref), [`DistProduct`](@ref)
"""

mutable struct FixedMarginalConstraint <: ReactiveMP.AbstractFormConstraint
    fixed_value :: Any
end

default_form_check_strategy(::FixedMarginalConstraint) = FormConstraintCheckLast()

is_point_mass_form_constraint(::FixedMarginalConstraint) = false

function constrain_form(constraint::FixedMarginalConstraint, something)
    if constraint.fixed_value !== nothing
        return Message(constraint.fixed_value, false, false)
    else
        return something
    end
end


text:
function create_structural_model_input(contracts::DataFrame, spot::DataFrame,
                                            structural_parameters::StructuralParameters)

    exogenous_variables       = handle_f0(spot)
    parameters                = handle_parameters(structural_parameters,
                                    exogenous_variables["f0_spot"])
    vec_contracts_by_day      = read_ContractsByDay(contracts)
    dict_parametric_functions = create_dict_parametric_functions()
    maximum_maturity          = get_maximum_maturity(parameters)

    return create_structural_model_input_struct(parameters,
                vec_contracts_by_day,
                exogenous_variables,
                dict_parametric_functions,
                maximum_maturity)
    end

function create_structural_model_input_struct(parameters::Dict,
          vec_contracts_by_day::Vector,
          exogenous_variables::Dict,
          dict_parametric_functions::Dict{String, Vector{Function}},
          maximum_maturity::Int64)

    timestamps          = date_of_trades(vec_contracts_by_day)
    exogenous_variables = handle_structural_model_structure(exogenous_variables, parameters, timestamps)
    X, X_var_names      = create_structural_model_X(exogenous_variables,
                    maximum_maturity,
                    dict_parametric_functions,
                    has_f0_spot = true)

    return StructuralModelInputs(parameters,
                        vec_contracts_by_day,
                        X), X_var_names

end

function handle_f0(spot::DataFrame)
    name       = "f0_spot"
    timestamps = convert(Vector{DateTime}, spot.timestamps)
    vals       = convert(Vector{Float64}, spot.vals)

    time_series = TimeSeriesInterface.TimeSeries(name,
                timestamps,
                vals)

    return Dict{String, TimeSeries{Float64}}("f0_spot" => time_series)
end

function handle_parameters(structural_parameters::StructuralParameters,
                                f0_spot::TimeSeries{Float64})

    parameters = Dict{String, Any}()
    parameters["interest_yield"]             = structural_parameters.interest_yield
    parameters["intercept"]                  = structural_parameters.intercept
    parameters["intercept_varying_maturity"] = structural_parameters.intercept_varying_maturity
    parameters["seasonality_trading_date"]   = structural_parameters.seasonality_trading_date
    parameters["seasonality_delivery_date"]  = structural_parameters.seasonality_delivery_date
    parameters["run_lasso"]                  = structural_parameters.run_lasso
    parameters["minimum_maturity"]           = structural_parameters.minimum_maturity
    parameters["λ"]                          = structural_parameters.λ
    parameters["f0_spot"]                    = f0_spot

    return parameters
end

function create_dict_parametric_functions()
    dict_parametric_functions = Dict{String,Vector{Function}}()
    dict_parametric_functions["seasonality_delivery_date"] = Function[
                                                                f_sin_t_j_1_harm,
                                                                f_cos_t_j_1_harm,
                                                                #f_sin_t_j_2_harm,
                                                                #f_cos_t_j_2_harm
                                                                ]
    dict_parametric_functions["seasonality_trading_date"] = Function[
                                                                f_sin_t_1_harm,
                                                                f_cos_t_1_harm,
                                                                #f_sin_t_2_harm,
                                                                #f_cos_t_2_harm
                                                                ]
    dict_parametric_functions["intercept_varying_maturity"] = Function[
                                                                f_sqrt_j_1,f_sqrt_j_2,f_sqrt_j_3
                                                                ]
    dict_parametric_functions["spot"] = Function[
                                            f_sqrt_j_0, f_sqrt_j_1,f_sqrt_j_2
                                            ]
    return dict_parametric_functions
end

function get_maximum_maturity(parameters::Dict)
    return parameters["minimum_maturity"]
end

function handle_structural_model_structure(exogenous_variables::Dict,
       parameters::Dict,
       timestamps::Vector)
    ones_aux_ts = TimeSeries(
                        "ones aux ts",
                        timestamps,
                        ones(length(timestamps))
                        )
    if parameters["seasonality_delivery_date"]
        exogenous_variables["seasonality_delivery_date"] = ones_aux_ts
    end
    if parameters["seasonality_trading_date"]
        exogenous_variables["seasonality_trading_date"] = ones_aux_ts
    end
    if parameters["intercept_varying_maturity"]
        exogenous_variables["intercept_varying_maturity"] = ones_aux_ts
    end
    if parameters["intercept"]
        exogenous_variables["intercept"] = ones_aux_ts
    end
    return exogenous_variables
end

function create_structural_model_X(structural_exogenous::Dict,
                                        maximum_maturity::Int,
                                        dict_parametric_functions::Dict;
                                        has_f0_spot = true)
    names_exogenous = get_alphabetically_ordered_keys(structural_exogenous)
    names_exogenous = handle_exogenous_names(names_exogenous)
    X = nothing
    X_columns_names = String[]
    for (n, exogenous) in enumerate(names_exogenous)
        series = structural_exogenous[exogenous]
        parametric_functions = handle_exogenous_parametric_functions(
                                        exogenous,
                                        dict_parametric_functions
                                        )
        x = create_exogenous(
                        series, maximum_maturity,
                        parametric_functions,
                        has_f0_spot = has_f0_spot
                        )
        x_columns_names = exogenous_name(
                                exogenous, parametric_functions
                                )
        X = append_X(X, x)
        X_columns_names = vcat(X_columns_names, x_columns_names)
    end
    return X, X_columns_names
end

function exogenous_name(exogenous::String, parametric_function::Function)
    if haskey(DICT_FUNCTIONS_NAMES, parametric_function)
        return exogenous * DICT_FUNCTIONS_NAMES[parametric_function]
    end
    return exogenous * "_customized_function"
end

function exogenous_name(exogenous::String, parametric_functions::Vector)
    names = String[]
    for (f, func) in enumerate(parametric_functions)
        push!(names, exogenous_name(exogenous, func))
    end
    return names
end

function get_alphabetically_ordered_keys(dict::Dict)
    ks = String[]
    for k in keys(dict)
        push!(ks, k)
    end
    return sort(ks)
end

function handle_exogenous_names(names_exogenous)
    names_exogenous = remove(names_exogenous, "f0_spot")
    return names_exogenous
end

function handle_exogenous_parametric_functions(exogenous::String,
           dict_parametric_functions::Dict)
    parametric_functions = Function[f_sqrt_j_0]
    if exogenous == "spot"
        parametric_functions = dict_parametric_functions["spot"]
    end
    if exogenous == "seasonality_delivery_date"
        parametric_functions = dict_parametric_functions["seasonality_delivery_date"]
    end
    if exogenous == "seasonality_trading_date"
        parametric_functions = dict_parametric_functions["seasonality_trading_date"]
    end
    if exogenous == "intercept_varying_maturity"
        parametric_functions = dict_parametric_functions["intercept_varying_maturity"]
    end
    return parametric_functions
end

function create_exogenous(series::TimeSeries,
                            maximum_maturity::Int,
                            parametric_function::Function;
                            has_f0_spot = true)

    T = length(series.timestamps)
    J = maximum_maturity
    x = Matrix{Float64}(undef, T, J + has_f0_spot)
    for t in 1:T
        for j in 1:J
            x[t,j] = series.vals[t] * parametric_function(series.timestamps[t], j)
        end
        if has_f0_spot
            x[t,J + 1] = series.vals[t] * parametric_function(series.timestamps[t], 0)
        end
    end

    return x
end

function create_exogenous(series::TimeSeries,
                            maximum_maturity::Int,
                            parametric_functions::Vector{Function};
                            has_f0_spot=false)
    T = length(series.timestamps)
    J = maximum_maturity
    num_coefs = length(parametric_functions)
    x = Array{Float64}(undef, T, J + has_f0_spot, num_coefs)

    for (f, func) in enumerate(parametric_functions)
        x[:,:, f] = create_exogenous(
                            series,
                            maximum_maturity,
                            func,
                            has_f0_spot=has_f0_spot
                            )
    end
    return x
end

function handle_exogenous_names(names_exogenous)
    names_exogenous = remove(names_exogenous, "f0_spot")
    return names_exogenous
end

function remove(vec_string::Vector{String}, string)
    if string in vec_string
        vec_string = setdiff(
        vec_string, vec_string[findall(x -> x == "f0_spot", vec_string)]
        )
    end
    return vec_string
end

function append_X(X,
                x::Array)
    if X == nothing
        return x
    end
    return cat(X, x, dims=3)
end

# Parametric Functions
f_sqrt_j_0(t::DateTime, j::Int) = sqrt(j)^0
f_sqrt_j_1(t::DateTime, j::Int) = sqrt(j)^1
f_sqrt_j_2(t::DateTime, j::Int) = sqrt(j)^2
f_sqrt_j_3(t::DateTime, j::Int) = sqrt(j)^3

f_sin_t_1_harm(t::DateTime, j::Int) = sin(2 * π * 1 * dayofyear(t) / 365)
f_cos_t_1_harm(t::DateTime, j::Int) = cos(2 * π * 1 * dayofyear(t) / 365)
#f_sin_t_2_harm(t::DateTime, j::Int) = sin(2 * π * 2 * dayofyear(t) / 365)
#f_cos_t_2_harm(t::DateTime, j::Int) = cos(2 * π * 2 * dayofyear(t) / 365)

f_sin_t_j_1_harm(t::DateTime, j::Int) = sin(2 * π * 1 * dayofyear(t + Day(j)) / 365)
f_cos_t_j_1_harm(t::DateTime, j::Int) = cos(2 * π * 1 * dayofyear(t + Day(j)) / 365)
#f_sin_t_j_2_harm(t::DateTime, j::Int) = sin(2 * π * 2 * dayofyear(t + Day(j)) / 365)
#f_cos_t_j_2_harm(t::DateTime, j::Int) = cos(2 * π * 2 * dayofyear(t + Day(j)) / 365)

DICT_FUNCTIONS_NAMES = Dict(
f_sqrt_j_0 => "_√j⁰",
f_sqrt_j_1 => "_√j¹",
f_sqrt_j_2 => "_√j²",
f_sqrt_j_3 => "_√j³",
f_sin_t_j_1_harm => "_sin_1_harmonic",
#f_sin_t_j_2_harm => "_sin_2_harmonic",
f_cos_t_j_1_harm => "_cos_1_harmonic",
#f_cos_t_j_2_harm => "_cos_2_harmonic",
f_sin_t_1_harm => "_sin_1_harmonic",
#f_sin_t_2_harm => "_sin_2_harmonic",
f_cos_t_1_harm => "_cos_1_harmonic",
#f_cos_t_2_harm => "_cos_2_harmonic"
)

function handle_data_frame_of_reconstructed_prices(contracts::DataFrame, reconstructed_forward_prices::Vector{ContractsByDay})

    reconstructed_contracts = deepcopy(Matrix(contracts))
    number_maturities       = length(reconstructed_forward_prices[1].contracts)
    number_observations     = length(reconstructed_forward_prices)
    number_cols_contracts   = size(contracts, 2)

    @assert number_observations == size(contracts, 1)

    for col in 2:3:number_cols_contracts - 2
        reconstructed_contracts[:, col] .= missing
    end

    for (t, contract_info) in enumerate(reconstructed_forward_prices)
        cont_maturity = 1
        for col in 2:3:number_cols_contracts - 2
            reconstructed_contracts[t, col] = contract_info.contracts[cont_maturity].price
            cont_maturity += 1
        end
    end

    return DataFrame(reconstructed_contracts, names(contracts))
end

max_stars_repo_name: dietzemarina/ForwardCurveSmoother

max_issues_repo_name: dietzemarina/ForwardCurveSmoother

max_forks_repo_name: dietzemarina/ForwardCurveSmoother


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using BenchmarkTools, ArraysOfArrays

inp = rand(5, 5, 3)
@show nestedview(inp, 2) == [inp[:, :, k] for k in axes(inp, 3)]
print("Comprehension:          ")
@btime out = [$inp[:, :, k] for k in axes($inp, 3)]
print("Comprehension with @view:")
@btime out = [@view($inp[:, :, k]) for k in axes($inp, 3)]
print("ArraysOfArrays.nestedview:")
@btime out = nestedview($inp, 2);

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
const example_solution_short = Solution(
    [-1.5, -1.2, -1.3],
    [[1, 2], [0, 1], [3, 4]],
    [0.25, 0.1, 0.2],
    [1, 1, 1],
    0.666
)

const example_degenerate_solution = Solution(
    [-1.5, -1.2, -1.5, -1.6, -1.1],
    [[1, 2, 3], [1, 3, 1], [3, 4, 1], [1, 4, 2], [4, 2, 1]],
    [0.23, 0.1, 0.23, 0.25, 0.05],
    [1, 1, 1, 1, 1],
    0.22
)

@testset "Bounding solution of size ≥ max_states gives solution of max_states length" begin
    bounded = bound_solution(example_solution_short, 2)
    @test length(bounded.energies) == 2
    @test length(bounded.states) == 2
    @test length(bounded.probabilities) == 2
    @test length(bounded.degeneracy) == 2

end

@testset "Bounding solution of size < max_states gives the same solution (up to permutation)" begin
    bounded = bound_solution(example_solution_short, 10)
    @test bounded.energies == [-1.5, -1.3, -1.2]
    @test bounded.probabilities == [0.25, 0.2, 0.1]
    @test bounded.states == [[1, 2], [3, 4], [0, 1]]
    @test bounded.degeneracy == [1, 1, 1]
end

@testset "Bounding solution clips at correct probability" begin
    @testset "when max_states=$(max_states)" for (max_states, expected_prob) ∈ [(1, 0.23), (2, 0.23), (3, 0.22)]
        bounded = bound_solution(example_degenerate_solution, max_states)
        @test bounded.largest_discarded_probability == expected_prob
        @test all(bounded.probabilities .≥ expected_prob)
    end
end

@testset "Bounding solution preserves correspondence between energies, states and probabilities" begin
    sorted_idx = sortperm(example_degenerate_solution.probabilities, rev=true)
    for max_states ∈ 1:5
        bounded = bound_solution(example_degenerate_solution, max_states)
        @test bounded.energies == example_degenerate_solution.energies[sorted_idx][1:max_states]
        @test bounded.probabilities == example_degenerate_solution.probabilities[sorted_idx][1:max_states]
        @test bounded.states == example_degenerate_solution.states[sorted_idx][1:max_states]
    end
end


text:
import Base.getindex

@inline getindex{A<:AbstractArray, IM1<:AbstractIndexManipulator, IM2<:AbstractIndexManipulator}(arr::A, ::Type{IM1}, ::Type{IM2}, x::Number...) = getindex(arr, Tuple{IM1}, IM2, x...)

@inline inplaceadd!{A<:AbstractArray, IM1<:AbstractIndexManipulator, IM2<:AbstractIndexManipulator}(arr::A, ::Type{IM1}, ::Type{IM2}, v::Number, x::Number...) = inplaceadd!(arr, Tuple{IM1}, IM2, v, x...)

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


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

## floating-point functions ##

copysign(x::Float64, y::Float64) = copysign_float(x, y)
copysign(x::Float32, y::Float32) = copysign_float(x, y)
copysign(x::Float32, y::Real) = copysign(x, Float32(y))
copysign(x::Float64, y::Real) = copysign(x, Float64(y))

flipsign(x::Float64, y::Float64) = bitcast(Float64, xor_int(bitcast(UInt64, x), and_int(bitcast(UInt64, y), 0x8000000000000000)))
flipsign(x::Float32, y::Float32) = bitcast(Float32, xor_int(bitcast(UInt32, x), and_int(bitcast(UInt32, y), 0x80000000)))
flipsign(x::Float32, y::Real) = flipsign(x, Float32(y))
flipsign(x::Float64, y::Real) = flipsign(x, Float64(y))

signbit(x::Float64) = signbit(bitcast(Int64, x))
signbit(x::Float32) = signbit(bitcast(Int32, x))
signbit(x::Float16) = signbit(bitcast(Int16, x))

"""
    maxintfloat(T)

The largest integer losslessly representable by the given floating-point DataType `T`.
"""
maxintfloat(::Type{Float64}) = 9007199254740992.
maxintfloat(::Type{Float32}) = Float32(16777216.)
maxintfloat(::Type{Float16}) = Float16(2048f0)
maxintfloat(x::T) where {T<:AbstractFloat} = maxintfloat(T)

"""
    maxintfloat(T, S)

The largest integer losslessly representable by the given floating-point DataType `T` that
also does not exceed the maximum integer representable by the integer DataType `S`.
"""
maxintfloat(::Type{S}, ::Type{T}) where {S<:AbstractFloat, T<:Integer} = min(maxintfloat(S), S(typemax(T)))
maxintfloat() = maxintfloat(Float64)

isinteger(x::AbstractFloat) = (x - trunc(x) == 0)

num2hex(x::Float16) = hex(bitcast(UInt16, x), 4)
num2hex(x::Float32) = hex(bitcast(UInt32, x), 8)
num2hex(x::Float64) = hex(bitcast(UInt64, x), 16)

function hex2num(s::AbstractString)
    if length(s) <= 4
        return bitcast(Float16, parse(UInt16, s, 16))
    end
    if length(s) <= 8
        return bitcast(Float32, parse(UInt32, s, 16))
    end
    return bitcast(Float64, parse(UInt64, s, 16))
end

"""
    round([T,] x, [digits, [base]], [r::RoundingMode])

Rounds `x` to an integer value according to the provided
[`RoundingMode`](@ref), returning a value of the same type as `x`. When not
specifying a rounding mode the global mode will be used
(see [`rounding`](@ref)), which by default is round to the nearest integer
([`RoundNearest`](@ref) mode), with ties (fractional values of 0.5) being
rounded to the nearest even integer.

# Examples
```jldoctest
julia> round(1.7)
2.0

julia> round(1.5)
2.0

julia> round(2.5)
2.0
```

The optional [`RoundingMode`](@ref) argument will change how the number gets
rounded.

`round(T, x, [r::RoundingMode])` converts the result to type `T`, throwing an
[`InexactError`](@ref) if the value is not representable.

`round(x, digits)` rounds to the specified number of digits after the decimal place (or
before if negative). `round(x, digits, base)` rounds using a base other than 10.

# Examples
```jldoctest
julia> round(pi, 2)
3.14

julia> round(pi, 3, 2)
3.125
```

!!! note
    Rounding to specified digits in bases other than 2 can be inexact when
    operating on binary floating point numbers. For example, the [`Float64`](@ref)
    value represented by `1.15` is actually *less* than 1.15, yet will be
    rounded to 1.2.

    # Examples
    ```jldoctest
    julia> x = 1.15
    1.15

    julia> @sprintf "%.20f" x
    "1.14999999999999991118"

    julia> x < 115//100
    true

    julia> round(x, 1)
    1.2
    ```
"""
round(T::Type, x)
round(x::Real, ::RoundingMode{:ToZero}) = trunc(x)
round(x::Real, ::RoundingMode{:Up}) = ceil(x)
round(x::Real, ::RoundingMode{:Down}) = floor(x)
# C-style round
function round(x::AbstractFloat, ::RoundingMode{:NearestTiesAway})
    y = trunc(x)
    ifelse(x==y,y,trunc(2*x-y))
end
# Java-style round
function round(x::AbstractFloat, ::RoundingMode{:NearestTiesUp})
    y = floor(x)
    ifelse(x==y,y,copysign(floor(2*x-y),x))
end
round(::Type{T}, x::AbstractFloat, r::RoundingMode) where {T<:Integer} = trunc(T,round(x,r))

# adapted from Matlab File Exchange roundsd: http://www.mathworks.com/matlabcentral/fileexchange/26212
# for round, og is the power of 10 relative to the decimal point
# for signif, og is the absolute power of 10
# digits and base must be integers, x must be convertable to float

function _signif_og(x, digits, base)
    if base == 10
        e = floor(log10(abs(x)) - digits + 1.)
        og = oftype(x, exp10(abs(e)))
    elseif base == 2
        e = exponent(abs(x)) - digits + 1.
        og = oftype(x, exp2(abs(e)))
    else
        e = floor(log(base, abs(x)) - digits + 1.)
        og = oftype(x, float(base) ^ abs(e))
    end
    return og, e
end

"""
    signif(x, digits, [base])

Rounds (in the sense of [`round`](@ref)) `x` so that there are `digits` significant digits, under a
base `base` representation, default 10.

# Examples
```jldoctest
julia> signif(123.456, 2)
120.0

julia> signif(357.913, 4, 2)
352.0
```
"""
function signif(x::Real, digits::Integer, base::Integer=10)
    digits < 1 && throw(DomainError(digits, "`digits` cannot be less than 1."))

    x = float(x)
    (x == 0 || !isfinite(x)) && return x
    og, e = _signif_og(x, digits, base)
    if e >= 0 # for numeric stability
        r = round(x/og)*og
    else
        r = round(x*og)/og
    end
    !isfinite(r) ? x : r
end

for f in (:round, :ceil, :floor, :trunc)
    @eval begin
        function ($f)(x::Real, digits::Integer, base::Integer=10)
            x = float(x)
            og = convert(eltype(x),base)^digits
            r = ($f)(x * og) / og

            if !isfinite(r)
                if digits > 0
                    return x
                elseif x > 0
                    return $(:ceil == f ? :(convert(eltype(x), Inf)) : :(zero(x)))
                elseif x < 0
                    return $(:floor == f ? :(-convert(eltype(x), Inf)) : :(-zero(x)))
                else
                    return x
                end
            end
            return r
        end
    end
end

# isapprox: approximate equality of numbers
"""
    isapprox(x, y; rtol::Real=sqrt(eps), atol::Real=0, nans::Bool=false, norm::Function)

Inexact equality comparison: `true` if `norm(x-y) <= atol + rtol*max(norm(x), norm(y))`. The
default `atol` is zero and the default `rtol` depends on the types of `x` and `y`. The keyword
argument `nans` determines whether or not NaN values are considered equal (defaults to false).

For real or complex floating-point values, `rtol` defaults to
`sqrt(eps(typeof(real(x-y))))`. This corresponds to requiring equality of about half of the
significand digits. For other types, `rtol` defaults to zero.

`x` and `y` may also be arrays of numbers, in which case `norm` defaults to `vecnorm` but
may be changed by passing a `norm::Function` keyword argument. (For numbers, `norm` is the
same thing as `abs`.) When `x` and `y` are arrays, if `norm(x-y)` is not finite (i.e. `±Inf`
or `NaN`), the comparison falls back to checking whether all elements of `x` and `y` are
approximately equal component-wise.

The binary operator `≈` is equivalent to `isapprox` with the default arguments, and `x ≉ y`
is equivalent to `!isapprox(x,y)`.

# Examples
```jldoctest
julia> 0.1 ≈ (0.1 - 1e-10)
true

julia> isapprox(10, 11; atol = 2)
true

julia> isapprox([10.0^9, 1.0], [10.0^9, 2.0])
true
```
"""
function isapprox(x::Number, y::Number; rtol::Real=rtoldefault(x,y), atol::Real=0, nans::Bool=false)
    x == y || (isfinite(x) && isfinite(y) && abs(x-y) <= atol + rtol*max(abs(x), abs(y))) || (nans && isnan(x) && isnan(y))
end

const ≈ = isapprox
≉(args...; kws...) = !≈(args...; kws...)

# default tolerance arguments
rtoldefault(::Type{T}) where {T<:AbstractFloat} = sqrt(eps(T))
rtoldefault(::Type{<:Real}) = 0
rtoldefault(x::Union{T,Type{T}}, y::Union{S,Type{S}}) where {T<:Number,S<:Number} = max(rtoldefault(real(T)),rtoldefault(real(S)))

# fused multiply-add
fma_libm(x::Float32, y::Float32, z::Float32) =
    ccall(("fmaf", libm_name), Float32, (Float32,Float32,Float32), x, y, z)
fma_libm(x::Float64, y::Float64, z::Float64) =
    ccall(("fma", libm_name), Float64, (Float64,Float64,Float64), x, y, z)
fma_llvm(x::Float32, y::Float32, z::Float32) = fma_float(x, y, z)
fma_llvm(x::Float64, y::Float64, z::Float64) = fma_float(x, y, z)
# Disable LLVM's fma if it is incorrect, e.g. because LLVM falls back
# onto a broken system libm; if so, use openlibm's fma instead
# 1.0000305f0 = 1 + 1/2^15
# 1.0000000009313226 = 1 + 1/2^30
# If fma_llvm() clobbers the rounding mode, the result of 0.1 + 0.2 will be 0.3
# instead of the properly-rounded 0.30000000000000004; check after calling fma
if (Sys.ARCH != :i686 && fma_llvm(1.0000305f0, 1.0000305f0, -1.0f0) == 6.103609f-5 &&
    (fma_llvm(1.0000000009313226, 1.0000000009313226, -1.0) ==
     1.8626451500983188e-9) && 0.1 + 0.2 == 0.30000000000000004)
    fma(x::Float32, y::Float32, z::Float32) = fma_llvm(x,y,z)
    fma(x::Float64, y::Float64, z::Float64) = fma_llvm(x,y,z)
else
    fma(x::Float32, y::Float32, z::Float32) = fma_libm(x,y,z)
    fma(x::Float64, y::Float64, z::Float64) = fma_libm(x,y,z)
end
function fma(a::Float16, b::Float16, c::Float16)
    Float16(fma(Float32(a), Float32(b), Float32(c)))
end

# This is necessary at least on 32-bit Intel Linux, since fma_llvm may
# have called glibc, and some broken glibc fma implementations don't
# properly restore the rounding mode
Rounding.setrounding_raw(Float32, Rounding.JL_FE_TONEAREST)
Rounding.setrounding_raw(Float64, Rounding.JL_FE_TONEAREST)

max_stars_repo_name: jdlangs/julia

max_issues_repo_name: jdlangs/julia

max_forks_repo_name: jdlangs/julia


text:
"""
    regulaFalsi(f,a,b)

adalah fungsi yang mencari nilai akar persamaan dari fungsi `f(x)` pada interval
di antara `a` dan `b`.

# Examples
```jldoctest
julia> f(x) = x*sin(x)-1;

julia> c, flag, M = regulaFalsi(f,0,2);

julia> c
1.1141571430336825

julia> flag
0

julia> M
4×5 Array{Float64,2}:
 0.0  0.0      1.09975  2.0      -0.0200192
 1.0  1.09975  1.12124  2.0       0.00983461
 2.0  1.09975  1.11416  1.12124   5.63036e-6
 3.0  1.09975  1.11416  1.11416   3.00226e-9
```
return solusi `c` dengan `flag` bernilai 0 jika metode regula falsi berhasil menemukan
solusi dan gagal jika tidak 0. Serta, matriks `M` yang berisi catatan proses tiap
iterasi `[k, a, c, b, f(c)]`

"""
function regulaFalsi(f,a,b)
    delta = 10^-7; maxi = 100; flag = 1;
    M = Array{Float64}(undef, 0, 5);
    fa = f(a); fb = f(b);
    if fa*fb > 0
        c = "error : fa fb harus beda tanda";
        flag = 2;
        return;
    end
    for k = 1:maxi
        c = b-fb*(b-a)/(fb-fa);
        fc = f(c);
        dx = min(c-a,b-c);
        M = [M ; [k-1 a c b fc] ];
        if fc == 0
            a = c;
            b = c;
        elseif fa*fc>0
            a = c;
            fa= fc;
        else
            b = c;
            fb= fc;
        end
        if abs(fc) < delta || abs(dx)< delta
            flag = 0; break;
        end
    end
    return c,flag,M
end


text:
module DataFitting

using Printf, PrettyTables
using Statistics, Distributions
using DataStructures
using LsqFit

import Base.push!
import Base.show
import Base.ndims
import Base.size
import Base.length
import Base.getindex
import Base.reshape
import Base.propertynames
import Base.getproperty
import Base.getindex
import Base.setindex!
import Base.iterate

export Domain, CartesianDomain, Measures,
    Prediction, Model, evaluate, parindex, thaw, freeze, fit!

include("domain.jl")

# ====================================================================
# Parameter
#
mutable struct Parameter
    val::Float64
    low::Float64              # lower limit value
    high::Float64             # upper limit value
    step::Float64
    free::Bool
    Parameter(value::Number) = new(float(value), -Inf, +Inf, NaN, true)
end

# ====================================================================
# A *component* is a generic implementation of a constituent part of a
# model.
#
# A component must inherit `AbstractComponent`, and implement the
# `ceval_data` and `evaluate` methods
abstract type AbstractComponent end

function getparams(comp::AbstractComponent)
    params = OrderedDict{Tuple{Symbol,Int}, Parameter}()
    for pname in fieldnames(typeof(comp))
        par = getfield(comp, pname)
        if isa(par, Parameter)
            params[(pname, 0)] = par
        elseif isa(par, Vector{Parameter})
            for i in 1:length(par)
                params[(pname, i)] = par[i]
            end
        end
    end
    params
end

# ====================================================================
# CompEval: a wrapper for a component evaluated on a specific domain
#
mutable struct CompEval{TDomain <: AbstractDomain, TComp <: AbstractComponent}
    domain::TDomain
    comp::TComp
    params::OrderedDict{Tuple{Symbol,Int}, Parameter}
    cdata
    counter::Int
    lastvalues::Vector{Float64}
    eval::Vector{Float64}
    ipar::Vector{Int}  # handled by Model

    function CompEval(domain::AbstractDomain, comp::AbstractComponent)
        params = getparams(comp)
        (cdata, len) = ceval_data(domain, comp)
        return new{typeof(domain), typeof(comp)}(
            domain, comp, params, cdata, 0,
            fill(NaN, length(params)),
            fill(NaN, len), Vector{Int}())
    end
end

# This is called to `update` to distinguish it from component's `evaluate`.
update(c::CompEval) = update(c, [par.val for par in values(c.params)])
function update(c::CompEval, pvalues::Vector{Float64})
    @assert length(c.params) == length(pvalues)

    # Do we actually need a new evaluation?
    if (any(c.lastvalues .!= pvalues)  ||  (c.counter == 0))
        c.lastvalues .= pvalues
        c.counter += 1
        @assert all(.!isnan.(pvalues))
        evaluate(c, pvalues...)
    end
    return c.eval
end

# ====================================================================
# Component fall back methods
ceval_data(domain::AbstractDomain, comp::AbstractComponent) =
    error("Component " * string(typeof(comp)) * " must implement its own method for `ceval_data`.")

evaluate(c::CompEval{TDomain, TComp}, args...) where {TDomain, TComp} =
    error("Component " * string(TComp) * " must implement its own method for `evaluate`.")

# ====================================================================
# Built-in components
#
include("components/ScalarParam.jl")
include("components/FuncWrap.jl")
include("components/OffsetSlope.jl")
include("components/Gaussian.jl")

# ====================================================================
# Parse a user defined structure or dictionary to extract all
# components
function extract_components(things...; prefix="")
    out = OrderedDict{Symbol, AbstractComponent}()
    for thing in things
        #println()
        #println("Thing $(typeof(thing))  (prefix = $(prefix))")
        if isa(thing, AbstractComponent)
            #println("Adding...")
            out[Symbol(prefix)] = thing
        else
            (length(prefix) > 0)  &&  (prefix *= "_")
            if isa(thing, AbstractDict)
                for (name, v) in thing
                    #println("Dict: Walk through $name :: $(typeof(v))")
                    merge!(out, extract_components(v; prefix=prefix * string(name)))
                end
            elseif isa(thing, Pair)
                #println("Pair: $(thing[1]), $(typeof(thing[2]))")
                if isa(thing[1], Symbol)
                    name = thing[1]
                    v = thing[2]
                    isa(v, Number)  &&  (v = ScalarParam(v))
                    if isa(v, AbstractComponent)
                        merge!(out, extract_components(v; prefix=prefix * string(name)))
                    end
                end
            elseif isstructtype(typeof(thing))
                for name in fieldnames(typeof(thing))
                    v = getfield(thing, name)
                    #println("Structure: Walk through $name :: $(typeof(v))")
                    merge!(out, extract_components(v; prefix=prefix * string(name)))
                end
            end
        end
    end
    return out
end

# ====================================================================
# A model prediction suitable to be compared to experimental data
mutable struct Prediction
    domain::AbstractDomain
    cevals::OrderedDict{Symbol, CompEval}
    eval::Vector{Float64}
    reduce_with_dict::Bool
    reducer::Union{Nothing, Function}
    counter::Int

    function Prediction(domain::AbstractDomain, things...;
                        prefix="", reduce=reduce)
        comps = extract_components(things...; prefix=prefix)
        cevals = OrderedDict{Symbol, CompEval}()
        for (name, comp) in comps
            cevals[name] = CompEval(domain, comp)
        end
        out = new(domain, cevals, Vector{Float64}(), false, reduce, 0)
        evaluate(out)  # TODO: is this correct?
        return out
    end
end

# Default reducer: add all components
reduce(domain::AbstractDomain, args...) = .+(args...)

# Reduce prediction by combining individual components
function reduce(pred::Prediction)
    if pred.reduce_with_dict
        d = Dict([(cname, ceval.eval) for (cname, ceval) in pred.cevals])
        expr = pred.reducer(pred.domain, d)
    else
        d = [ceval.eval for (cname, ceval) in pred.cevals]
        expr = pred.reducer(pred.domain, d...)
    end
    if length(pred.eval) == 0
        append!(pred.eval, expr)
    else
        pred.eval .= expr
    end
    pred.counter += 1
end

function evaluate(pred::Prediction)
    for (name, ceval) in pred.cevals
        update(ceval)
    end
    reduce(pred)
end

# ====================================================================
# Global model, actually a collection of `Prediction`s.
mutable struct Model
    preds::Vector{Prediction}
    comps::OrderedDict{Symbol, AbstractComponent}
    cfree::OrderedDict{Symbol, Bool}
    params::OrderedDict{Tuple{Symbol, Symbol, Int}, Parameter}
    pvalues::Vector{Float64}
    actual::Vector{Float64}
    buffer::Vector{Float64}
    partransform::Function
end

function Model(v::Vector{Prediction})
    model = Model(v, OrderedDict{Symbol, AbstractComponent}(),
                  OrderedDict{Symbol, Bool}(),
                  OrderedDict{Tuple{Symbol, Symbol, Int}, Parameter}(),
                  Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), default_partransform)
    evaluate(model)
    return model
end

Model(p::Prediction) = Model([p])
Model(args...; kw...) = Model(Prediction(args...; kw...))

function evaluate(model::Model)
    @assert length(model.preds) >= 1

    # Save list of previously free components
    cfree = deepcopy(model.cfree)

    # Collect components and parameters
    empty!(model.comps)
    empty!(model.cfree)
    empty!(model.params)
    for pred in model.preds
        for (cname, ceval) in pred.cevals
            model.comps[cname] = ceval.comp
            model.cfree[cname] = get(cfree, cname, true)
            for (pname, par) in ceval.params
                cpname = (cname, pname[1], pname[2])
                model.params[cpname] = par
            end
        end
    end

    # Populate CompEval.ipar and evaluate all predictions
    ndata = 0
    cpnames = keys(model.params)
    for pred in model.preds
        for (cname, ceval) in pred.cevals
            empty!(ceval.ipar)
            for (pname, par) in ceval.params
                cpname = (cname, pname[1], pname[2])
                push!(ceval.ipar, findfirst(cpnames .== Ref(cpname)))
            end
            update(ceval)
        end
        reduce(pred)
        ndata += length(pred.eval)
    end

    model.pvalues = [par.val for par in values(model.params)]
    model.actual = deepcopy(model.pvalues)
    model.buffer = Vector{Float64}(undef, ndata)
    quick_evaluate(model)
    return model
end

default_partransform(model::Model, pvalues::Vector{Float64}, actual::Vector{Float64}) = nothing

# This is supposed to be called from `fit!`, not by user
function quick_evaluate(model::Model)
    model.actual .= model.pvalues  # copy all values by default
    model.partransform(model, model.pvalues, model.actual)

    for pred in model.preds
        for (cname, ceval) in pred.cevals
            update(ceval, model.actual[ceval.ipar])
        end
    end
    for pred in model.preds
        reduce(pred)
    end
    nothing
end

function Base.push!(m::Model, p::Prediction)
    push!(m.preds, p)
    evaluate(model)
    return model
end

Base.getindex(m::Model, i::Int) = m.preds[i].eval
Base.getindex(m::Model, cname::Symbol) = m.comps[cname]

parindex(model::Model, cname::Symbol, pname::Symbol, i::Int=0) =
    findfirst(keys(model.params) .== Ref((cname, pname, i)))

function freeze(model::Model, cname::Symbol)
    @assert cname in keys(model.cfree) "Component $c is not defined"
    model.cfree[cname] = false
    model
end

function thaw(model::Model, cname::Symbol)
    @assert cname in keys(model.cfree) "Component $c is not defined"
    model.cfree[cname] = true
    model
end

# ====================================================================
# Fit results
#
struct BestFitPar
    val::Float64
    unc::Float64
    free::Bool
    calc::Float64  # value after transformation
end

struct BestFitComp
    params::OrderedDict{Symbol, Union{BestFitPar, Vector{BestFitPar}}}
    BestFitComp() = new(OrderedDict{Symbol, Union{BestFitPar, Vector{BestFitPar}}}())
end

Base.propertynames(comp::BestFitComp) = keys(getfield(comp, :params))
Base.getproperty(comp::BestFitComp, p::Symbol) = getfield(comp, :params)[p]
Base.getindex(comp::BestFitComp, p::Symbol) = getfield(comp, :params)[p]
Base.length(comp::BestFitComp) = length(getfield(comp, :params))
Base.iterate(comp::BestFitComp, args...) = iterate(getfield(comp, :params), args...)
Base.setindex!(comp::BestFitComp, x, p::Symbol) = getfield(comp, :params)[p] = x

struct BestFitResult
    comps::OrderedDict{Symbol, BestFitComp}
    ndata::Int
    dof::Int
    cost::Float64
    status::Symbol      #:Optimal, :NonOptimal, :Warn, :Error
    log10testprob::Float64
    elapsed::Float64
end

Base.getindex(res::BestFitResult, cname::Symbol) = res.comps[cname]

# ====================================================================
function data1D(model::Model, data::Vector{T}) where T<:AbstractMeasures
    out = Vector{Measures_1D}()
    for i in 1:length(model.preds)
        pred = model.preds[i]
        @assert(length(data[i]) == length(pred.eval),
                "Length of dataset $i do not match corresponding model prediction.")
        push!(out, flatten(data[i], pred.domain))
    end
    return out
end

function residuals1d(model::Model, data1d::Vector{Measures_1D})
    c1 = 1
    for i in 1:length(model.preds)
        pred = model.preds[i]
        c2 = c1 + length(pred.eval) - 1
        model.buffer[c1:c2] .= ((pred.eval .- data1d[i].val) ./ data1d[i].unc)
        c1 = c2 + 1
    end
    return model.buffer
end

# ====================================================================
abstract type AbstractMinimizer end

using LsqFit
mutable struct lsqfit <: AbstractMinimizer
end

function minimize(minimizer::lsqfit, func::Function, params::Vector{Parameter})
    ndata = length(func(getfield.(params, :val)))
    bestfit = LsqFit.curve_fit((dummy, pvalues) -> func(pvalues),
                               1.:ndata, fill(0., ndata),
                               getfield.(params, :val),
                               lower=getfield.(params, :low),
                               upper=getfield.(params, :high))
    status = :NonOptimal
    (bestfit.converged)  &&  (status = :Optimal)
    error = LsqFit.margin_error(bestfit, 0.6827)
    return (status, getfield.(Ref(bestfit), :param), error)
end

macro enable_CMPFit()
    return esc(:(
        import DataFitting.minimize;

        mutable struct cmpfit <: DataFitting.AbstractMinimizer;
        config::CMPFit.Config;
        cmpfit() = new(CMPFit.Config());
        end;

        function minimize(minimizer::cmpfit, func::Function, params::Vector{DataFitting.Parameter});
        guess = getfield.(params, :val);
        low   = getfield.(params, :low);
        high  = getfield.(params, :high);
        parinfo = CMPFit.Parinfo(length(guess));
        for i in 1:length(guess);
        llow  = isfinite(low[i])   ?  1  :  0;
        lhigh = isfinite(high[i])  ?  1  :  0;
        parinfo[i].limited = (llow, lhigh);
        parinfo[i].limits  = (low[i], high[i]);
        end;
        bestfit = CMPFit.cmpfit((pvalues) -> func(pvalues),
                                guess, parinfo=parinfo, config=minimizer.config);
        return (:Optimal, getfield.(Ref(bestfit), :param), getfield.(Ref(bestfit), :perror));
        end;
    ))
end

fit!(model::Model, data::T; kw...) where T<:AbstractMeasures =
    fit!(model, [data]; kw...)

function fit!(model::Model, data::Vector{T};
              minimizer=lsqfit()) where T<:AbstractMeasures
    elapsedTime = Base.time_ns()
    evaluate(model)

    free = Vector{Bool}()
    for (cpname, par) in model.params
        push!(free, par.free  &&  model.cfree[cpname[1]])
    end
    ifree = findall(free)
    @assert length(ifree) > 0 "No free parameter in the model"

    # Flatten empirical data
    data1d = data1D(model, data)

    # Evaluate normalized residuals starting from free parameter values
    function pval2resid(pvalues_free::Vector{Float64})
        model.pvalues[ifree] .= pvalues_free  # update parameter values
        quick_evaluate(model)
        return residuals1d(model, data1d)
    end

    (status, best_val, best_unc) = minimize(minimizer, pval2resid,
                                            collect(values(model.params))[ifree])

    model.pvalues[ifree] .= best_val
    setfield!.(values(model.params), :val, model.pvalues)
    uncerts = fill(NaN, length(model.pvalues))
    uncerts[ifree] .= best_unc

    # Prepare output
    quick_evaluate(model)  # ensure best fit values are used
    comps = OrderedDict{Symbol, BestFitComp}()
    for cname in keys(model.comps)
        comps[cname] = BestFitComp()
    end
    i = 1
    for (cpname, par) in model.params
        cname = cpname[1]
        pname = cpname[2]
        parid = cpname[3]
        bfpar = BestFitPar(model.pvalues[i], uncerts[i],
                           (i in ifree), model.actual[i])
        if parid == 0
            comps[cname][pname] = bfpar
        else
            if parid == 1
                comps[cname][pname] = [bfpar]
            else
                push!(comps[cname][pname], bfpar)
            end
        end
        i += 1
    end
    cost = sum(abs2, model.buffer)
    dof = length(model.buffer) - length(ifree)

    result = BestFitResult(comps, length(model.buffer), dof, cost, status,
                           logccdf(Chisq(dof), cost) * log10(exp(1)),
                           float(Base.time_ns() - elapsedTime) / 1.e9)
    return result
end

include("show.jl")

end


text:
export @with

using GeneralizedGenerated

struct TypelevelExpr{T}
    expr::Expr
    TypelevelExpr(expr::Expr) = new{to_type(expr)}(expr)
end

with(m::Module, nt::NamedTuple, ex::Expr) = with(m, nt, TypelevelExpr(ex))

@gg function with(m::Module, nt::NamedTuple{N,T}, ::TypelevelExpr{E}) where {N,T,E}
    ex = from_type(E)
    q = quote end
    for x in N
        xname = QuoteNode(x)
        push!(q.args, :($x = Base.getproperty(nt, $xname)))
    end
    push!(q.args, ex)
    @under_global :m q
end

@gg function with(m::Module, nt1::NamedTuple{N1}, nt2::NamedTuple{N2}, ::TypelevelExpr{E}) where {N1,N2,E}
    ex = from_type(E)
    q = quote end
    for x in N1
        xname = QuoteNode(x)
        push!(q.args, :($x = Base.getproperty(nt1, $xname)))
    end
    for x in N2
        xname = QuoteNode(x)
        push!(q.args, :($x = Base.getproperty(nt2, $xname)))
    end
    push!(q.args, ex)
    @under_global :m q
end

"""
    @with(ctx..., body)

Compute `body` using the context `ctx`. `ctx` is typically a named tuple, but the only requirement is that it supports `Base.getproperty`.

This macro works by creating a type-level representation of `body` and passing this to a "generalized generated" function that dispatches on the types. There's a default for named tuples:

```
julia> nt = (x=2, y=3)
(x = 2, y = 3)

julia> @with nt begin
       x^2 + y
       end
7
```

This can be extended to other types and a different number of arguments. For example,
```
julia> using NestedTuples: with

julia> NestedTuples.with(nt1, nt2, ex) = with(nt1, ex) + with(nt2, ex)

julia> @with nt nt begin
       x^2 + y
       end
14
```
or
```
julia> NestedTuples.with(nt1, v::Vector, ex) = with(nt1, ex) .+ v

julia> @with nt [1,2,3] begin
       x^2 + y
       end
3-element Vector{Int64}:
  8
  9
 10
```
"""
macro with(args...)
    ctx = esc.(args[1:end-1])
    ex = args[end]

    tle = TypelevelExpr(ex)
    quote
        with($__module__, $(ctx...), $tle)
    end
end


text:
using ExpFamily
using Base.Test

srand(123)

C  = [1.0 0.5; 0.5 1.0]
P  = 2/3 * [2 -1; -1 2]
m  = [1.0;-0.5]

gNP = GaussianNatParam(mean=m, cov=C)
gMP = GaussianMeanParam(mean=m, cov=C)

@test isapprox(gNP.theta1, P*m)
@test isapprox(gNP.theta2, -P)
@test isapprox(gMP.mu1, m)
@test isapprox(gMP.mu2, (m*m'+C)/2)

@test isvalid(gNP)
@test isvalid(gMP)

mpFromNP  = meanparam(gNP)
mpFromNP2 = meanparam(gNP, 2.0)
npFromMP  = natparam(gMP)

@test isapprox(mpFromNP.mu1, gMP.mu1)
@test isapprox(mpFromNP.mu2, gMP.mu2)
@test isapprox(npFromMP.theta1, gNP.theta1)
@test isapprox(npFromMP.theta2, gNP.theta2)
@test isapprox(cov(mpFromNP2), 2cov(mpFromNP))

@test isapprox(mean(gNP), m)
@test isapprox(mean(gMP), m)
@test isapprox(cov(gNP), C)
@test isapprox(cov(gMP), C)
@test isapprox(var(gMP), diag(C))
@test isapprox(var(gNP), diag(C))
@test isapprox(std(gMP), sqrt.(diag(C)))
@test isapprox(std(gNP), sqrt.(diag(C)))

@test_throws DomainError GaussianNatParam(mean=m,cov=-C,check=true)
@test_throws DomainError GaussianMeanParam(mean=m,cov=-C,check=true)

dim = 4
gNPones = ones(GaussianNatParam, dim)
gMPones = ones(GaussianMeanParam, dim)

@test isapprox(mean(gNPones), zeros(dim))
@test isapprox(cov(gNPones),  eye(dim))
@test isapprox(mean(gMPones), zeros(dim))
@test isapprox(cov(gMPones),  eye(dim))

@test length(gNPones) == dim
@test length(gMPones) == dim

rNP = rand(gNPones,5000)
rMP = rand(gMPones,5000)

@test isapprox(mean(rNP,2), zeros(dim), atol=1e-1)
@test isapprox(mean(rMP,2), zeros(dim), atol=1e-1)
@test isapprox(cov(rNP,2),  eye(dim),   atol=1e-1)
@test isapprox(cov(rMP,2),  eye(dim),   atol=1e-1)

gNP2 = gNP+gNP
gMP2 = gMP+gMP

@test isapprox(gNP2.theta1, 2gNP.theta1)
@test isapprox(gNP2.theta2.data, 2gNP.theta2.data)
@test isapprox(gMP2.mu1, 2gMP.mu1)
@test isapprox(gMP2.mu2.data, 2gMP.mu2.data)

m1  = randn(dim)
m2  = randn(dim)
C1  = randn(dim,dim)
C1 *= C1'
C2  = randn(dim,dim)
C2 *= C2'

gNPa = GaussianNatParam(mean=m1, cov=C1)
gNPb = GaussianNatParam(mean=m2, cov=C2)
gMPa = GaussianMeanParam(mean=m1, cov=C1)
gMPb = GaussianMeanParam(mean=m2, cov=C2)

gNPtest = gNPa + gNPb - natparam(gMPb)
gMPtest = gMPa + gMPb - meanparam(gNPb)

@test isapprox(gNPtest.theta1, gNPa.theta1)
@test isapprox(gNPtest.theta2, gNPa.theta2)
@test isapprox(gMPtest.mu1, gMPa.mu1)
@test isapprox(gMPtest.mu2, gMPa.mu2)

@test isapprox(gMPtest, gMPa)

@test isapprox(GaussianNatParam(vec(gNPa),dim),gNPa)
@test isapprox(GaussianMeanParam(vec(gMPa),dim),gMPa)

@test isapprox(2gNPa + gNPa*2 - 3gNPa, gNPa)
@test isapprox(2gMPa + gMPa*2 - 3gMPa, gMPa)
@test isapprox(gNPa/2 + gNPa/2, gNPa)
@test isapprox(gMPa/2 + gMPa/2, gMPa)

gNP = GaussianNatParam(mean=m, cov=C)
@test isapprox(project(gNP), gNP)

nP = -diagm([1.5, 1e7])
m  = [1.0;-0.5]
lm = 1e-6

gNPproj = project(GaussianNatParam(-nP*m, nP), minvar=lm)
gMPproj = project(GaussianMeanParam(mean=m, cov=inv(-nP)), minvar=lm)

nPthresh = -diagm([1.5, 1./lm])

@test isapprox(gNPproj.theta1, -nPthresh*m)
@test isapprox(gNPproj.theta2, -diagm([1.5, 1./lm]))
@test isapprox(gMPproj.mu1, m)
@test isapprox(cov(gMPproj), diagm([1./1.5, lm]))


text:
using PyCall
using Zygote: @adjoint
pushfirst!(PyVector(pyimport("sys")."path"), @__DIR__)
cctc = pyimport("c_ctc")

function collapse(seq)
  s = [x for x in seq if x != 62]
  if isempty(s) return s end
  s = [seq[1]]
  for ch in seq[2:end]
    if ch != s[end] && ch != 62
      push!(s, ch)
    end
  end
  return s
end

function c_ctc_(yhat, y)
  y = mapslices(argmax, y, dims=1)
  y = collapse(y)
  y = reshape(y, 1, length(y))

  l, g = cctc.c_ctc(yhat, Int32.(y .- 1))
  l = l[1]
  return l, g
end

function c_ctc(yhat, y)
  return c_ctc_(yhat, y)[1]
end

@adjoint function c_ctc(yhat, y)
  l, g = c_ctc_(yhat, y)
  return l, Δ -> (Δ .* g, Δ)
end

max_stars_repo_name: maetshju/flux-ctc-grad

max_issues_repo_name: maetshju/flux-ctc-grad

max_forks_repo_name: maetshju/flux-ctc-grad


text:
function analyzer(mat::Dict, outdir::AbstractString; dx=1, makeplots=true, isExt=false)

  R1map = mat["R10"]
  S0map = mat["S0"]
  modelmap = mat["modelmap"]
  Ct = mat["Ct"]
  Kt = mat["Kt"]
  ve = mat["ve"]
  vp = mat["vp"]
  resid = mat["resid"]
  q = quantile(S0map[:], 0.99)
  S0map[S0map .> q] .= q
  back = (S0map .- minimum(S0map)) / (maximum(S0map) - minimum(S0map))
  mask = convert(Array{Bool,2}, mat["mask"])

  # compare to known truths
  if ~isExt
    u = ones(div(10,dx),div(10,dx))
    Kt_truth = repeat([0.01*u; 0.02*u; 0.05*u; 0.1*u; 0.2*u; 0.35*u], 1, 5)
    ve_truth = repeat([0.01*u; 0.05*u; 0.1*u; 0.2*u; 0.5*u]', 6, 1)
  else
    u = ones(div(180,dx),div(10,dx))
    Kt_truth = hcat(0.01u, 0.02u, 0.05u, 0.1u, 0.2u)

    u = ones(div(10,dx),div(50,dx))
    u = vcat(0.1u, 0.2u, 0.5u)
    ve_truth = vcat(u, u, u, u, u, u)

    u = ones(div(30,dx),div(50,dx))
    vp_truth = vcat(0.001u, 0.005u, 0.01u, 0.02u, 0.05u, 0.1u)
  end

  Kt_error = clamp.(100.0*(Kt - Kt_truth) ./ (Kt_truth .+ eps()), -100.0, 100.0)
  ve_error = clamp.(100.0*(ve - ve_truth) ./ (ve_truth .+ eps()), -100.0, 100.0)
  cccKt = ccc(Kt_truth, Kt)
  cccve = ccc(ve_truth, ve)
  cccvp = 0.0 # ccc for ExtTofts is calculated in next block
  printstyled("Kt\n\tRMSE:\t$(sqrt(norm(Kt_error)^2 / length(Kt_error))) %\n", color=:green)
  printstyled("\terrmax:\t$(maximum(abs.(Kt_error)))\n", color=:green)
  printstyled("\tCCC:\t$cccKt\n", color=:green)
  printstyled("ve\n\tRMSE:\t$(sqrt(norm(ve_error)^2 / length(ve_error))) %\n", color=:green)
  printstyled("\terrmax:\t$(maximum(abs.(ve_error)))\n", color=:green)
  printstyled("\tCCC:\t$cccve\n", color=:green)

  if isExt
    vp_error = clamp.(100.0*(vp - vp_truth) ./ (vp_truth .+ eps()), -100.0, 100.0)
    cccvp = ccc(vp_truth, vp)
    printstyled("vp\n\tRMSE:\t$(sqrt(norm(vp_error)^2 / length(vp_error))) %\n", color=:green)
    printstyled("\terrmax:\t$(maximum(abs.(vp_error))) %\n", color=:green)
    printstyled("\tCCC:\t$cccvp\n", color=:green)
  end

  if !makeplots
    return (cccKt, cccve, cccvp)
  end

  if ~isExt
    ytpos = collect((0+floor(Integer, 5/dx)):div(10,dx):(div(60,dx)-1))
    xtpos = collect((0+floor(Integer, 5/dx)):div(10,dx):(div(50,dx)-1))
    ytlabels = [string(x) for x in [0.01,0.02,0.05,0.1,0.2,0.35]]
    xtlabels = [string(x) for x in [0.01,0.05,0.1,0.2,0.5]]
    # Size of figures that contain 2D maps and their x/y labels
    mapWidth = 4.5
    mapHeight = 4.5
    mapLabelX = "\$v_\\mathrm{e}\$"
    mapLabelY = "\$K^\\mathrm{trans}\$"
  else
    ytpos = collect((div(10,dx)+floor(Integer, 5/dx)):div(30,dx):(div(180,dx)-1))
    xtpos = collect((0+floor(Integer, 5/dx)):div(10,dx):(div(50,dx)-1))
    ytlabels = [string(x) for x in [0.001, 0.005, 0.01, 0.02, 0.05, 0.1]]
    xtlabels = [string(x) for x in [0.01,0.02,0.05,0.1,0.2]]
    mapWidth = 3.5
    mapHeight = 6
    mapLabelX = "\$K^\\mathrm{trans}\$"
    mapLabelY = "\$v_\\mathrm{p}\$"
  end

  println("Plotting results ...")

  # AIF
  figure(figsize=(4.5,4.5))
  clf()
  plot(mat["t"], mat["Cp"], "ko-")
  xlabel("time (min)")
  yticks(collect(0:2:10)) # This produces an error
  ylim(0,10)
  ylabel("[Gd-DTPA] (mM)")
  title("arterial input function, \$C_p\$")
  savefig("$outdir/aif.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(modelmap, interpolation="nearest", cmap="cubehelix")
  title("model used")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/modelmap.pdf",bbox_inches="tight")

  # PARAMETER MAPS
  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(Kt, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=maximum(Kt_truth))
  title("\$K^\\mathrm{trans}\$ (min\$^{-1}\$)")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/Kt.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(ve, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=maximum(ve_truth))
  title("\$v_\\mathrm{e}\$")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/ve.pdf",bbox_inches="tight")

  if isExt
    figure(figsize=(mapWidth, mapHeight))
    clf()
    imshow(vp, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=maximum(vp_truth))
    title("\$v_p\$")
    xticks(xtpos, xtlabels, fontsize=8)
    yticks(ytpos, ytlabels)
    xlabel(mapLabelX)
    ylabel(mapLabelY)
    colorbar()
    savefig("$outdir/vp.pdf",bbox_inches="tight")
  end

  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(resid, interpolation="nearest", cmap="cubehelix", vmin=0)
  title("residual")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/resid.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  m = maximum(abs.(Kt_error))
  imshow(Kt_error, interpolation="nearest", cmap="PiYG", vmin=-m, vmax=m)
  title("% error in \$K^\\mathrm{trans}\$")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/Kt_error.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  m = maximum(abs.(ve_error))
  imshow(ve_error, interpolation="nearest", cmap="PiYG", vmin=-m, vmax=m)
  title("% error in \$v_\\mathrm{e}\$")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/ve_error.pdf",bbox_inches="tight")

  if isExt
    figure(figsize=(mapWidth, mapHeight))
    clf()
    m = maximum(abs.(vp_error))
    imshow(vp_error, interpolation="nearest", cmap="PiYG", vmin=-m, vmax=m)
    title("% error in \$v_p\$")
    xticks(xtpos, xtlabels, fontsize=8)
    yticks(ytpos, ytlabels)
    xlabel(mapLabelX)
    ylabel(mapLabelY)
    colorbar()
    savefig("$outdir/vp_error.pdf",bbox_inches="tight")
  end

  return (cccKt, cccve, cccvp)
end

function analyze(n, mat::Dict, outdir::AbstractString; kwargs...)
  if n == 4
    return analyzer(mat, outdir; isExt=true, kwargs...)
  elseif n == 6
    return analyzer(mat, outdir; isExt=false, kwargs...)
  end
end

function validate(n, outdir::AbstractString; kwargs...)
  @assert n == 4 || n == 6 "n must be 4 or 6"
  cd("$(@__DIR__)/../test/q$n")

  # Create noisy data
  makeQibaNoisy(n)

  println("Running analysis of noise-free QIBA v$n data ...")
  isdir("$outdir/results") || mkdir("$outdir/results")
  results = fitdata(datafile="qiba$n.mat",outfile="$outdir/results/results.mat",save=false)
  ccc = analyze(n, results, "$outdir/results", dx=10; kwargs...)

  println("Running analysis of noisy QIBA v$n data ...")
  isdir("$outdir/results_noisy") || mkdir("$outdir/results_noisy")
  results = fitdata(datafile="qiba$(n)noisy.mat", outfile="$outdir/results_noisy/results.mat",save=false)

  cccnoisy = analyze(n, results, "$outdir/results_noisy"; kwargs...)
  println("Validation complete. If saved, results can be found in $outdir.")
  return (ccc, cccnoisy)
end

validate(n; kwargs...) = validate(n, "$(@__DIR__)/../test/q$n"; kwargs...)
function validate(; kwargs...)
  ccc6, cccnoisy6 = validate(6; kwargs...)
  ccc4, cccnoisy4 = validate(4; kwargs...)
  return (ccc6, cccnoisy6, ccc4, cccnoisy4)
end

function makeQibaNoisy(n; nRep=10, doOverwrite=true, noiseSigma=-1.0)
# Purpose: Reads in noiseless QIBA data and outputs noisy version
# Each voxel is also replicated 10 times

  # Location of noiseless data
  inFile = "$(@__DIR__)/../test/q$n/qiba$n.mat"

  # Define the output file location/name
  outFileName = basename(inFile)[1:end-4] * "noisy.mat"
  outFile = joinpath( dirname(inFile), outFileName )

  if ( isfile(outFile) && ~doOverwrite )
    # If output file already exists and we do not want to overwrite ...
    return # ... then end it right here
  end

  println("Producing noisy version of $inFile")
  # Load data
  matData = matread(inFile)

  # We can now make desired modifications to loaded data:
  matData["mask"] = Array{UInt8,2}(repeat(matData["mask"], inner=[nRep, nRep]))

  # QIBA6 and QIBA4 have different needs for T1-mapping
  if (n==4)
    # First dim of T1data is the flip angles which won't be repeated
    matData["T1data"] = repeat(matData["T1data"], inner=[1, nRep, nRep])
  else
    matData["S0"] = repeat(matData["S0"], inner=[nRep, nRep])
    matData["R10"] = repeat(matData["R10"], inner=[nRep, nRep])
  end

  # Special treatment for DCE data
  # First, use an easier-to-type variable name
  dceDat = matData["DCEdata"]
  # Repeat each element the desired number of times
  dceDat = repeat(dceDat, inner=[1, nRep, nRep]) # First dim is time, not repeated
  # Default noise = arbitraryWeight * baselineSignal / sqrt(2)
  if (noiseSigma < 0)
    noiseSigma = 0.2 * dceDat[1,1,1] / sqrt(2)
  end
  # Add complex noise
  Random.seed!(8702572558940514935) # Fixed arbitrary seed for reproducible noise
  dceDat = dceDat + noiseSigma * ( randn(size(dceDat)) + im*randn(size(dceDat)) )
  # Take the magntude of the complex signal
  dceDat = abs.(dceDat)
  # Replace noiseless data with noisy data
  matData["DCEdata"] = dceDat

  # Save modified data to disk
  matwrite(outFile, matData)
  println("Noisy data saved to $outFile")

  return
end


text:
# Van der Waals radii in angstroms
const van_der_waals_radii = Dict{Int64, Float64}(MolecularGraph.ATOM_VANDERWAALS_RADII)

function vdwradius(atom::SDFileAtom)
    return van_der_waals_radii[atomnumber(atom.symbol)]
end

function vdwradii(mol::Union{UndirectedGraph,SubgraphView})
    return Dict(zip(nodeset(mol), [vdwradius(nodeattr(mol,idx)) for idx in nodeset(mol)]))
end

"""
    vdwradius!(mol)

Add an attribute `:atomradius` to `mol`, a dictionary with the atom indices as keys and the Van der Waals
radii as values.
"""
vdwradii!(mol::UndirectedGraph) = attributes(mol)[:vdwradii] = vdwradii(mol)

function vdwradii!(mol::SubgraphView)
    pmol = mol.graph
    pvdw = get!(attributes(pmol), :vdwradii, vdwradii!(pmol))
    return filter(pr -> pr.first ∈ nodeset(mol), pvdw)
end

const const_volume_coeff = (2/(9π))^(1/3)

function vdwvolume_sigma(mol::Union{UndirectedGraph,SubgraphView})
    return Dict(zip(nodeset(mol),[const_volume_coeff * van_der_waals_radii[atomnumber(nodeattr(mol,idx).symbol)]^2 for idx in nodeset(mol)]))
end


text:
#module OneD
#using  ..Utils

#export timeOfFlight!

function trajectory!(p::AbstractPenetrator;dt=1e-6)
    ma = machNumber(p.velocity_m)
    cd0 = p.coeffAero(ma)
    ρ0 = density(0.0)
    drag = -0.125 * ρ0 * pi * p.calibre^2 * cd0 * p.velocity_m^2
    acceleration = drag / p.mass
    p.velocity_m += acceleration * dt
    p.position += p.velocity_m * dt
end

function timeOfFlight!(p::AbstractPenetrator, target::AbstractTarget;dt=1e-6)
#function flight(p::Projectile, ρ::Float64, dt::Float64)
    t = 0.0
    old_position = p.position
    old_velocity = p.velocity_m
    ρ = target.position

    while ρ > p.position
        old_position = p.position
        old_velocity = p.velocity_m
        t += dt
        #@process trajectory_projectile(sim, p, dt)
        trajectory!(p)
        #println(p.position)
    end
    flight_time = t-dt + dt * (ρ - old_position) / (p.position - old_position)
    p.position = old_position
    p.velocity_m = old_velocity + (p.velocity_m - old_velocity) * (flight_time - t + dt) / dt

    return flight_time

end

#end


text:
const MMI = MLJModelInterface

using PersistenceDiagrams: AbstractVectorizer

abstract type RipsererModel <: MMI.Unsupervised end

const PointLike{N} = AbstractVector{NTuple{N,MMI.Continuous}}
const DistanceMatrix = AbstractMatrix{MMI.Continuous}
const ImageLike{N} = AbstractArray{MMI.Continuous,N}

function _transform(model::RipsererModel, verbosity::Int, X)
    result = (; [Symbol(:dim_, i) => PersistenceDiagram[] for i in 0:(model.dim_max)]...)
    if verbosity > 0
        progbar = Progress(MMI.nrows(X); desc="Computing persistent homology...")
    end
    for x in X
        r = ripserer(_filtration(model, x); _ripserer_args(model)...)
        for i in 1:(model.dim_max + 1)
            push!(result[i], r[i])
        end
        verbosity > 0 && next!(progbar)
    end
    return MMI.table(result)
end

function MMI.fit(model::RipsererModel, verbosity::Int, X)
    diagrams = _transform(model, verbosity, X)
    vectorizers, _, report = MMI.fit(model.vectorizer, verbosity, diagrams)
    return (vectorizers, nothing, report)
end

function MMI.transform(model::RipsererModel, vectorizers, X)
    diagrams = _transform(model, 0, X)
    return MMI.transform(model.vectorizer, vectorizers, diagrams)
end

MMI.output_scitype(::Type{<:RipsererModel}) = MMI.Table(MMI.Continuous)

"""
    RipsPersistentHomology

Compute Vietoris-Rips persistent homology and convert the results to a table of continuous
values.

!!! warning Warning
    Computing Vietoris-Rips persistent homology may be CPU and memory intensive even for
    modestly-sized data sets. Consider using [`AlphaPersistentHomology`](@ref) for
    low-dimensional data.

# Hyperparameters

* `vectorizer = PersistenceImageVectorizer()`: `AbstractVectorizer` used to transform
  persistence diagrams to continuous vectors. See the [PersistenceDiagrams.jl
  Documentation](https://mtsch.github.io/PersistenceDiagrams.jl/dev/mlj/) for more
  information.

* `dim_max = 1`: compute persistent homology up to this dimension.

* `modulus = 2`: compute persistent homology with coefficients in the prime field of integers
  mod `modulus`.

* `threshold = nothing`: compute persistent homology up to diameter smaller than threshold.

* `cutoff = 0`: only keep intervals with `persistence(interval) > cutoff`.

* `sparse = false`: use a sparse Rips filtration.

* `collapse = false`: use the [`EdgeCollapsedRips`](@ref) filtration. This is usually a good
  choice for computations with a high `dim_max`.

# See also

* [`ripserer`](@ref)
* [`Rips`](@ref)
* [`EdgeCollapsedRips`](@ref)

"""
MMI.@mlj_model mutable struct RipsPersistentHomology <: RipsererModel
    vectorizer::AbstractVectorizer = PersistenceImageVectorizer()
    dim_max::Int = 1::(_ ≥ 0)
    modulus::Int = 2::(is_prime(_))
    threshold::Union{Nothing,Float64} = nothing::(isnothing(_) || _ > 0)
    cutoff::Float64 = 0::(_ ≥ 0)
    sparse::Bool = false
    collapse::Bool = false
end

function _filtration(model::RipsPersistentHomology, x)
    if isnothing(model.threshold)
        kwargs = (; threshold=model.threshold, sparse=model.sparse)
    else
        kwargs = (; sparse=model.sparse)
    end
    if model.collapse
        return EdgeCollapsedRips(x; kwargs...)
    else
        return Rips(x; kwargs...)
    end
end

function _ripserer_args(model::RipsPersistentHomology)
    return (; dim_max=model.dim_max, cutoff=model.cutoff, modulus=model.modulus)
end

function MMI.input_scitype(::Type{<:RipsPersistentHomology})
    return Union{
        MMI.Table(PointLike{N}) where {N},
        MMI.Table(DistanceMatrix),
        AbstractVector{PointLike{N}} where {N},
        AbstractVector{DistanceMatrix},
    }
end

"""
    AlphaPersistentHomology

Compute alpha complex persistent homology and convert the results to a table of continuous
values.

!!! warning Warning
    Using high-dimensional data with this model may be computationaly expensive. Consider
    using [`RipsPersistentHomology`](@ref).

!!! warning
    This model uses [MiniQhull.jl](https://github.com/gridap/MiniQhull.jl). Please see
    the installation instructions if fitting causes errors. MiniQhull currently has
    problems running on Windows. See [this
    issue](https://github.com/gridap/MiniQhull.jl/issues/5) for more info.

# Hyperparameters

* `vectorizer = PersistenceImageVectorizer()`: `AbstractVectorizer` used to transform
  persistence diagrams to continuous vectors. See the [PersistenceDiagrams.jl
  Documentation](https://mtsch.github.io/PersistenceDiagrams.jl/dev/mlj/) for more
  information.

* `dim_max = 1`: compute persistent homology up to this dimension.

* `modulus = 2`: compute persistent homology with coefficients in the prime field of integers
  mod `modulus`.

* `threshold = nothing`: compute persistent homology up to diameter smaller than threshold.

* `cutoff = 0`: only keep intervals with `persistence(interval) > cutoff`.

# See also

* [`ripserer`](@ref)
* [`Alpha`](@ref)

"""
MMI.@mlj_model mutable struct AlphaPersistentHomology <: RipsererModel
    vectorizer::AbstractVectorizer = PersistenceImageVectorizer()
    dim_max::Int = 1::(_ ≥ 0)
    modulus::Int = 2::(is_prime(_))
    threshold::Union{Nothing,Float64} = nothing::(isnothing(_) || _ > 0)
    cutoff::Float64 = 0::(_ ≥ 0)
end

function _filtration(model::AlphaPersistentHomology, data)
    if isnothing(model.threshold)
        return Alpha(data)
    else
        return Alpha(data; threshold=model.threshold)
    end
end

function _ripserer_args(model::AlphaPersistentHomology)
    return (; dim_max=model.dim_max, cutoff=model.cutoff, modulus=model.modulus)
end

function MMI.input_scitype(::Type{<:AlphaPersistentHomology})
    return Union{MMI.Table(PointLike{N}),AbstractVector{PointLike{N}}} where {N}
end

"""
    CubicalPersistentHomology

Compute cubical persistent homology and convert the results to a table of continuous
values.

# Hyperparameters

* `vectorizer = PersistenceImageVectorizer()`: `AbstractVectorizer` used to transform
  persistence diagrams to continuous vectors. See the [PersistenceDiagrams.jl
  Documentation](https://mtsch.github.io/PersistenceDiagrams.jl/dev/mlj/) for more
  information.

* `dim_max = 1`: compute persistent homology up to this dimension.

* `threshold = nothing`: compute persistent homology up to diameter smaller than threshold.

* `cutoff = 0`: only keep intervals with `persistence(interval) > cutoff`.

* `negate = false`: negate the image before computation.

# See also

* [`ripserer`](@ref)
* [`Cubical`](@ref)

"""
MMI.@mlj_model mutable struct CubicalPersistentHomology <: RipsererModel
    vectorizer::AbstractVectorizer = PersistenceImageVectorizer()
    dim_max::Int = 1::(_ ≥ 0)
    threshold::Union{Nothing,Float64} = nothing::(isnothing(_) || _ > 0)
    cutoff::Float64 = 0::(_ ≥ 0)
    negate::Bool = false
end

function _filtration(model::CubicalPersistentHomology, data)
    data = model.negate ? -data : data
    if isnothing(model.threshold)
        return Cubical(data)
    else
        return Cubical(data; threshold=model.threshold)
    end
end

function _ripserer_args(model::CubicalPersistentHomology)
    return (; dim_max=model.dim_max, cutoff=model.cutoff)
end

function MMI.input_scitype(::Type{<:CubicalPersistentHomology})
    return Union{
        MMI.Table(ImageLike{N}) where {N},
        MMI.Table(MMI.Image),
        AbstractVector{ImageLike{N}} where {N},
        MMI.AbstractVector{MMI.Image},
    }
end

MMI.metadata_pkg.(
    (RipsPersistentHomology, AlphaPersistentHomology, CubicalPersistentHomology),
    name="Ripserer",
    uuid="aa79e827-bd0b-42a8-9f10-2b302677a641",
    url="https://github.com/mtsch/Ripserer.jl",
    license="MIT",
    julia=true,
    is_wrapper=false,
)