text:
# This file is part of the TaylorIntegration.jl package; MIT licensed

using TaylorIntegration, Elliptic
using LinearAlgebra: norm
using Test
using Logging
import Logging: Warn

@testset "Testing `jettransport.jl`" begin

    local _order = 28
    local _abstol = 1.0E-20

    max_iters_reached() = "Maximum number of integration steps reached; exiting.\n"

    f(x, p, t) = x^2
    g(x, p, t) = 0.3x

    @testset "Test TaylorN jet transport (t0, tmax): 1-dim case" begin
        p = set_variables("ξ", numvars=1, order=5)
        x0 = 3.0 #"nominal" initial condition
        x0TN = x0 + p[1] #jet transport initial condition
        t0=0.0
        tmax=0.3
        tvTN, xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            f, x0TN, t0, tmax, _order, _abstol, maxsteps=1))
        @test size(tvTN) == (2,)
        @test size(xvTN) == (2,)
        tvTN, xvTN = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0TN, t0, tmax, _order, _abstol))
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0, t0, tmax, _order, _abstol))
        exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution
        δsol = exactsol(tvTN[end], x0TN, t0)-xvTN[end]
        δcoeffs = map(y->y[1], map(x->x.coeffs, δsol.coeffs))
        @test isapprox(δcoeffs, zeros(6), atol=1e-10, rtol=0)
        @test isapprox(x0, evaluate(xvTN[1]))
        @test isapprox(xv[1], evaluate(xvTN[1]))
        @test isapprox(xv[end], evaluate(xvTN[end]))
        xvTN_0 = evaluate.(xvTN)
        @test norm(exactsol.(tvTN, x0, t0)-xvTN_0, Inf) < 1E-12
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            x0_disp = x0+disp
            dv = map(x->[disp], tvTN) #a vector of identical displacements
            tv_disp, xv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0_disp, t0, tmax, _order, _abstol))
            xvTN_disp = evaluate.(xvTN, dv)
            @test norm(exactsol.(tvTN, x0_disp, t0)-xvTN_disp, Inf) < 1E-12 #analytical vs jet transport
            @test norm(x0_disp-evaluate(xvTN[1], [disp]), Inf) < 1E-12
            @test norm(xv_disp[1]-evaluate(xvTN[1], [disp]), Inf) < 1E-12
            @test norm(xv_disp[end]-evaluate(xvTN[end], [disp]), Inf) < 1E-12
        end

        y0 = 1.0 #"nominal" initial condition
        u0 = 0.0 #initial time
        y0TN = y0 + p[1] #jet transport initial condition
        uvTN, yvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            g, y0TN, u0, 10/0.3, _order, _abstol, maxsteps=1))
        @test size(uvTN) == (2,)
        @test size(yvTN) == (2,)
        uvTN, yvTN = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0TN, u0, 10/0.3, _order, _abstol))
        uv, yv = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0, u0, 10/0.3, _order, _abstol))
        exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0))
        δsol_g = exactsol_g(uvTN[end], y0TN, u0)-yvTN[end]
        δcoeffs_g = map(y->y[1], map(x->x.coeffs, δsol_g.coeffs))
        @test isapprox(δcoeffs_g, zeros(6), atol=1e-10, rtol=0)
        @test isapprox(y0, evaluate(yvTN[1]), atol=1e-10, rtol=0)
        @test isapprox(yv[1], evaluate(yvTN[1]), atol=1e-10, rtol=0)
        @test isapprox(yv[end], evaluate(yvTN[end]), atol=1e-10, rtol=0)
        yvTN_0 = evaluate.(yvTN)
        @test norm(exactsol_g.(uvTN, y0, u0)-yvTN_0, Inf) < 1E-9
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            y0_disp = y0+disp
            dv = map(x->[disp], uvTN) #a vector of identical displacements
            uv_disp, yv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                g, y0_disp, u0, 10/0.3, _order, _abstol))
            yvTN_disp = evaluate.(yvTN, dv)
            @test norm(exactsol_g.(uvTN, y0_disp, u0)-yvTN_disp, Inf) < 1E-9 #analytical vs jet transport
            @test norm(y0_disp-yvTN[1]([disp]), Inf) < 1E-9
            @test norm(yv_disp[1]-yvTN[1]([disp]), Inf) < 1E-9
            @test norm(yv_disp[end]-yvTN[end]([disp]), Inf) < 1E-9
        end
    end

    @testset "Test Taylor1 jet transport (t0, tmax): 1-dim case" begin
        p = Taylor1([0.0,1.0], 5)
        x0 = 3.0 #"nominal" initial condition
        x0T1 = x0 + p #jet transport initial condition
        t0 = 0.0
        tmax = 0.3
        tvT1, xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            f, x0T1, t0, tmax, _order, _abstol, maxsteps=1))
        @test size(tvT1) == (2,)
        @test size(xvT1) == (2,)
        tvT1, xvT1 = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0T1, t0, tmax, _order, _abstol))
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0, t0, tmax, _order, _abstol))
        exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution
        δsol = exactsol(tvT1[end], x0T1, t0)-xvT1[end] #analytical vs jet transport diff at end of integration
        δcoeffs = δsol.coeffs
        @test isapprox(δcoeffs, zeros(6), atol=1e-10, rtol=0)
        @test isapprox(x0, evaluate(xvT1[1]))
        @test isapprox(xv[1], evaluate(xvT1[1]))
        @test isapprox(xv[end], evaluate(xvT1[end]))
        xvT1_0 = evaluate.(xvT1)
        @test norm(exactsol.(tvT1, x0, t0)-xvT1_0, Inf) < 1E-12
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            x0_disp = x0+disp
            tv_disp, xv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0_disp, t0, tmax, _order, _abstol))
            xvT1_disp = evaluate.(xvT1, disp)
            @test norm(exactsol.(tvT1, x0_disp, t0)-xvT1_disp, Inf) < 1E-12 #analytical vs jet transport
            @test norm(x0_disp-evaluate(xvT1[1], disp), Inf) < 1E-12
            @test norm(xv_disp[1]-evaluate(xvT1[1], disp), Inf) < 1E-12
            @test norm(xv_disp[end]-evaluate(xvT1[end], disp), Inf) < 1E-12
        end

        y0 = 1.0 #"nominal" initial condition
        u0 = 0.0 #initial time
        y0T1 = y0 + p #jet transport initial condition
        uvT1, yvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            g, y0T1, t0, 10/0.3, _order, _abstol, maxsteps=1)) #warmup lap
        # test maxsteps break
        @test size(uvT1) == (2,)
        @test size(yvT1) == (2,)
        uvT1, yvT1 = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0T1, t0, 10/0.3, _order, _abstol)) #Taylor1 jet transport integration
        uv, yv = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0, t0, 10/0.3, _order, _abstol)) #reference integration
        exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0))  #the analytical solution
        δsol_g = exactsol_g(uvT1[end], y0T1, t0)-yvT1[end] #analytical vs jet transport diff at end of integration
        δcoeffs_g = δsol_g.coeffs
        @test isapprox(δcoeffs_g, zeros(6), atol=1e-10, rtol=0)
        @test isapprox(y0, evaluate(yvT1[1]), atol=1e-10, rtol=0)
        @test isapprox(yv[1], evaluate(yvT1[1]), atol=1e-10, rtol=0)
        @test isapprox(yv[end], evaluate(yvT1[end]), atol=1e-10, rtol=0)
        yvT1_0 = evaluate.(yvT1)
        @test norm(exactsol_g.(uvT1, y0, u0)-yvT1_0, Inf) < 1E-9
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            y0_disp = y0+disp
            uv_disp, yv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                g, y0_disp, u0, 10/0.3, _order, _abstol))
            yvT1_disp = evaluate.(yvT1, disp)
            @test norm(exactsol_g.(uvT1, y0_disp, u0)-yvT1_disp, Inf) < 1E-9 #analytical vs jet transport
            @test norm(y0_disp-evaluate(yvT1[1], disp), Inf) < 1E-9
            @test norm(yv_disp[1]-evaluate(yvT1[1], disp), Inf) < 1E-9
            @test norm(yv_disp[end]-evaluate(yvT1[end], disp), Inf) < 1E-9
        end
    end

    @testset "Test Taylor1 jet transport (trange): 1-dim case" begin
        p = Taylor1([0.0,1.0], 5)
        x0 = 3.0 #"nominal" initial condition
        x0T1 = x0 + p #jet transport initial condition
        tv = 0.0:0.05:0.33
        xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            f, x0T1, tv, _order, _abstol, maxsteps=1))
        @test size(xvT1) == (7,)
        ta = vec(tv)
        xvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            f, x0T1, ta, _order, _abstol, maxsteps=1))
        @test size(xvT1) == (7,)
        xvT1 = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0T1, tv, _order, _abstol))
        xvT11 = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0T1, ta, _order, _abstol))
        @test xvT11 == xvT1
        xv = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0, tv, _order, _abstol))
        exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution
        δsol = exactsol(tv[end], x0T1, tv[1])-xvT1[end]
        δcoeffs = δsol.coeffs
        @test length(tv) == length(xvT1)
        @test isapprox(δcoeffs, zeros(6), atol=1e-10, rtol=0)
        xv_analytical = exactsol.(tv, x0, tv[1])
        xvT1_0 = xvT1.()
        @test isapprox(xv_analytical, xvT1_0, atol=1e-10, rtol=0)
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            xv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0+disp, tv, _order, _abstol))
            xv_disp2 = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0+disp, ta, _order, _abstol))
            @test xv_disp == xv_disp2
            xvT1_disp = xvT1.(disp)
            @test norm(exactsol.(tv, x0+disp, tv[1])-xvT1_disp, Inf) < 1E-12 #analytical vs jet transport
            @test norm(xv_disp-xvT1_disp, Inf) < 1E-12 # integration vs jet transport
        end

        y0 = 1.0 #"nominal" initial condition
        y0T1 = y0 + p #jet transport initial condition
        uv = 0.0:1/0.3:10/0.3
        yvT1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            g, y0T1, uv, _order, _abstol, maxsteps=1))
        @test size(yvT1) == (11,)
        yvT1 = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0T1, uv, _order, _abstol))
        yv = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0, uv, _order, _abstol))
        exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0))
        δsol_g = exactsol_g(uv[end], y0T1, uv[1])-yvT1[end]
        δcoeffs_g = δsol_g.coeffs
        @test length(uv) == length(yvT1)
        @test isapprox(δcoeffs_g, zeros(6), atol=1e-10, rtol=0)
        yv_analytical = exactsol_g.(uv, y0, uv[1])
        yvT1_0 = evaluate.(yvT1)
        @test isapprox(yv_analytical, yvT1_0, atol=1e-10, rtol=0)
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            yv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                g, y0+disp, uv, _order, _abstol))
            yvT1_disp = evaluate.(yvT1, disp)
            @test norm(exactsol_g.(uv, y0+disp, uv[1])-yvT1_disp, Inf) < 1E-9 #analytical vs jet transport
            @test norm(yv_disp-yvT1_disp, Inf) < 1E-9 # integration vs jet transport
        end
    end

    @testset "Test TaylorN jet transport (trange): 1-dim case" begin
        p = set_variables("ξ", numvars=1, order=5)
        x0 = 3.0 #"nominal" initial condition
        x0TN = x0 + p[1] #jet transport initial condition
        tv = 0.0:0.05:0.33
        xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            f, x0TN, tv, _order, _abstol, maxsteps=1))
        @test size(xvTN) == (7,)
        ta = vec(tv)
        xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            f, x0TN, ta, _order, _abstol, maxsteps=1))
        @test size(xvTN) == (7,)
        xvTN = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0TN, tv, _order, _abstol))
        xvTN_ = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0TN, ta, _order, _abstol))
        @test xvTN == xvTN_
        xv = (@test_logs min_level=Logging.Warn taylorinteg(
            f, x0, tv, _order, _abstol))
        exactsol(t, x0, t0) = x0/(1.0-x0*(t-t0)) #the analytical solution
        δsol = exactsol(tv[end], x0TN, tv[1])-xvTN[end]
        δcoeffs = map(y->y[1], map(x->x.coeffs, δsol.coeffs))
        @test length(tv) == length(xvTN)
        @test isapprox(δcoeffs, zeros(6), atol=1e-10, rtol=0)
        xv_analytical = exactsol.(tv, x0, tv[1])
        xvTN_0 = evaluate.(xvTN)
        @test isapprox(xv_analytical, xvTN_0, atol=1e-10, rtol=0)
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            dv = map(x->[disp], tv) #a vector of identical displacements
            xv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0+disp, tv, _order, _abstol))
            xv_disp_ = (@test_logs min_level=Logging.Warn taylorinteg(
                f, x0+disp, ta, _order, _abstol))
            @test xv_disp == xv_disp_
            xvTN_disp = evaluate.(xvTN, dv)
            @test norm(exactsol.(tv, x0+disp, tv[1])-xvTN_disp, Inf) < 1E-12 #analytical vs jet transport
            @test norm(xv_disp-xvTN_disp, Inf) < 1E-12 # integration vs jet transport
        end

        y0 = 1.0 #"nominal" initial condition
        y0TN = y0 + p[1] #jet transport initial condition
        uv = 0.0:1/0.3:10/0.3
        yvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            g, y0TN, uv, _order, _abstol, maxsteps=1))
        @test size(yvTN) == (11,)
        yvTN = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0TN, uv, _order, _abstol))
        yv = (@test_logs min_level=Logging.Warn taylorinteg(
            g, y0, uv, _order, _abstol))
        exactsol_g(u, y0, u0) = y0*exp(0.3(u-u0))
        δsol_g = exactsol_g(uv[end], y0TN, uv[1])-yvTN[end]
        δcoeffs_g = map(y->y[1], map(x->x.coeffs, δsol_g.coeffs))
        @test length(uv) == length(yvTN)
        @test isapprox(δcoeffs_g, zeros(6), atol=1e-10, rtol=0)
        yv_analytical = exactsol_g.(uv, y0, uv[1])
        yvTN_0 = evaluate.(yvTN)
        @test isapprox(yv_analytical, yvTN_0, atol=1e-10, rtol=0)
        for i in 1:5
            disp = 0.001*rand() #a small, random displacement
            dv = map(x->[disp], uv) #a vector of identical displacements
            yv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                g, y0+disp, uv, _order, _abstol))
            yvTN_disp = evaluate.(yvTN, dv)
            @test norm(exactsol_g.(uv, y0+disp, uv[1])-yvTN_disp, Inf) < 1E-9 #analytical vs jet transport
            @test norm(yv_disp-yvTN_disp, Inf) < 1E-9 # integration vs jet transport
        end
    end

    let
        function harmosc!(dx, x, p, t) #the harmonic oscillator ODE
            dx[1] = x[2]
            dx[2] = -x[1]*x[3]^2
            dx[3] = zero(x[1])
            nothing
        end

        @testset "Test Taylor1 jet transport (t0,tmax): harmonic oscillator" begin
            t = Taylor1([0.0, 1.0], 10)
            ω0 = 1.0
            x0 = [0.0,ω0,ω0]
            x0T1 = x0+[0t,t,t]
            tv1, xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
                harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=1))
            @test size(tv1) == (2,)
            @test size(xv1) == (2, 3)
            tv1, xv1 = (@test_logs min_level=Logging.Warn taylorinteg(
                harmosc!, x0T1, 0.0, 100pi, _order, _abstol, maxsteps=2000))
            y0 = evaluate.(xv1)
            x1(t,δω) = sin((ω0+δω)*t)
            x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t)
            @test norm(y0[:,1]-x1.(tv1,0.0), Inf) < 1E-11
            @test norm(y0[:,2]-x2.(tv1,0.0), Inf) < 1E-11
            for i in 1:5
                δω=0.001*rand()
                x0_disp = x0+[0.0,δω,δω]
                tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
                    harmosc!, x0_disp, 0.0, 100pi, _order, _abstol, maxsteps=2000))
                y1 = evaluate.(xv1, δω)
                @test norm(y1[:,1]-x1.(tv1,δω),Inf) < 1E-11
                @test norm(y1[:,2]-x2.(tv1,δω),Inf) < 1E-11
                @test norm(x0_disp-evaluate.(xv1[1,:], δω), Inf) < 1E-11
                @test norm(xv[1,:]-evaluate.(xv1[1,:], δω), Inf) < 1E-11
                @test norm(xv[end,:]-evaluate.(xv1[end,:], δω), Inf) < 1E-11
            end
        end

        @testset "Test Taylor1 jet transport (trange): harmonic oscillator" begin
            t = Taylor1([0.0, 1.0], 10)
            ω0 = 1.0
            x0 = [0.0,ω0,ω0]
            x0T1 = x0+[0t,t,t]
            tv = 0.0:0.25*(2pi):100pi
            xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
                harmosc!, x0T1, tv, _order, _abstol, maxsteps=1))
            @test length(tv) == 201
            @test size(xv1) == (201, 3)
            ta = vec(tv)
            xv1 = (@test_logs (Warn, max_iters_reached()) taylorinteg(
                harmosc!, x0T1, ta, _order, _abstol, maxsteps=1))
            @test length(ta) == 201
            @test size(xv1) == (201, 3)
            xv1 = (@test_logs min_level=Logging.Warn taylorinteg(
                harmosc!, x0T1, tv, _order, _abstol, maxsteps=2000))
            xv1_ = (@test_logs min_level=Logging.Warn taylorinteg(
                harmosc!, x0T1, ta, _order, _abstol, maxsteps=2000))
            @test xv1 == xv1_
            y0 = evaluate.(xv1)
            x1(t,δω) = sin((ω0+δω)*t)
            x2(t,δω) = (ω0+δω)*cos((ω0+δω)*t)
            @test norm(y0[:,1]-x1.(tv,0.0), Inf) < 1E-11
            @test norm(y0[:,2]-x2.(tv,0.0), Inf) < 1E-11
            for i in 1:5
                δω=0.001*rand()
                x0_disp = x0+[0.0,δω,δω]
                xv = (@test_logs min_level=Logging.Warn taylorinteg(
                    harmosc!, x0_disp, tv, _order, _abstol, maxsteps=2000))
                xv_ = (@test_logs min_level=Logging.Warn taylorinteg(
                    harmosc!, x0_disp, ta, _order, _abstol, maxsteps=2000))
                @test xv == xv_
                y1 = evaluate.(xv1, δω)
                @test norm(y1[:,1]-x1.(tv,δω), Inf) < 1E-11
                @test norm(y1[:,2]-x2.(tv,δω), Inf) < 1E-11
                @test norm(y1-xv, Inf) < 1E-11
            end
        end
    end

    let
        function harmosc!(dx, x, p, t) #the harmonic oscillator ODE
            dx[1] = x[2]
            dx[2] = -x[1]
            nothing
        end

        @testset "Test TaylorN jet transport (t0,tmax): harmonic oscillator" begin
            p = set_variables("ξ", numvars=2, order=5)
            x0 = [-1.0,0.45]
            x0TN = x0 + p
            tvTN, xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
                harmosc!, x0TN, 0.0, 10pi, _order, _abstol, maxsteps=1))
            @test length(tvTN) == 2
            @test size(xvTN) == (length(tvTN), length(x0TN))
            tvTN, xvTN = (@test_logs min_level=Logging.Warn taylorinteg(
                harmosc!, x0TN, 0.0, 10pi, _order, _abstol))
            tv , xv = (@test_logs min_level=Logging.Warn taylorinteg(
                harmosc!, x0  , 0.0, 10pi, _order, _abstol))
            x_analyticsol(t,x0,p0) = p0*sin(t)+x0*cos(t)
            p_analyticsol(t,x0,p0) = p0*cos(t)-x0*sin(t)
            x_δsol = x_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,1]
            x_δcoeffs = map(y->y[1], map(x->x.coeffs, x_δsol.coeffs))
            p_δsol = p_analyticsol(tvTN[end], x0TN[1], x0TN[2])-xvTN[end,2]
            p_δcoeffs = map(y->y[1], map(x->x.coeffs, p_δsol.coeffs))

            @test (length(tvTN), length(x0)) == size(xvTN)
            @test isapprox(x_δcoeffs, zeros(6), atol=1e-10, rtol=0)
            @test isapprox(x0, map( x->evaluate(x), xvTN[1,:]))
            @test isapprox(xv[1,:], map( x->evaluate(x), xvTN[1,:])) # nominal solution must coincide with jet evaluated at ξ=(0,0) at initial time
            @test isapprox(xv[end,:], map( x->evaluate(x), xvTN[end,:])) #nominal solution must coincide with jet evaluated at ξ=(0,0) at final time

            xvTN_0 = map( x->evaluate(x), xvTN ) # the jet evaluated at the nominal solution
            @test isapprox(xv[end,:], xvTN_0[end,:]) # nominal solution must coincide with jet evaluated at ξ=(0,0) at final time
            @test isapprox(xvTN_0[1,:], xvTN_0[end,:]) # end point must coincide with a full period
        end
    end

    function pendulum!(dx, x, p, t) #the simple pendulum ODE
        dx[1] = x[2]
        dx[2] = -sin(x[1])
        nothing
    end

    @testset "Test TaylorN jet transport (trange): simple pendulum" begin
        varorder = 2 #the order of the variational expansion
        p = set_variables("ξ", numvars=2, order=varorder) #TaylorN steup
        q0 = [1.3, 0.0] #the initial conditions
        q0TN = q0 + p #parametrization of a small neighbourhood around the initial conditions
        # T is the librational period == 4Elliptic.K(sin(q0[1]/2)^2)
        T = 4Elliptic.K(sin(q0[1]/2)^2) # equals 7.019250311844546
        t0 = 0.0 #the initial time
        tmax = T #the final time
        integstep = 0.25*T #the time interval between successive evaluations of the solution vector

        #the time range
        tr = t0:integstep:tmax;
        #xv is the solution vector representing the propagation of the initial condition q0 propagated until time T
        xv = (@test_logs min_level=Logging.Warn taylorinteg(
            pendulum!, q0, tr, _order, _abstol, maxsteps=100))
        #xvTN is the solution vector representing the propagation of the initial condition q0 plus variations (ξ₁,ξ₂) propagated until time T
        #note that q0 is a Vector{Float64}, but q0TN is a Vector{TaylorN{Float64}}
        #but apart from that difference, we're calling `taylorinteg` essentially with the same parameters!
        #thus, jet transport is reduced to a beautiful application of Julia's multiple dispatch!
        xvTN = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            pendulum!, q0TN, tr, _order, _abstol, maxsteps=1))
        @test size(xvTN) == (5,2)
        xvTN = (@test_logs min_level=Logging.Warn taylorinteg(
            pendulum!, q0TN, tr, _order, _abstol, maxsteps=100))

        xvTN_0 = map( x->evaluate(x, [0.0, 0.0]), xvTN ) # the jet evaluated at the nominal solution

        @test isapprox(xvTN_0[1,:], xvTN_0[end,:]) #end point must coincide with a full period
        @test isapprox(xv, xvTN_0) #nominal solution must coincide with jet evaluated at ξ=(0,0)

        #testing another jet transport method:
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            pendulum!, q0, t0, tmax, _order, _abstol, maxsteps=100))
        @test isapprox(xv[1,:], xvTN_0[1,:]) #nominal solution must coincide with jet evaluated at ξ=(0,0) at initial time
        @test isapprox(xv[end,:], xvTN_0[end,:]) #nominal solution must coincide with jet evaluated at ξ=(0,0) at final time

        # a small displacement
        disp = 0.0001 #works even with 0.001, but we're giving it some margin

        #compare the jet solution evaluated at various variation vectors ξ, wrt to full solution, at each time evaluation point
        xv_disp = (@test_logs (Warn, max_iters_reached()) taylorinteg(
            pendulum!, q0+[disp,0.0], tr, _order, _abstol, maxsteps=1))
        @test size(xv_disp) == (5,2)
        for i in 1:10
            # generate a random angle
            ϕ = 2pi*rand()
            # generate a random point on a circumference of radius disp
            ξ = disp*[cos(ϕ), sin(ϕ)]
            #propagate in time full solution with initial condition q0+ξ
            xv_disp = (@test_logs min_level=Logging.Warn taylorinteg(
                pendulum!, q0+ξ, tr, _order, _abstol, maxsteps=100))
            #evaluate jet at q0+ξ
            xvTN_disp = map( x->evaluate(x, ξ), xvTN )
            #the propagated exact solution at q0+ξ should be approximately equal to the propagated jet solution evaluated at the same point:
            @test isapprox( xvTN_disp, xv_disp )
        end
    end
end


text:
"trying to demonstrate an unwanted behaviour in Julia"

function makeRand ()
    thing::Int = randomDevice()
    throw (BoundsError())
    return (thing)
    end

function checker (number::Int)
    done=0
    try
        while ((rem(makeRand(),100)!=0)&&done==0)
          done=1
          print("done")
          return (true)
        end
            print ("not")
            return (false)
        catch e
        println (e)
        return(true)
    end
end

ext: jl


text:
ch = NEAT.Chromosome(g, 1, 2, NEAT.Recurrent(),:ConnectionGene)
@test typeof(ch) == NEAT.Chromosome
@test typeof(ch.id) == Int64
@test typeof(ch.inputCnt) == Int64
@test typeof(ch.outputCnt) == Int64
@test typeof(ch.node_gene_type) == NEAT.Recurrent
@test typeof(ch.conn_gene_type) == Symbol
# @test typeof(ch.connection_genes) == Symbol
# @test typeof(ch.node_genes) == Symbol
@test typeof(ch.fitness) == Float64
@test typeof(ch.species_id) == Int64
@test typeof(ch.parent1_id) == Int64
@test typeof(ch.parent2_id) == Int64

g = NEAT.Global(config)
NEAT.incChromeId!(g);NEAT.incChromeId!(g)
@test g.chromosomeCnt == 2

ch = NEAT.create_minimally_connected(g)
@test length(ch.node_genes) == g.cf.input_nodes + g.cf.output_nodes
@test length(ch.connection_genes) == g.cf.output_nodes

ch = NEAT.create_fully_connected(g)
@test length(ch.node_genes) == g.cf.input_nodes + g.cf.output_nodes
@test length(ch.connection_genes) == g.cf.input_nodes * g.cf.output_nodes
@test NEAT.maxInnov(ch.connection_genes).innovNumber == 3

nodeCnt = length(ch.node_genes)
ng, cg = NEAT.mutate_add_node!(ch,g,NEAT.Recurrent())
@test cg.enable == false
@test length(ch.connection_genes) == g.cf.input_nodes * g.cf.output_nodes + 2
@test length(ch.node_genes) == nodeCnt + 1

connectCnt = length(ch.connection_genes)
NEAT.mutate_add_connection!(ch, g, NEAT.Recurrent())
@test length(ch.connection_genes) ==connectCnt + 1

@test NEAT.size(ch) == (1,5)

ch1 = NEAT.create_fully_connected(g)
ng, cg = NEAT.mutate_add_node!(ch1,g ,NEAT.Recurrent())
ng, cg = NEAT.mutate_add_node!(ch1,g ,NEAT.Recurrent())
ng, cg = NEAT.mutate_add_node!(ch1,g ,NEAT.Recurrent())
@test length(ch1.node_genes) == 7
@test length(ch1.connection_genes) == 9

ch2 = NEAT.create_fully_connected(g)
ng, cg = NEAT.mutate_add_node!(ch2,g ,NEAT.Recurrent())
ng, cg = NEAT.mutate_add_node!(ch2,g ,NEAT.Recurrent())
@test length(ch2.node_genes) == 6
@test length(ch2.connection_genes) == 7

@test NEAT.distance(NEAT.
    create_fully_connected(g),NEAT.create_fully_connected(g),g.cf) > 0

child = NEAT.Chromosome(g, ch1.id, ch2.id, NEAT.Recurrent(),:ConnectionGene)
NEAT.inherit_genes!(g, child, ch1, ch2)
@test length(child.node_genes) == 7 # test that it is larger of two
@test length(child.connection_genes) == 9 # test that it is larger of two

child = NEAT.crossover(g, ch1, ch2)
@test length(child.node_genes) == 7 # test that it is larger of two
@test length(child.connection_genes) == 9 # test that it is larger of two

#----------  FeedForward Test  ----------

ch = NEAT.Chromosome(g, 1, 2, NEAT.FeedForward(),:ConnectionGene)
@test typeof(ch) == NEAT.Chromosome
@test typeof(ch.node_gene_type) == NEAT.FeedForward
@test typeof(ch.node_order) == Vector{Int64}
@test length(ch.node_order) == 0

g = NEAT.Global(config)
g.cf.feedforward = true
ch2 = NEAT.create_fully_connected(g)
@test typeof(ch2.node_gene_type) == NEAT.FeedForward

ng, cg = NEAT.mutate_add_node!(ch2, g, NEAT.FeedForward())
@test length(ch2.node_order) == 1
@test ch2.node_order[1] == ng.id

ng, cg = NEAT.mutate_add_node!(ch2, g, NEAT.FeedForward())
@test length(ch2.node_order) == 2

ng, cg = NEAT.mutate_add_node!(ch2, g, NEAT.FeedForward())
@test length(ch2.node_order) == 3
# IMPORTANT TO DO: Write function to test for cyclical graph recurrence

NEAT.mutate_add_connection!(ch2, g, NEAT.FeedForward())
@test length(ch2.connection_genes) == 10
NEAT.mutate_add_connection!(ch2, g, NEAT.FeedForward())
@test length(ch2.connection_genes) == 11

g.cf.feedforward = true
ch2 = NEAT.create_unconnected(g)
NEAT.add_hidden_nodes!(g,ch2, 4,NEAT.FeedForward())
@test length(ch2.node_order) == 4


text:
function xavier(in_dim, out_dim)
    return randn(in_dim, out_dim) ./ sqrt(in_dim)
end

function he(in_dim, out_dim)
    return randn(in_dim, out_dim) ./ sqrt(2 / in_dim)
end


text:
# Classical Multidimensional Scaling

## convert Gram matrix to Distance matrix

"""
    gram2dmat!(D, G)

Convert a Gram matrix `G` to a distance matrix, and write the results to `D`.
"""
function gram2dmat!(D::AbstractMatrix{DT}, G::AbstractMatrix) where {DT<:Real}
    # argument checking
    m = size(G, 1)
    n = size(G, 2)
    m == n || error("D should be a square matrix.")
    size(D) == (m, n) ||
        throw(DimensionMismatch("Sizes of D and G do not match."))

    # implementation
    for j = 1:n
        for i = 1:j-1
            @inbounds D[i,j] = D[j,i]
        end
        D[j,j] = zero(DT)
        for i = j+1:n
            @inbounds D[i,j] = sqrt(G[i,i] + G[j,j] - 2 * G[i,j])
        end
    end
    return D
end

"""
    gram2dmat(G)

Convert a Gram matrix `G` to a distance matrix.
"""
gram2dmat(G::AbstractMatrix{T}) where {T<:Real} = gram2dmat!(similar(G, T), G)

## convert Distance matrix to Gram matrix

"""
    dmat2gram!(G, D)

Convert a distance matrix `D` to a Gram matrix, and write the results to `G`.
"""
function dmat2gram!(G::AbstractMatrix{GT}, D::AbstractMatrix) where GT
    # argument checking
    n = LinearAlgebra.checksquare(D)
    size(G) == (n, n) ||
        throw(DimensionMismatch("Sizes of G and D do not match."))

    # implementation
    u = zeros(GT, n)
    s = 0.0
    for j = 1:n
        s += (u[j] = sum(abs2, view(D,:,j)) / n)
    end
    s /= n

    for j = 1:n
        for i = 1:j-1
            @inbounds G[i,j] = G[j,i]
        end
        for i = j:n
            @inbounds G[i,j] = (u[i] + u[j] - abs2(D[i,j]) - s) / 2
        end
    end
    return G
end

momenttype(T) = typeof((zero(T) * zero(T) + zero(T) * zero(T))/ 2)

"""
    dmat2gram(D)

Convert a distance matrix `D` to a Gram matrix.
"""
dmat2gram(D::AbstractMatrix{T}) where {T<:Real} = dmat2gram!(similar(D, momenttype(T)), D)

## Classical MDS

"""
*Classical Multidimensional Scaling* (MDS), also known as Principal Coordinates Analysis (PCoA),
is a specific technique in this family that accomplishes the embedding in two steps:

1. Convert the distance matrix to a Gram matrix. This conversion is based on
the following relations between a distance matrix ``D`` and a Gram matrix ``G``:

```math
\\mathrm{sqr}(\\mathbf{D}) = \\mathbf{g} \\mathbf{1}^T + \\mathbf{1} \\mathbf{g}^T - 2 \\mathbf{G}
```

Here, ``\\mathrm{sqr}(\\mathbf{D})`` indicates the element-wise square of ``\\mathbf{D}``,
and ``\\mathbf{g}`` is the diagonal elements of ``\\mathbf{G}``. This relation is
itself based on the following decomposition of squared Euclidean distance:

```math
\\| \\mathbf{x} - \\mathbf{y} \\|^2 = \\| \\mathbf{x} \\|^2 + \\| \\mathbf{y} \\|^2 - 2 \\mathbf{x}^T \\mathbf{y}
```

2. Perform eigenvalue decomposition of the Gram matrix to derive the coordinates.

*Note:*  The Gramian derived from ``D`` may have non-positive or degenerate
eigenvalues.  The subspace of non-positive eigenvalues is projected out
of the MDS solution so that the strain function is minimized in a
least-squares sense.  If the smallest remaining eigenvalue that is used
for the MDS is degenerate, then the solution is not unique, as any
linear combination of degenerate eigenvectors will also yield a MDS
solution with the same strain value.
"""
struct MDS{T<:Real} <: NonlinearDimensionalityReduction
    d::Real                  # original dimension
    X::AbstractMatrix{T}     # fitted data, X (d x n)
    λ::AbstractVector{T}     # eigenvalues in feature space, (k x 1)
    U::AbstractMatrix{T}     # eigenvectors in feature space, U (n x k)
end

## properties
"""
    size(M)

Returns tuple where the first value is the MDS model `M` input dimension,
*i.e* the dimension of the observation space, and the second value is the output
dimension, *i.e* the dimension of the embedding.
"""
size(M::MDS) = (M.d, size(M.U,2))

"""
    projection(M)

Get the MDS model `M` eigenvectors matrix (of size ``(n, p)``) of the embedding space.
The eigenvectors are arranged in descending order of the corresponding eigenvalues.
"""
projection(M::MDS) = M.U

"""
    eigvecs(M::MDS)

Get the MDS model `M` eigenvectors matrix.
"""
eigvecs(M::MDS) = projection(M)

"""
    eigvals(M::MDS)

Get the eigenvalues of the MDS model `M`.
"""
eigvals(M::MDS) = M.λ

"""
    loadings(M::MDS)

Get the loading of the MDS model `M`.
"""
loadings(M::MDS) = sqrt.(M.λ)' .* M.U

## use

"""
    predict(M, x::AbstractVector)

Calculate the out-of-sample transformation of the observation `x` for the MDS model `M`.
Here, `x` is a vector of length `d`.
"""
function predict(M::MDS, x::AbstractVector{T}; distances=false) where {T<:Real}
    d = if isnan(M.d) # model has only distance matrix
        @assert distances "Cannot transform points if model was fitted with a distance matrix. Use point distances."
        size(x, 1) != size(M.X, 1) && throw(
            DimensionMismatch("Point distances should be calculated to all original points"))
        x
    else
        if distances
            size(x, 1) != size(M.X, 2) && throw(
                DimensionMismatch("Point distances should be calculated to all original points."))
            x
        else
            size(x, 1) != size(M.X, 1) && throw(
                DimensionMismatch("Points and original data must have same dimensionality."))
            pairwise((x,y)->norm(x-y), M.X, x)
        end
    end

    # get distance matrix
    D = isnan(M.d) ? M.X : pairwise((x,y)->norm(x-y), eachcol(M.X), symmetric=true)
    d = d.^2

    # b = 0.5*(ones(n,n)*d./n - d + D*ones(n,1)./n - ones(n,n)*D*ones(n,1)./n^2)
    mD = mean(D.^2, dims=2)
    b = (d  .- mean(d, dims=1) .- mD .+ mean(mD)) ./ -2

    # sqrt(λ)⁻¹U'b
    λ = vcat(M.λ, zeros(T, size(M)[2] - length(M.λ)))
    return M.U' * b ./ sqrt.(λ)
end

"""
    predict(M)

Returns a coordinate matrix of size ``(p, n)`` for the MDS model `M`, where each column
is the coordinates for an observation in the embedding space.
"""
function predict(M::MDS{T}) where {T<:Real}
    d, p = size(M)
    # if there are non-positive missing eigval then pad with zeros
    λ = vcat(M.λ, zeros(T, p - length(M.λ)))
    return diagm(0=>sqrt.(λ)) * M.U'
end

## show

function show(io::IO, M::MDS)
    d, p = size(M)
    print(io, "Classical MDS(indim = $d, outdim = $p)")
end

## interface functions

"""
    fit(MDS, X; kwargs...)

Compute an embedding of `X` points by classical multidimensional scaling (MDS).
There are two calling options, specified via the required keyword argument `distances`:

    mds = fit(MDS, X; distances=false, maxoutdim=size(X,1)-1)

where `X` is the data matrix. Distances between pairs of columns of `X` are computed using the Euclidean norm.
This is equivalent to performing PCA on `X`.

    mds = fit(MDS, D; distances=true, maxoutdim=size(D,1)-1)

where `D` is a symmetric matrix `D` of distances between points.
"""
function fit(::Type{MDS}, X::AbstractMatrix{T};
             maxoutdim::Int = size(X,1)-1,
             distances::Bool) where T<:Real

    # get distance matrix and space dimension
    D, d = if !distances
        pairwise((x,y)->norm(x-y), eachcol(X), symmetric=true), size(X,1)
    else
        X, NaN
    end
    G = dmat2gram(D)

    n = size(D, 1)
    m = min(maxoutdim, n) #Actual number of eigenpairs wanted

    E = eigen!(Hermitian(G))

    #Sometimes dmat2gram produces a negative definite matrix, and the sign just
    #needs to be flipped. The heuristic to check for this robustly is to check
    #if there is a negative eigenvalue of magnitude larger than the largest
    #positive eigenvalue, and flip the sign of eigenvalues if necessary.
    mineig, maxeig = extrema(E.values)
    if mineig < 0 && abs(mineig) > abs(maxeig)
        #do flip
        ord = sortperm(E.values)
        λ = -E.values[ord[1:m]]
    else
        ord = sortperm(E.values; rev=true)
        λ = E.values[ord[1:m]]
    end

    for i in 1:m
        if λ[i] <= 0
            #Keeping all remaining eigenpairs would not change solution (if 0)
            #or make the answer _worse_ (if <0).
            #The least squares solution would want to throw all these away.
            @warn("Gramian has only $(i-1) positive eigenvalue(s)")
            m = i-1
            ord = ord[1:m]
            λ = λ[1:m]
            break
        end
    end

    #Check if the last considered eigenvalue is degenerate
    if m>0
        nevalsmore = sum(abs.(E.values[ord[m+1:end]] .- λ[m]) .< n*eps())
        nevals = sum(abs.(E.values .- λ[m]) .< n*eps())
        if nevalsmore > 1
            @warn("The last eigenpair is degenerate with $(nevals-1) others; $nevalsmore were ignored. Answer is not unique")
        end
    end
    U = E.vectors[:, ord[1:m]]

    #Add trailing zero coordinates if dimension of embedding space (p) exceeds
    #number of eigenpairs used (m)
    if m < maxoutdim
        U = [U zeros(T, n, maxoutdim-m)]
    end

    return MDS(d, X, λ, U)
end

"""
    stress(M)

Get the stress of the MDS mode `M`.
"""
function stress(M::MDS)
    # calculate distances if original data was stored
    DX = isnan(M.d) ? M.X : pairwise((x,y)->norm(x-y), eachcol(M.X), symmetric=true)
    DY = pairwise((x,y)->norm(x-y), eachcol(predict(M)), symmetric=true)
    n = size(DX,1)
    return sqrt(2*sum((DX - DY).^2)/sum(DX.^2));
end


text:
using JuMP
using SCIP
model = Model(SCIP.Optimizer)
# set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ON)
# set_optimizer_attribute(model, "maxNodes", 5)
month_num = 8
product_num = 1
month = 1:month_num
product = 1:product_num
# paramerters
q = [3000] # product fix cost
p = [60] # cost per product
d = [200, 200, 200, 200, 200, 200, 500, 500]' # demand of product by month
s_init = [0] # initial stock
h = [3] # stock cost per product month

@variable(model, x[i in product, j in month] >= 0)
@variable(model, s[i in product, j in month] >= 0)
@variable(model, y[i in product, j in month], Bin)

# 库存平衡
@constraint(model ,dem_sat1[i in product], s_init[i] + x[i, 1] == d[i, 1] + s[i, 1])
@constraint(model, dem_sat[i in product, j in month[2:end]], s[i, j - 1] + x[i, j] == d[i, j] + s[i, j] )
# 最大产量约束
@constraint(model, vub[i in product, j in month], x[i, j] <= sum(d[i, k]  for k in j:month_num) * y[i, j])
# 目标函数
@expression(model ,cost, sum(p[i] * x[i, j] + q[i] * y[i, j] for i in product, j in month) + sum(h[i] * s[i, j] for i in product, j in month[1:end - 1]) + sum(h[i] * s[i, month[end]] / 2 for i in product))
@objective(model, Min, cost)
optimize!(model)
# 显示结果
cost_val = value.(cost)
x_val = value.(x)
s_val = value.(s)
@show cost_val;
@show x_val;
@show s_val;

max_stars_repo_name: Behemoth-s/ProductPlanningByMIP

max_issues_repo_name: Behemoth-s/ProductPlanningByMIP

max_forks_repo_name: Behemoth-s/ProductPlanningByMIP


text:
module BiomolecularStructures
	# include sub-modules
	include(joinpath(dirname(@__FILE__), "KABSCH/", "kabsch.jl"))
	include(joinpath(dirname(@__FILE__), "PDB/", "pdb.jl"))
	include(joinpath(dirname(@__FILE__), "PLOT/", "plot.jl"))
	include(joinpath(dirname(@__FILE__), "WebBLAST/", "hit.jl"))
	include(joinpath(dirname(@__FILE__), "WebBLAST/", "WebBLAST.jl"))
	include(joinpath(dirname(@__FILE__), "MODELLER/", "modeller.jl"))
	include(joinpath(dirname(@__FILE__), "MAFFT/", "mafft.jl"))
end


text:
using Documenter: doctest
using MultilineStrings
using MultilineStrings: MultilineStrings, @m_str, indent, interpolate, multiline
using Test
using YAML: YAML

function yaml_block(str, block_scalar)
    yaml = "example: $block_scalar\n$(indent(str, 2))"
    YAML.load(yaml)["example"]
end

const TEST_STRINGS = [
    # Modified example from: https://yaml-multiline.info/
    "yaml-multiline" => """
        Several lines of text,
        with some "quotes" of various 'types'
        and also two blank lines:

        plus another line at the end.

        """,
    "whitespace" => "a\nb\n c \nd\n",
    "no ending newline" => "foo",
    "starting newline" => "\nbar",
]

# Validate `yaml_block` function
for (test, str) in TEST_STRINGS
    @assert yaml_block(str, "|+") == str
end

@testset "MultilineStrings.jl" begin
    @testset "multiline" begin
        @testset "string: $test" for (test, str) in TEST_STRINGS
            expected_lk = yaml_block(str, "|+")
            expected_lc = yaml_block(str, "|")
            expected_ls = yaml_block(str, "|-")

            expected_fk = yaml_block(str, ">+")
            expected_fc = yaml_block(str, ">")
            expected_fs = yaml_block(str, ">-")

            @testset "literal" begin
                @test multiline(str, :literal, :keep) == expected_lk
                @test multiline(str, :literal, :clip) == expected_lc
                @test multiline(str, :literal, :strip) == expected_ls

                @test multiline(str, style=:literal, chomp=:keep) == expected_lk
                @test multiline(str, style=:literal, chomp=:clip) == expected_lc
                @test multiline(str, style=:literal, chomp=:strip) == expected_ls

                @test multiline(str, "lk") == expected_lk
                @test multiline(str, "lc") == expected_lc
                @test multiline(str, "ls") == expected_ls

                @test multiline(str, "|+") == expected_lk
                @test multiline(str, "|-") == expected_ls
            end

            @testset "folding" begin
                @test multiline(str, :folded, :keep) == expected_fk
                @test multiline(str, :folded, :clip) == expected_fc
                @test multiline(str, :folded, :strip) == expected_fs

                @test multiline(str, style=:folded, chomp=:keep) == expected_fk
                @test multiline(str, style=:folded, chomp=:clip) == expected_fc
                @test multiline(str, style=:folded, chomp=:strip) == expected_fs

                @test multiline(str, "fk") == expected_fk
                @test multiline(str, "fc") == expected_fc
                @test multiline(str, "fs") == expected_fs

                @test multiline(str, ">+") == expected_fk
                @test multiline(str, ">-") == expected_fs
            end

            @testset "default chomp" begin
                @test multiline(str, style=:literal) == expected_ls
                @test multiline(str, style=:folded) == expected_fs

                @test multiline(str, "l") == expected_ls
                @test multiline(str, "f") == expected_fs

                @test multiline(str, "|") == expected_lc
                @test multiline(str, ">") == expected_fc
            end

            @testset "default style" begin
                @test multiline(str, chomp=:keep) == expected_fk
                @test multiline(str, chomp=:clip) == expected_fc
                @test multiline(str, chomp=:strip) == expected_fs

                @test multiline(str, "k") == expected_fk
                @test multiline(str, "c") == expected_fc
                @test multiline(str, "s") == expected_fs

                @test multiline(str, "+") == expected_fk
                @test multiline(str, "-") == expected_fs
            end

            @testset "default style/chomp" begin
                @test multiline(str) == expected_fs
                @test multiline(str, "") == expected_fs
            end
        end

        @testset "invalid indicators" begin
            @test_throws ArgumentError multiline("", "fs_")  # Too many indicators
            @test_throws ArgumentError multiline("", "sf")   # Order matters
            @test_throws ArgumentError multiline("", "_s")   # Invalid style
            @test_throws ArgumentError multiline("", "f_")   # Invalid chomp
            @test_throws ArgumentError multiline("", "_")    # Invalid style/chomp
        end
    end

    @testset "interpolate" begin
        @test interpolate("x") == "x"
        @test interpolate("\$x") == Expr(:string, :x)
        @test interpolate("\$(x)") == Expr(:string, :x)
        @test interpolate("\$(\"x\")") == Expr(:string, "x")

        @test interpolate("<\$x>") == Expr(:string, "<", :x, ">")
        @test interpolate("<\$(x)>") == Expr(:string, "<", :x, ">")

        # The quoting in these examples can result in exceptions being raised during parsing
        # if handled incorrectly.
        @test interpolate("\"\\n\"") == "\"\\n\""
        @test interpolate("\$(join([\"a\", \"b\"], \", \"))") == Expr(:string, :(join(["a", "b"], ", ")))
    end

    @testset "@m_str" begin
        @testset "quoting" begin
            # Use of double-quotes could cause failure if not handled properly:
            # `syntax: incomplete: invalid string syntax`
            @test m"\"\n\"" == "\" \""
        end

        @testset "string-interpolation" begin
            # If processing would accidentally take place in the interpolated code then
            # we could see "a b " as the result.
            @test m"""$(join(("a", "b") .* "\n", ""))"""fc == "a b\n"
        end
    end

    include("indent.jl")

    doctest(MultilineStrings)
end


text:
## Prepare layers & data
import Pkg
Pkg.activate(".")
@time include(joinpath("..", "..", "..", "src", "required.jl")) # best way for Intellisense
figdir = joinpath("docs", "qcbs", "2020", "fig")

# Load temperature
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
temp_full = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
# temp5 = SimpleSDMPredictor(WorldClim, BioClim, 1; resolution = 5.0, coords...)

# Clip to QC
coords_qc = (left = -82.5, right = -51.0,
             bottom = 41.0, top = 71.0)
temp_qc = temp_full[coords_qc]
# Verify dimensions
plot(temp_qc)
# Set Greenland values to nothing
xmin = SimpleSDMLayers._match_longitude(temp_qc, -55.0)
xmax = SimpleSDMLayers._match_longitude(temp_qc, coords_qc.right)
ymin = SimpleSDMLayers._match_latitude(temp_qc, 60.0)
ymax = SimpleSDMLayers._match_latitude(temp_qc, coords_qc.top)
temp_qc.grid[ymin:ymax, xmin:xmax] .= nothing
plot(temp_qc)

## Create maps for QC
qc_inferno = plot(temp_qc)
qc_lightgrey = plot(temp_qc, c = :lightgrey)
qc_viridis = plot(temp_qc, c = :viridis)
qc_white = plot(temp_qc, c = :white)

## Remove background
# Custom function
function removebackground!(p; kw...)
    plot!(p,
          frame = :origin,
          colorbar = :none,
          guides = "",
          grid = :none,
          bg = :transparent,
          bg_inside = :transparent;
          kw...
          )
    return p
end

# Remove background
removebackground!(qc_inferno)
removebackground!(qc_lightgrey)
removebackground!(qc_viridis)
removebackground!(qc_white)

# Save results
savefig(plot(qc_inferno,   dpi = 200), joinpath(figdir, "qc_no-bg_inferno.png"))
savefig(plot(qc_lightgrey, dpi = 200), joinpath(figdir, "qc_no-bg_lightgrey.png"))
savefig(plot(qc_viridis,   dpi = 200), joinpath(figdir, "qc_no-bg_viridis.png"))
savefig(plot(qc_white,     dpi = 200), joinpath(figdir, "qc_no-bg_white.png"))

# Paint it white
plot!(qc_inferno, bg_inside = :white)
plot!(qc_lightgrey, bg_inside = :white)
plot!(qc_viridis, bg_inside = :white)

# Save results
savefig(plot(qc_inferno,   dpi = 200), joinpath(figdir, "qc_inside-bg_inferno.png"))
savefig(plot(qc_lightgrey, dpi = 200), joinpath(figdir, "qc_inside-bg_lightgrey.png"))
savefig(plot(qc_viridis,   dpi = 200), joinpath(figdir, "qc_inside-bg_viridis.png"))

## Draw rectangle over subareas
# Get coordinates
coords_NE = (left = -80.0, right = -60.0, bottom = 40.0, top = 50.0)
coords_SW = (left = -120.0, right = -100.0, bottom = 30.0, top = 40.0)

# Rectangle coordinates
function rectangle_from_coords(xb,yb,xt,yt)
    [
        xb  yb
        xt  yb
        xt  yt
        xb  yt
        xb  yb
        NaN NaN
    ]
end
rect_NE = rectangle_from_coords(coords_NE.left,  coords_NE.bottom,
                                coords_NE.right, coords_NE.top)
rect_SW = rectangle_from_coords(coords_SW.left,  coords_SW.bottom,
                                coords_SW.right, coords_SW.top)

# Plot subarea rectangles
rectplot = plot(temp_full, c = :lightgrey, cb = :none)
plot!(rectplot, rect_NE[:, 1], rect_NE[:, 2], label = "NE subarea")
plot!(rectplot, rect_SW[:, 1], rect_SW[:, 2], label = "SW subarea")
plot!(rectplot, bg_outside = :transparent)
savefig(plot(rectplot, dpi = 200), joinpath(figdir, "subarea_map.png"))

## Remove background on analysis figures
# Predefine desired outcome
outcome = "bart"
# Load analysis script
include(joinpath("..", "..", "..", "src", "04_full-extent.jl"))

# Remove backgrounds
plot!(richness_plot, bg_outside = :transparent, title = "")
plot!(lcbdtr_plot, bg_outside = :transparent, title = "")
plot!(rel2d_plot, bg_outside = :transparent, title = "")

# Save results
savefig(plot(richness_plot, dpi = 200), joinpath(figdir, "res_richness.png"))
savefig(plot(lcbdtr_plot, dpi = 200), joinpath(figdir, "res_lcbd.png"))
savefig(plot(rel2d_plot, dpi = 200), joinpath(figdir, "res_relationship.png"))

## Distribution figures
# Single species distributions
spplot1 = plot_layer(distributions[1][coords_qc], c = :BuPu)
spplot2 = plot_layer(distributions[17][coords_qc], c = :BuPu)
spplot3 = plot_layer(distributions[22][coords_qc], c = :BuPu)
# Remove backgrounds
removebackground!(spplot1)
removebackground!(spplot2)
removebackground!(spplot3)
# Save
savefig(plot(spplot1, dpi = 200), joinpath(figdir, "sp_bart_1.png"))
savefig(plot(spplot2, dpi = 200), joinpath(figdir, "sp_bart_2.png"))
savefig(plot(spplot3, dpi = 200), joinpath(figdir, "sp_bart_3.png"))

## Raw distribution figures
# Load data
distributions_bart = copy(distributions)
@load joinpath("data", "jld2", "raw-distributions.jld2") distributions
distributions_raw = copy(distributions)
distributions = distributions_bart
# Plot
spplot1_raw = plot_layer(distributions_raw[1][coords_qc], c = :BuPu)
spplot2_raw = plot_layer(distributions_raw[17][coords_qc], c = :BuPu)
spplot3_raw = plot_layer(distributions_raw[22][coords_qc], c = :BuPu)
# Remove backgrounds
removebackground!(spplot1_raw)
removebackground!(spplot2_raw)
removebackground!(spplot3_raw)
# Save
savefig(plot(spplot1_raw, dpi = 200), joinpath(figdir, "sp_raw_1.png"))
savefig(plot(spplot2_raw, dpi = 200), joinpath(figdir, "sp_raw_2.png"))
savefig(plot(spplot3_raw, dpi = 200), joinpath(figdir, "sp_raw_3.png"))

## Subarea comparison figure
include(joinpath("..", "..", "..", "src", "05_subareas.jl"))
plot!(combined_plot, bg_outside = :transparent)
savefig(plot(combined_plot, dpi = 200), joinpath(figdir, "subarea_comparison.png"))

## Relationship plot
@load joinpath("data", "jld2", "raw-distributions.jld2") distributions
outcome = "raw"
include(joinpath("..", "..", "..", "src", "04_full-extent.jl"))
plot(richness, lcbd, colour = :transparent, smooth = true)

# Arrange data
resdf = DataFrame([richness, lcbd])
rename!(resdf, :x1 => :richness, :x2 => :lcbd)
filter!(x -> !isnothing(x.richness) && !isnothing(x.lcbd), resdf)
resdf.richness = Array{Float64}(resdf.richness)
resdf.lcbd = Array{Float64}(resdf.lcbd)

# Fit loess model
using Loess
function loess_vals(x, y; kw...)
    model = loess(x, y; kw...)
    xloess = range(extrema(x)...; step = 0.1)
    yloess = Loess.predict(model, xloess)
    return xloess, yloess
end
# Plot
xs, ys = loess_vals(resdf.richness, resdf.lcbd)
plot(xs, ys)

# Blank relationship plot
x = 0:0.1:10.0
y = @__dot__ 0.001*x^2 - 0.041*x - 3.394 # Heino 2017
y = @__dot__ 0.002*x^2 - 0.041*x - 3.394 # Heino 2017
y = @__dot__ 0.025*x^2 - 0.442*x - 2.690 # Silva 2018
relplot = plot(x, y,
               c = :black,
               legend = false,
               ylabel = "LCBD",
               xlabel = "Richness",
               )
plot!(labelfontsize = 16,
      ylim = (-3.62, Inf),
      ticks = false,
      thickness_scaling = 2)
# plot!(xlim = (-0.25, 10.0),
#       ylim = (-4.75, Inf),
#       aspect_ratio = 4.5
#       ) # for Silva 2018
savefig(plot(relplot, dpi = 50, bg = :transparent),
        joinpath(figdir, "relationship-plot.png"))

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
# This file was generated, do not modify it. # hide
using PyPlot
ioff() # hide
figure(figsize=(8,6))

ym1 = y .== -1
ym2 = .!ym1
plot(X[ym1, 1], X[ym1, 2], ls="none", marker="o")
plot(X[ym2, 1], X[ym2, 2], ls="none", marker="x")

savefig(joinpath(@OUTPUT, "ISL-lab-9-g1.svg")) # hide


text:
export linearsolve, LinSolver, Precond, LinSolveOpts, solvesystem

using JOLI

const LinSolver = Set([:fgmres,:lufact])
const Precond = Set([:vgmres,:mlgmres,:identity])

mutable struct LinSolveOpts
    tol::Float64
    maxit::Int64
    maxinnerit::Int64
    solver::Symbol
    #precond::Union{Symbol,Function,LinSolveOpts,joAbstractOperator}
    precond::Any
    outputfreq::Int64
    tag::String
    function LinSolveOpts(;tol::Float64=1e-6,
                          maxit::Int64=10000,
                          maxinnerit::Int64=20,
                          solver::Symbol=:fgmres,
                          precond::Union{Symbol,Function}=:identity,
                          outputfreq::Int64=0,
                          tag::String="")
        solver in LinSolver || throw(ArgumentError("Invalid solver $(solver), valid options are $(LinSolver)"))
        if typeof(precond)==Symbol
            precond in Precond || throw(ArgumentError("Invalid preconditioner: $(precond)"))
        end
        new(tol,maxit,maxinnerit,solver,precond,outputfreq,tag)
    end
end

# Generic linear solution interface
#
function linearsolve(op,b,x0,lsopts::LinSolveOpts;forw_mode::Bool=true)
    if length(x0)==0
        x0 = similar(b) .* 0.
    end
    size(op,1)==size(op,2) || throw(ArgumentError("A must be square"))
    size(b,1)==size(op,1) || throw(ArgumentError("b has $(size(b,1)) elements, op is $(size(op,1)) x $(size(op,2))"))
    size(x0,1)==size(op,1) || throw(ArgumentError("x0 dimension mismatch"))
    if forw_mode
        A = op
    else
        A = op'
    end

    if typeof(lsopts.precond)==LinSolveOpts
        prec_opts = lsopts.precond

        if prec_opts.solver==:fgmres
            P = x->FGMRES(A,x,similar(x) .* 0.,m=prec_opts.maxinnerit,maxit=prec_opts.maxit,tol=prec_opts.tol)[1]
        end
    elseif typeof(lsopts.precond)<:Function
        P = x->lsopts.precond(x,similar(x) .* 0.,forw_mode)
    elseif typeof(lsopts.precond)==Symbol
        if lsopts.precond==:identity
            P = nothing
        else
            throw(ArgumentError("Unrecognized preconditioner $(lsopts.precond)"))
        end
    elseif typeof(lsopts.precond)<:joAbstractOperator
        if forw_mode
            P = x->lsopts.precond*x
        else
            P = x->lsopts.precond'*x
        end
    else
        P = nothing
    end

    if lsopts.solver==:fgmres
        (y,res) = FGMRES(A,b,x0,m=lsopts.maxinnerit,maxiter=lsopts.maxit,tol=lsopts.tol,precond=P,outputfreq=lsopts.outputfreq)
        return y
    elseif lsopts.solver==:lufact
        y = A\b
        return y
    else
        throw(ArgumentError("Unrecognized solver $(lsopts.solver)"))
    end
end

function solvesystem(H,lsopts::LinSolveOpts)
    return (b,x,mode::Bool)->linearsolve(H,b,x,lsopts,forw_mode=mode)
end


text:
mutable struct Bag{T <: Vector{<:AbstractPiece}} <: AbstractBag
  player::AbstractPlayer
  pieces::T
  previews::Int
end

function Bag(cur_player::AbstractPlayer=Player())
  cur_bag = Bag(
    cur_player,
    AbstractPiece[],
    3
  )

  cur_player.bag = cur_bag

  cur_bag
end


text:
mutable struct Round{T <: Vector{<:AbstractLog}} <: AbstractRound
  player::AbstractPlayer
  logs::T

  batch_count::Int
  log_count::Int

  is_keeping_score::Bool
  is_making_call::Bool

  server_id::Int
  client_password::AbstractString
end

function Round(cur_player::AbstractPlayer=Player())
  cur_password = string(UUIDs.uuid4())

  cur_batch_count = 0
  cur_log_count = 1

  cur_request = Dict(
    "round" => Dict(
      "password" => cur_password,
      "batch_count" => cur_batch_count
    )
  )

  cur_is_keeping_score = false

  cur_response = nothing

  try
    cur_response = HTTP.post(
      rounds_endpoint,
      json_headers,
      JSON.json(cur_request)
    )

    cur_is_keeping_score = ( cur_response.status == 201 )
  catch cur_error
    any(
      cur_error_type -> isa(cur_error, cur_error_type),
      api_errors
    ) || rethrow(cur_error)

    tetris_js(
      cur_player.game.scope,
      JSString("""
        \$(".js-offline").removeClass("hidden");
      """)
    )
  end

  cur_server_id = cur_is_keeping_score ?
    api_fetch(cur_response, "id") : -1

  cur_round = Round(
    cur_player,
    AbstractLog[],
    cur_batch_count,
    cur_log_count,
    cur_is_keeping_score,
    false,
    cur_server_id,
    cur_password
  )

  cur_player.round = cur_round

  cur_round
end


text:
using Distributions
using Base.Test

for d in (Truncated(Normal(0, 1), -1, 1),
          Truncated(Normal(3, 10), 7, 8),
          Truncated(Normal(-5, 1), -Inf, -10))
    @test all(insupport(d, rand(d, 1000)))
end

d = Truncated(Normal(0, 1), -0.1, +0.1)

@test pdf(d, 0.0) > pdf(Normal(0, 1), 0.0)
@test pdf(d, -1.0) == 0.0
@test pdf(d, +1.0) == 0.0

@test logpdf(d, 0.0) > logpdf(Normal(0, 1), 0.0)
@test isinf(logpdf(d, -1.0))
@test isinf(logpdf(d, +1.0))

@test cdf(d, -1.0) == 0.0
@test cdf(d, -0.09) < cdf(Normal(0, 1), -0.09)
@test cdf(d, 0.0) == 0.5
@test cdf(d, +0.09) > cdf(Normal(0, 1), +0.09)
@test cdf(d, +1.0) == 1.0

@test quantile(d, 0.01) > -0.1
@test abs(quantile(d, 0.5) - 0.0) < 1e-8
@test quantile(d, 0.99) < +0.1

@test abs(cdf(d, quantile(d, 0.01)) - 0.01) < 1e-8
@test abs(cdf(d, quantile(d, 0.50)) - 0.50) < 1e-8
@test abs(cdf(d, quantile(d, 0.99)) - 0.99) < 1e-8


text:
### reinterpret

@pure samesize(::Type{T}, ::Type{S}) where {T,S} = sizeof(T) == sizeof(S)

## Color->Color
# function reinterpretc(::Type{CV1}, a::Array{CV2,1}) where {CV1<:Colorant,CV2<:Colorant}
#     CV = ccolor(CV1, CV2)
#     l = (length(a)*sizeof(CV2))÷sizeof(CV1)
#     l*sizeof(CV1) == length(a)*sizeof(CV2) || throw(ArgumentError("sizes are incommensurate"))
#     reshape(reinterpret(CV, a), (l,))
# end
function reinterpretc(::Type{CV1}, a::AbstractArray{CV2}) where {CV1<:Colorant,CV2<:Colorant}
    CV = ccolor(CV1, CV2)
    if samesize(CV, CV2)
        return reshape(reinterpret(CV, a), size(a))
    end
    throw(ArgumentError("result shape not specified"))
end

## Color->T
function reinterpretc(::Type{T}, a::AbstractArray{CV}) where {T<:Number,CV<:Colorant}
    if samesize(T, CV)
        return reinterpret(T, a)
    end
    axs = axes(a)
    if sizeof(CV) == sizeof(T)*_len(CV)
        return reinterpret(T, reshape(a, Base.OneTo(1), axs...))
    end
    throw(ArgumentError("result shape not specified"))
end
reinterpretc(::Type{T}, a::AbstractArray{CV,0}) where {T<:Number,CV<:Colorant} =
    reinterpret(T, reshape(a, 1))

_len(::Type{C}) where {C} = _len(C, eltype(C))
_len(::Type{C}, ::Type{Any}) where {C} = error("indeterminate type")
_len(::Type{C}, ::Type{T}) where {C,T} = sizeof(C) ÷ sizeof(T)

## T->Color
# We have to distinguish two forms of call:
#   form 1: reinterpretc(RGB{N0f8}, img)
#   form 2: reinterpretc(RGB, img)
function reinterpretc(::Type{CV}, a::AbstractArray{T}) where {CV<:Colorant,T<:Number}
    @noinline throwdm(::Type{C}, ind1) where C =
        throw(DimensionMismatch("indices $ind1 are not consistent with color type $C"))
    CVT = ccolor_number(CV, T)
    if samesize(CVT, T)
        return reinterpret(CVT, a)
    end
    axs = axes(a)
    if axs[1] == Base.OneTo(sizeof(CVT) ÷ sizeof(eltype(CVT)))
        return reshape(reinterpret(CVT, a), tail(axs))
    end
    throwdm(CV, axs[1])
end

# ccolor_number converts form 2 calls to form 1 calls
ccolor_number(::Type{CV}, ::Type{T}) where {CV<:Colorant,T<:Number} =
    ccolor_number(CV, eltype(CV), T)
ccolor_number(::Type{CV}, ::Type{CVT}, ::Type{T}) where {CV,CVT<:Number,T} = CV # form 1
ccolor_number(::Type{CV}, ::Type{Any}, ::Type{T}) where {CV<:Colorant,T} = CV{T} # form 2

### convert
#
# The main contribution here is "concretizing" the colorant type: allow
#    convert(RGB, a)
# rather than requiring
#    convert(RGB{N0f8}, a)
# Where possible the raw element type of the source is retained.
Base.convert(::Type{Array{C}},   img::Array{C,n}) where {C<:Color1,n} = img
Base.convert(::Type{Array{C,n}}, img::Array{C,n}) where {C<:Color1,n} = img
Base.convert(::Type{Array{C}},   img::Array{C,n}) where {C<:Colorant,n} = img
Base.convert(::Type{Array{C,n}}, img::Array{C,n}) where {C<:Colorant,n} = img

function Base.convert(::Type{Array{Cdest}},
                      img::AbstractArray{Csrc,n}) where {Cdest<:Colorant,n,Csrc<:Colorant}
    convert(Array{Cdest,n}, img)
end

function Base.convert(::Type{Array{Cdest,n}},
                      img::AbstractArray{Csrc,n}) where {Cdest<:Colorant,n,Csrc<:Colorant}
    copyto!(Array{ccolor(Cdest, Csrc)}(undef, size(img)), img)
end

function Base.convert(::Type{Array{Cdest}},
                      img::BitArray{n}) where {Cdest<:Color1,n}
    convert(Array{Cdest,n}, img)
end

function Base.convert(::Type{Array{Cdest}},
                      img::AbstractArray{T,n}) where {Cdest<:Color1,n,T<:Real}
    convert(Array{Cdest,n}, img)
end

function Base.convert(::Type{Array{Cdest,n}},
                      img::BitArray{n}) where {Cdest<:Color1,n}
    copyto!(Array{ccolor(Cdest, Gray{Bool})}(undef, size(img)), img)
end

function Base.convert(::Type{Array{Cdest,n}},
                      img::AbstractArray{T,n}) where {Cdest<:Color1,n,T<:Real}
    copyto!(Array{ccolor(Cdest, Gray{T})}(undef, size(img)), img)
end

function Base.convert(::Type{OffsetArray{Cdest,n,A}},
                      img::AbstractArray{Csrc,n}) where {Cdest<:Colorant,n, A <:AbstractArray,Csrc<:Colorant}
    copyto!(OffsetArray{ccolor(Cdest, Csrc)}(undef, axes(img)),img)
end

function Base.convert(::Type{OffsetArray{C,n,A}},
                      img::AbstractArray{C,n}) where {C<:Colorant,n, A <:AbstractArray}
    img
end

# for docstrings in the operations below
shortname(::Type{T}) where {T<:FixedPoint} = (io = IOBuffer(); FixedPointNumbers.showtype(io, T); String(take!(io)))
shortname(::Type{T}) where {T} = string(T)

# float32, float64, etc. Used for conversions like
#     Array{RGB{N0f8}} -> Array{RGB{Float32}},
# since
#    convert(Array{RGB{Float32}}, A)
# is annoyingly verbose for such a common operation.
for (fn,T) in (#(:float16, Float16),   # Float16 currently has promotion problems
               (:float32, Float32), (:float64, Float64),
               (:n0f8, N0f8), (:n6f10, N6f10),
               (:n4f12, N4f12), (:n2f14, N2f14), (:n0f16, N0f16))
    @eval begin
        ($fn)(::Type{C}) where {C<:Colorant} = base_colorant_type(C){$T}
        ($fn)(::Type{S}) where {S<:Number  } = $T
        ($fn)(c::Colorant) = convert(($fn)(typeof(c)), c)
        ($fn)(n::Number)   = convert(($fn)(typeof(n)), n)
        @deprecate ($fn)(A::AbstractArray{C}) where {C<:Colorant} ($fn).(A)
        fname = $(Expr(:quote, fn))
        Tname = shortname($T)
@doc """
    $fname.(img)

converts the raw storage type of `img` to `$Tname`, without changing the color space.
""" $fn

    end
end

"""
    float(x::Colorant)
    float(T::Type{<:Colorant})

convert the storage type of pixel `x` to a floating point data type while
preserving the `Colorant` information.

If the input is Type `T`, then it is equivalent to [`floattype`](@ref).
"""
float(x::Colorant) = floattype(typeof(x))(x)
float(::Type{T}) where T <: Colorant = floattype(T)


text:
# Rounding
function round(x::Decimal; digits::Int=0, normal::Bool=false)
    shift = BigInt(digits) + x.q
    if shift > BigInt(0) || shift < x.q
        (normal) ? x : normalize(x, rounded=true)
    else
        c = Base.round(x.c / BigInt(10)^(-shift))
        d = Decimal(x.s, BigInt(c), x.q - shift)
        (normal) ? d : normalize(d, rounded=true)
    end
end

function trunc(x::Decimal; digits::Int=0, normal::Bool=false)
    shift = BigInt(digits) + x.q
    if shift > BigInt(0) || shift < x.q
        (normal) ? x : normalize(x, rounded=true)
    else
        c = Base.trunc(x.c / BigInt(10)^(-shift))
        d = Decimal(x.s, BigInt(c), x.q - shift)
        (normal) ? d : normalize(d, rounded=true)
    end
end


text:
export count

import Base: count
import Base: show

"""
    count()

Creates a count operator, which counts the number of
emissions on the source and emits that number when the source completes.

# Producing

Stream of type `<: Subscribable{Int}`

# Examples
```jldoctest
using Rocket

source = from([ i for i in 1:42 ])
subscribe!(source |> count(), logger())
;

# output

[LogActor] Data: 42
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`RightTypedOperator`](@ref), [`ProxyObservable`](@ref), [`logger`](@ref)
"""
count() = CountOperator()

struct CountOperator <: RightTypedOperator{Int} end

function on_call!(::Type{L}, ::Type{Int}, operator::CountOperator, source) where L
    return proxy(Int, source, CountProxy())
end

struct CountProxy <: ActorProxy end

actor_proxy!(::Type, proxy::CountProxy, actor::A) where A = CountActor{A}(actor)

mutable struct CountActor{A} <: Actor{Any}
    actor   :: A
    current :: Int

    CountActor{A}(actor::A) where A = new(actor, 0)
end

on_next!(actor::CountActor, data) = begin actor.current += 1 end
on_error!(actor::CountActor, err) = begin error!(actor.actor, err) end
on_complete!(actor::CountActor)   = begin next!(actor.actor, actor.current); complete!(actor.actor) end

Base.show(io::IO, ::CountOperator) = print(io, "CountOperator()")
Base.show(io::IO, ::CountProxy)    = print(io, "CountProxy()")
Base.show(io::IO, ::CountActor)    = print(io, "CountActor()")


text:
using DomainSets
import BasisFunctions: resize
## Functions to provide oversampled grids for FrameFun approximations.

resize(grid::ProductGrid, n) = ProductGrid(map(resize, elements(grid), n)...)

equispacedgrid(domain::AbstractInterval, L::Int) =
    EquispacedGrid(L, infimum(domain), supremum(domain))

equispacedgrid(domain::ProductDomain, L::Tuple) =
    cartesianproduct(map(equispacedgrid, elements(domain), L))

equispacedgrid(domain::Domain, L::Tuple) = subgrid(equispacedgrid(boundingbox(domain), L), domain)

function oversampling_grid(dict::Dictionary, L)
    if hasinterpolationgrid(dict)
        try
            return resize(interpolation_grid(dict), L)
        catch
            return interpolation_grid(resize(dict, L))
        end
    else
        # If the dictionary does not have an associated grid, we compute an
        # equispaced grid on its support.
        equispacedgrid(support(dict), L)
    end
end

oversampling_grid(dict::TensorProductDict, L) = ProductGrid(map(oversampling_grid, elements(dict), L)...)
oversampling_grid(dict::TensorProductDict, L::CartesianIndex) = oversampling_grid(dict, L.I)
oversampling_grid(dict::TensorProductDict, L::Int) = error("Expects a tuple or a CartesianIndex")

oversampling_grid(dict::BasisFunctions.CompositeDict, L) = oversampling_grid(element(dict,1), L)

oversampling_grid(dict::WeightedDict, L) = oversampling_grid(superdict(dict), L)

oversampling_grid(dict::MappedDict, L) = apply_map(oversampling_grid(superdict(dict), L), mapping(dict))

oversampling_grid(dict::OperatedDict, L) = oversampling_grid(superdict(dict), L)

# Make an initial guess for the sampling parameter L. The problem is that we have to
# determine the type of L. We try to infer this from the dimension of the dictionary and/or
# of its support.
initialguess(ap, M) = initialguess(dictionary(ap), M)
initialguess(dict::Dictionary1d, M::Int) = M
initialguess(dict::BasisFunctions.CompositeDict, M) = initialguess(element(dict,1), M)
initialguess(dict::BasisFunctions.CompositeDict, M::Int) = initialguess(element(dict,1), M)
initialguess(dict::Dictionary, M) = _initialguess(dict, M, support(dict), size(dict))
_initialguess(dict::Dictionary, M, domain::Domain, size) = size
_initialguess(dict::Dictionary, M, domain::Domain1d, size::Tuple{Int}) = M
function _initialguess(dict::Dictionary, M, domain::Domain2d, size::Tuple{Int})
    n = round(Int,sqrt(size[1]))
    (n,n)
end

# TODO Should everything be implemented for Dictionaries and Platforms or do we shift to platforms only
match_and_correct_sampling_parameter(strategy::SamplingStrategy, platform::Platform, param, M; dict=dictionary(platform, param), options...) =
    match_and_correct_sampling_parameter(strategy, platform, param, M, initialguess(dict, M); dict=dict, options...)
match_and_correct_sampling_parameter(::SamplingStrategy, dict::Dictionary, M; options...) =
    match_and_correct_sampling_parameter(dict, M, initialguess(dict, M); options...)

function match_and_correct_sampling_parameter(strategy::SamplingStrategy, platform::Platform, param, M, L_init; dict=dictionary(platform, param), options...)
    L_trial = match_sampling_parameter(dict, M, L_init)
    # The above does not always return a suitable L since it does not account for the platform/Dictionary
    # correct the result to a suitable one if necesarry
    correct_sampling_parameter(strategy, platform, param, L_trial; dict=dict, options...)
end

function match_and_correct_sampling_parameter(dict::Dictionary, M, L_init; options...)
    L_trial = match_sampling_parameter(dict, M, L_init)
    # The above does not always return a suitable L since it does not account for the platform/Dictionary
    # correct the result to a suitable one if necesarry
    correct_sampling_parameter(dict, L_trial; options...)
end

correct_sampling_parameter(dict::Dictionary, L_trial; options...) = L_trial
correct_sampling_parameter(strategy::SamplingStrategy, platform::Platform, param, L_trial; options...) = L_trial
correct_sampling_parameter(strategy::SamplingStrategy, p::ProductPlatform, param, L; options...) =
    tuple(map(x->correct_sampling_parameter(strategy, x...; options...), zip(elements(p), param, L))...)

function match_sampling_parameter(samplingobject::Dictionary, M, L_init)
    objective(L) = length(oversampling_grid(samplingobject, L))-M
    monotonic_bisect(objective, L_init)
end

struct BisectionType end
Base.broadcastable(bt::BisectionType) = Ref(bt)

bisect_smaller(::BisectionType, L::Int) = L>>1
bisect_smaller(::BisectionType, L::Number) = L/2
bisect_larger(::BisectionType, L::Number) = 2L

bisect_average(bt::BisectionType, L1, L2) = bisect_average(bt, promote(L1,L2)...)
bisect_average(::BisectionType, L1::Int, L2::Int) = ceil(Int, (L1+L2)/2)
bisect_average(::BisectionType, L1::T, L2::T) where {T} = convert(T, (L1+L2)/2)

bisect_smaller(bt::BisectionType, L::Tuple) = bisect_smaller.(bt, L)
bisect_larger(bt::BisectionType, L::Tuple) = bisect_larger.(bt, L)
bisect_average(bt::BisectionType, L1::Tuple, L2::Tuple) = bisect_average.(bt, L1, L2)

bisect_small(l1) = bisect_smaller(BisectionType(), l1)
bisect_larger(l1) = bisect_larger(BisectionType(), l2)
bisect_average(l1, l2) = bisect_average(BisectionType(), l1, l2)

"""
Find the root of a monotonically increasing function with a given starting value
`L_init` and using bisection.

The value is increased or decreased first (using `bisect_larger` and
`bisect_smaller`) until an interval is found for which the objective function
switches sign. Next, the smallest value of `L` for which the objective function
becomes positive is located using bisection (with `bisect_average` to average out
two values of `L`).

The user can optionally supply an argument `bt`, which will be passed as the
first argument to the `bisect_X` routines. Thus, the variable `L` can have any
type as long as making it larger/smaller or computing averages can be meaningfully
defined.

The bisection is run until convergence. This ensures that if `L` has a discrete
nature, the optimal solution is always found exactly, in finite time.
"""
function monotonic_bisect(objective, L_init; maxiterations = 20, bt = BisectionType())
    local L1
    Z = objective(L_init)
    if Z ≈ 0
        return L_init
    end
    iterations = 0
    if Z > 0
        while (Z > 0) && (iterations < maxiterations)
            iterations += 1
            L1 = L_init
            L_init = bisect_smaller(bt, L1)
            Z = objective(L_init)
        end
        iterations == maxiterations && error("Maximal number of iterations reached in monotonic_bisect.")
        monotonic_bisect(objective, L_init, L1, bt)
    else
        while (Z < 0) && (iterations <  maxiterations)
            iterations += 1
            L1 = L_init
            L_init = bisect_larger(bt, L1)
            Z = objective(L_init)
        end
        iterations == maxiterations && error("Maximal number of iterations reached in monotonic_bisect.")
        monotonic_bisect(objective, L1, L_init, bt)
    end
end

function monotonic_bisect(objective, Lmin, Lmax, bt)
    Lmin, Lmax, bt
    Lmid = bisect_average(bt, Lmin, Lmax)
    if (Lmid == Lmax) || (Lmid == Lmin)
        return Lmid
    end
    Z = objective(Lmid)
    if Z ≈ 0
        return Lmid
    elseif Z > 0
        monotonic_bisect(objective, Lmin, Lmid, bt)
    else
        monotonic_bisect(objective, Lmid, Lmax, bt)
    end
end


text:
module m1
using FinEtools
using FinEtoolsFlexBeams.FESetCorotBeamModule: FESetL2CorotBeam
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionCircle
using Test
function test()
    cs = CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
    fes = FESetL2CorotBeam([1 2; 2 3; 3 4], cs)
    @test count(fes) == 3
    true
end
end
using .m1
m1.test()

module m2
using FinEtools
using FinEtoolsFlexBeams.FESetCorotBeamModule: FESetL2CorotBeam
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionCircle
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member
using Test
function test()
    L = 42
    xyz = [0 0 0;
    0 L/4 L*1/4;
    L/4 L/4 L*2/4;
    L/4 0 L*3/4;
    0 0 L]
    nL = 20

    cs = CrossSectionCircle(s -> 5.9910, s -> [0.0, 0.0, 1.0])
    fens, fes = frame_member(xyz, nL, cs)
    @test count(fes) == 20
    true
end
end
using .m2
m2.test()


text:
function ImGui_ImplOpenGL3_CreateFontsTexture(ctx::Context)
    # build texture atlas
    fonts = unsafe_load(igGetIO().Fonts)
    pixels = Ptr{Cuchar}(C_NULL)
    width, height = Cint(0), Cint(0)
    @c ImFontAtlas_GetTexDataAsRGBA32(fonts, &pixels, &width, &height, C_NULL)

    # upload texture to graphics system
    last_texture = GLint(0)
    @c glGetIntegerv(GL_TEXTURE_BINDING_2D, &last_texture)
    ctx.FontTexture = GLuint(0)
    @c glGenTextures(1, &ctx.FontTexture)
    glBindTexture(GL_TEXTURE_2D, ctx.FontTexture)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glPixelStorei(GL_UNPACK_ROW_LENGTH, 0)
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA, width, height, 0, GL_RGBA, GL_UNSIGNED_BYTE, pixels)

    # store our identifier
    ImFontAtlas_SetTexID(fonts, ImTextureID(Int(ctx.FontTexture)))

    # restore state
    glBindTexture(GL_TEXTURE_2D, last_texture)

    return true
end

function ImGui_ImplOpenGL3_DestroyFontsTexture(ctx::Context)
    io::Ptr{ImGuiIO} = igGetIO()
    @c glDeleteTextures(1, &ctx.FontTexture)
    ImFontAtlas_SetTexID(unsafe_load(io.Fonts), ImTextureID(0))
    ctx.FontTexture = 0
    return true
end


text:
"""
    ScaledHotVector{T}

Represents a vector of at most one value different from 0.
"""
struct ScaledHotVector{T} <: AbstractVector{T}
    active_val::T
    val_idx::Int
    len::Int
end

Base.size(v::ScaledHotVector) = (v.len,)

@inline function Base.getindex(v::ScaledHotVector{T}, idx::Integer) where {T}
    @boundscheck if !(1 ≤ idx ≤ length(v))
        throw(BoundsError(v, idx))
    end
    if v.val_idx != idx
        return zero(T)
    end
    return v.active_val
end

Base.sum(v::ScaledHotVector) = v.active_val

function LinearAlgebra.dot(v1::ScaledHotVector, v2::AbstractVector)
    return v1.active_val * v2[v1.val_idx]
end

function LinearAlgebra.dot(v1::ScaledHotVector{<:Number}, v2::SparseArrays.SparseVector{<:Number})
    return v1.active_val * v2[v1.val_idx]
end

LinearAlgebra.dot(v1::AbstractVector, v2::ScaledHotVector) = dot(v2, v1)

# warning, no bound check
function LinearAlgebra.dot(v1::ScaledHotVector, v2::ScaledHotVector)
    if length(v1) != length(v2)
        throw(DimensionMismatch("v1 and v2 do not have matching sizes"))
    end
    return v1.active_val * v2.active_val * (v1.val_idx == v2.val_idx)
end

function Base.:*(v::ScaledHotVector, x::Number)
    return ScaledHotVector(v.active_val * x, v.val_idx, v.len)
end

Base.:*(x::Number, v::ScaledHotVector) = v * x

function Base.:+(x::ScaledHotVector, y::AbstractVector)
    if length(x) != length(y)
        throw(DimensionMismatch())
    end
    yc = Base.copymutable(y)
    @inbounds yc[x.val_idx] += x.active_val
    return yc
end

Base.:+(y::AbstractVector, x::ScaledHotVector) = x + y

function Base.:+(x::ScaledHotVector{T1}, y::ScaledHotVector{T2}) where {T1,T2}
    n = length(x)
    T = promote_type(T1, T2)
    if n != length(y)
        throw(DimensionMismatch())
    end
    res = spzeros(T, n)
    @inbounds res[x.val_idx] = x.active_val
    @inbounds res[y.val_idx] = y.active_val
    return res
end

Base.:-(x::ScaledHotVector{T}) where {T} = ScaledHotVector{T}(-x.active_val, x.val_idx, x.len)

Base.:-(x::AbstractVector, y::ScaledHotVector) = +(x, -y)
Base.:-(x::ScaledHotVector, y::AbstractVector) = +(x, -y)

Base.:-(x::ScaledHotVector, y::ScaledHotVector) = +(x, -y)

Base.similar(v::ScaledHotVector{T}) where {T} = spzeros(T, length(v))

function Base.convert(::Type{Vector{T}}, v::ScaledHotVector) where {T}
    vc = zeros(T, v.len)
    vc[v.val_idx] = v.active_val
    return vc
end

"""
    RankOneMatrix{T, UT, VT}

Represents a rank-one matrix `R = u * vt'`.
Composes like a charm.
"""
struct RankOneMatrix{T,UT<:AbstractVector,VT<:AbstractVector} <: AbstractMatrix{T}
    u::UT
    v::VT
end

function RankOneMatrix(u::UT, v::VT) where {UT,VT}
    T = promote_type(eltype(u), eltype(v))
    return RankOneMatrix{T,UT,VT}(u, v)
end

# not checking indices
Base.@propagate_inbounds function Base.getindex(R::RankOneMatrix, i, j)
    @boundscheck (checkbounds(R.u, i); checkbounds(R.v, j))
    @inbounds R.u[i] * R.v[j]
end

Base.size(R::RankOneMatrix) = (length(R.u), length(R.v))
function Base.:*(R::RankOneMatrix, v::AbstractVector)
    temp = fast_dot(R.v, v)
    return R.u * temp
end

function Base.:*(R::RankOneMatrix, M::AbstractMatrix)
    temp = R.v' * M
    return RankOneMatrix(u, temp')
end

function Base.:*(R1::RankOneMatrix, R2::RankOneMatrix)
    # middle product
    temp = fast_dot(R1.v, R2.u)
    return RankOneMatrix(R1.u * temp, R2.v)
end

Base.Matrix(R::RankOneMatrix) = R.u * R.v'
Base.collect(R::RankOneMatrix) = Matrix(R)
Base.copymutable(R::RankOneMatrix) = Matrix(R)
Base.copy(R::RankOneMatrix) = RankOneMatrix(copy(R.u), copy(R.v))

function Base.convert(::Type{<:RankOneMatrix{T,Vector{T},Vector{T}}}, R::RankOneMatrix) where {T}
    return RankOneMatrix(convert(Vector{T}, R.u), convert(Vector{T}, R.v))
end

function LinearAlgebra.dot(
    R::RankOneMatrix{T1},
    S::SparseArrays.AbstractSparseMatrixCSC{T2},
) where {T1<:Real,T2<:Real}
    (m, n) = size(R)
    T = promote_type(T1, T2)
    if (m, n) != size(S)
        throw(DimensionMismatch("Size mismatch"))
    end
    s = zero(T)
    if m * n == 0
        return s
    end
    rows = SparseArrays.rowvals(S)
    vals = SparseArrays.nonzeros(S)
    @inbounds for j in 1:n
        for ridx in SparseArrays.nzrange(S, j)
            i = rows[ridx]
            v = vals[ridx]
            s += v * R.u[i] * R.v[j]
        end
    end
    return s
end

Base.@propagate_inbounds function Base.:-(a::RankOneMatrix, b::RankOneMatrix)
    @boundscheck size(a) == size(b) || throw(DimensionMismatch())
    r = similar(a)
    @inbounds for j in 1:size(a, 2)
        for i in 1:size(a, 1)
            r[i, j] = a.u[i] * a.v[j] - b.u[i] * b.v[j]
        end
    end
    return r
end

Base.@propagate_inbounds function Base.:+(a::RankOneMatrix, b::RankOneMatrix)
    @boundscheck size(a) == size(b) || throw(DimensionMismatch())
    r = similar(a)
    @inbounds for j in 1:size(a, 2)
        for i in 1:size(a, 1)
            r[i, j] = a.u[i] * a.v[j] + b.u[i] * b.v[j]
        end
    end
    return r
end


text:
abstract type AbstractClock end

mutable struct Clock <: AbstractClock
    t0::Float64
end

mutable struct ClockTimer <: AbstractClock
    t0::Float64
    t1::Float64
    done::Bool
end

ClockTimer() = Clock(time(), Inf, false)
Clock() = Clock(time())

get_time(c::AbstractClock) = time() - c.t0
set_time!(c::AbstractClock, t0) = (c.t0 = t0)
reset!(c::AbstractClock) = set_time!(c, time())
set_endtime!(c::ClockTimer, t1) = (c.t1 = t1)

function get_time!(c::ClockTimer)
    t1 = time() - c.t0
    if t1 >= c.t1
        c.done = true
    end
    t1
end

function reset!(c::ClockTimer)
    c.t0 = time()
    c.done = false
end

function is_done!(c::ClockTimer)
    get_time!(c)
    c.done
end