text:
import Term: Panel, TextBox

@testset "\e[34mPanel - no content" begin
    for style in ("default", "red", "on_blue")
        testpanel(
            Panel(;fit=true, style=style), 3, 2
        )

        testpanel(
            Panel(), 88, 2
        )

        testpanel(
            Panel(; width=12, height=4, style=style), 12, 4
        )
    end

end

@testset "\e[34mPANEL - fit - measure" begin
    for style in ("default", "red", "on_blue")
        for justify in (:left, :center, :right)
            # ----------------------------- text only content ---------------------------- #
            testpanel(
                Panel("t"; fit=true, style=style), 7, 3
            )
            testpanel(
                Panel("test"; fit=true, style=style), 10, 3
            )
            testpanel(
                Panel("1234\n123456789012"; fit=true, style=style), 18, 4
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아"; fit=true, style=style), 28, 3
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아\n1234567890123456789012"; fit=true, style=style), 28, 4
            )
            testpanel(
                Panel("."^500; fit=true, style=style), displaysize(stdout)[2], nothing
            )

            # ------------------------------- nested panels ------------------------------ #
            testpanel(
                Panel(
                    Panel("test"; fit=true, style=style);
                fit=true, style=style),
                16, 5
            )

            testpanel(
                Panel(
                    Panel(Panel("."; fit=true, style=style); fit=true, style=style);
                fit=true, style=style), 19, 7
            )

            # @test_nothrow Panel(
            #         Panel("."^250; fit=true, style=style); fit=true, style=style
            #     )
        end
    end
end

@testset "\e[34mPANEL - nofit - measure" begin
    for style in ("default", "red", "on_blue")
        for justify in (:left, :center, :right)
            # ----------------------------- text only content ---------------------------- #
            testpanel(
                Panel("t"; style=style), 88, 3
            )
            testpanel(
                Panel("test"; style=style), 88, 3
            )
            testpanel(
                Panel("1234\n123456789012"; style=style), 88, 4
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아"; style=style), 88, 3
            )
            testpanel(
                Panel("나랏말싸미 듕귁에 달아\n1234567890123456789012"; style=style), 88, 4
            )
            testpanel(
                Panel("."^1500; style=style), 88, 21
            )

            # ------------------------------- nested panels ------------------------------ #
            testpanel(
                Panel(
                    Panel("test");
                fit=true, style=style), 94, 5
            )

            testpanel(
                Panel(
                    Panel(Panel("."); fit=true, style=style);
                fit=true, style=style), 100, 7
            )

            testpanel(
                Panel(
                    Panel("."^250); fit=true, style=style
                ), 94, 8
            )

            testpanel(
                Panel(
                    Panel("test"; style=style);
            ), 94, 5
            )

            testpanel(
                Panel(
                    Panel(Panel("."; style=style); style=style);
            ), 100, 7
            )

            testpanel(
                Panel(
                    Panel("."^250; style=style);
                ), 94, 8
            )

            testpanel(
                Panel(
                    Panel("t1"; style=style),
                    Panel("t2"; style=style),
                ), 94, 8
            )

            testpanel(
                Panel(
                    Panel("test", width=22);  width=30, height=8
                ), 30, 8
            )

            testpanel(
                Panel(
                    Panel("test", width=42);  width=30, height=8
                ), 48, 8
            )

            testpanel(
                Panel(
                    Panel("test", width=42,height=12);  width=30, height=8
                ), 48, 14
            )
        end
    end
end

@testset "PANEL - compare to string" begin
    pts = """
Lorem[red] ipsum dolor s[/red]it amet, consectetur adipiscing elit,
ed do eiusmod tempor [bold blue]incididu[underline]nt ut labore et dolore magna aliqua. Ut enim ad minim
veniam, quis nos[/underline]trud exercitation ullamco laboris nisi ut aliquip ex
ea commodo consequat. Duis aute [/bold blue][red on_black]irure dolor in reprehenderit
in voluptate velit esse cillum dolore eu fugiat nulla
pariatur. Excepteur sint occaecat[/red on_black] cupidatat non proident,
sunt in [green]culpa qui officia[/green] deserunt mollit anim
id est laborum."""

     pts = replace(pts, "\n"=>"")

    pts2 = replace("""
Lorem ipsum [bold]dolor sit[/bold] amet, consectetur adipiscing elit,
ed do e[red]iusmod tempor incididunt[/red] ut [bold]labore et [underline]dolore[/underline] magna aliqua.[/bold] Ut enim ad minim
veniam, quis[green] nostrud exercitation [on_black]ullamco laboris nisi ut aliquip ex [/on_black]
ea commodo consequat.[blue] Duis aute irure dolor in[/blue] reprehenderit
in voluptate velit[/green] esse [italic]cillum dolore[/italic][red] eu[/red][italic green] fugiat [/italic green]nulla
pariatur. Excepteur[red] sint[/red][blue] occaecat cupidatat [/blue]non proident,
sunt in culpa qui [italic]officia[/italic] deserunt mollit anim
id est laborum.""", "\n"=>"")

    if dotest
        @test string(Panel(pts; style="red", width=44)) == "\e[31m╭──────────────────────────────────────────╮\e[39m\n\e[31m│\e[39m  Lorem\e[31m ipsum dolor s\e[39mit amet, consectetu\e[0m  \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  r adipiscing elit,ed do eiusmod tempor\e[0m  \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34mincididu\e[4mnt ut labore et dolore magna\e[0m    \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[4maliqua. Ut enim ad minimveniam, quis\e[0m    \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[4mnos\e[24m\e[34mtrud exercitation ullamco laboris\e[0m    \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34mnisi ut aliquip ex ea commodo consequa\e[0m  \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34mt. Duis aute \e[31m\e[40mirure dolor in reprehende\e[0m  \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34m\e[31m\e[40mrit in voluptate velit esse cillum\e[0m      \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34m\e[31m\e[40mdolore eu fugiat nulla pariatur.\e[0m        \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34m\e[31m\e[40mExcepteur sint occaecat cupidatat\e[0m       \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34mnon proident, sunt in \e[32mculpa qui offic\e[0m   \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34m\e[32mia deserunt mollit anim id est laborum\e[0m  \e[0m\e[31m│\e[39m\e[0m\n\e[31m│\e[39m  \e[1m\e[34m\e[34m\e[49m.\e[0m                                       \e[0m\e[31m│\e[39m\e[0m\n\e[31m╰──────────────────────────────────────────╯\e[39m"

        @test string(Panel(pts; width=60)) == "\e[22m╭──────────────────────────────────────────────────────────╮\e[22m\n\e[22m│\e[22m  Lorem\e[31m ipsum dolor s\e[39mit amet, consectetur adipiscing\e[0m      \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  elit,ed do eiusmod tempor \e[1m\e[34mincididu\e[4mnt ut labore et\e[0m       \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[4mdolore magna aliqua. Ut enim ad minimveniam, quis nos\e[24m\e[34m\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[34mtrud exercitation ullamco laboris nisi ut aliquip ex\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[34mea commodo consequat. Duis aute \e[31m\e[40mirure dolor in repreh\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[34m\e[31m\e[40menderit in voluptate velit esse cillum dolore eu fugia\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[34m\e[31m\e[40mt nulla pariatur. Excepteur sint occaecat cupidatat\e[0m     \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[34mnon proident, sunt in \e[32mculpa qui officia deserunt\e[0m        \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[1m\e[34m\e[34m\e[49mmollit anim id est laborum.\e[0m                             \e[0m\e[22m│\e[22m\e[0m\n\e[22m╰──────────────────────────────────────────────────────────╯\e[22m"
    end

    @test string(Panel(pts2; width=44)) == "\e[22m╭──────────────────────────────────────────╮\e[22m\n\e[22m│\e[22m  Lorem ipsum \e[1mdolor sit\e[22m amet, consectetu\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  r adipiscing elit,ed do e\e[31miusmod tempor\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[31mincididunt\e[39m\e[22m ut \e[1mlabore et \e[4mdolore\e[24m\e[1m magna\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[1maliqua.\e[22m Ut enim ad minimveniam, quis\e[32m\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32mnostrud exercitation \e[40mullamco laboris\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[40mnisi ut aliquip ex \e[49m\e[32mea commodo consequa\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32mt.\e[34m Duis aute irure dolor in\e[39m\e[32m reprehende\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[32mrit in voluptate velit\e[39m\e[22m esse \e[3mcillum\e[0m      \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3mdolore\e[23m\e[22m\e[31m eu\e[39m\e[22m\e[3m\e[32m fugiat \e[22mnulla pariatur.\e[0m        \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32m\e[22mExcepteur\e[31m sint\e[39m\e[22m\e[34m occaecat cupidatat \e[39m\e[22mnon\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32m\e[22mproident, sunt in culpa qui \e[3mofficia\e[23m\e[22m\e[0m     \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32mdeserunt mollit anim id est laborum.\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m╰──────────────────────────────────────────╯\e[22m"

    @test string(Panel(pts2; width=51)) == "\e[22m╭─────────────────────────────────────────────────╮\e[22m\n\e[22m│\e[22m  Lorem ipsum \e[1mdolor sit\e[22m amet, consectetur adipi\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  scing elit,ed do e\e[31miusmod tempor incididunt\e[39m\e[22m\e[0m     \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22mut \e[1mlabore et \e[4mdolore\e[24m\e[1m magna aliqua.\e[22m Ut enim\e[0m      \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1mad minimveniam, quis\e[32m nostrud exercitation \e[40mul\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[40mlamco laboris nisi ut aliquip ex \e[49m\e[32mea commodo\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32mconsequat.\e[34m Duis aute irure dolor in\e[39m\e[32m reprehen\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[32mderit in voluptate velit\e[39m\e[22m esse \e[3mcillum dolore\e[23m\e[22m\e[31m\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[31meu\e[39m\e[22m\e[3m\e[32m fugiat \e[22mnulla pariatur. Excepteur\e[31m sint\e[39m\e[22m\e[34m\e[0m       \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32m\e[34moccaecat cupidatat \e[39m\e[22mnon proident, sunt in\e[0m       \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32m\e[22mculpa qui \e[3mofficia\e[23m\e[22m deserunt mollit anim id est\e[0m  \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32mlaborum.\e[0m                                       \e[0m\e[22m│\e[22m\e[0m\n\e[22m╰─────────────────────────────────────────────────╯\e[22m"

    @test string(Panel(pts2; width=67)) == "\e[22m╭─────────────────────────────────────────────────────────────────╮\e[22m\n\e[22m│\e[22m  Lorem ipsum \e[1mdolor sit\e[22m amet, consectetur adipiscing elit,ed\e[0m     \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  do e\e[31miusmod tempor incididunt\e[39m\e[22m ut \e[1mlabore et \e[4mdolore\e[24m\e[1m magna aliqu\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[1ma.\e[22m Ut enim ad minimveniam, quis\e[32m nostrud exercitation \e[40mullamco\e[0m   \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[40mlaboris nisi ut aliquip ex \e[49m\e[32mea commodo consequat.\e[34m Duis aute\e[0m     \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[34mirure dolor in\e[39m\e[32m reprehenderit in voluptate velit\e[39m\e[22m esse \e[3mcillum\e[0m    \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3mdolore\e[23m\e[22m\e[31m eu\e[39m\e[22m\e[3m\e[32m fugiat \e[22mnulla pariatur. Excepteur\e[31m sint\e[39m\e[22m\e[34m occaecat\e[0m       \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32m\e[34mcupidatat \e[39m\e[22mnon proident, sunt in culpa qui \e[3mofficia\e[23m\e[22m deserunt\e[0m     \e[0m\e[22m│\e[22m\e[0m\n\e[22m│\e[22m  \e[22m\e[1m\e[32m\e[32m\e[22m\e[3m\e[32mmollit anim id est laborum.\e[0m                                    \e[0m\e[22m│\e[22m\e[0m\n\e[22m╰─────────────────────────────────────────────────────────────────╯\e[22m"
end

@testset "\e[34mPanel + renderables" begin
    testpanel(
        Panel(
            RenderableText("x".^5)
        ), 88, 3
    )

    testpanel(
        Panel(
            RenderableText("x".^500)
        ), 88, nothing
    )

    testpanel(
        Panel(
            RenderableText("x".^5); fit=true
        ), 11, 3
    )

    testpanel(
        Panel(
            RenderableText("x".^500); fit=true
        ), displaysize(stdout)[2], nothing
    )

end

@testset "\e[34mPANEL - titles" begin
    for fit in (true, false)
        for justify in (:left, :center, :right)
            for style in ("red", "bold", "on_green")
                testpanel(
                    Panel("."^50, title="test",
                            title_style=style,
                            title_justify=justify,
                            subtitle="subtest",
                            subtitle_style=style,
                            subtitle_justify=justify,
                            fit=fit
                            ),
                    fit ? nothing : 88,
                    nothing
                )

                testpanel(
                    Panel(
                        Panel("."^50, title="test",
                                title_style=style,
                                title_justify=justify,
                                subtitle="subtest",
                                subtitle_style=style,
                                subtitle_justify=justify,
                                  fit=fit,
                                )
                    ),
                    fit ? nothing : 94,
                    nothing
                )
            end
        end
    end
end

@testset "\e[34mTBOX" begin
    w = displaysize(stdout)[2]

    for justify in (:left, :center, :right)

        testpanel(
            TextBox(
                "nofit"^25;
                width=1000,
                justify=justify
            ), w - 4, nothing)

        testpanel(
            TextBox(
                "truncate"^25;
                width=100,
                fit=:truncate,
                justify=justify
        ), 100, 3)

        testpanel(
            TextBox(
                "truncate"^25;
                width=100,
                justify=justify
            ), 100, 5)

        testpanel(
            TextBox(
                "truncate"^8;
                fit=:fit,
                justify=justify
            ), 68, 4)

        # testpanel(
        #     TextBox(
        #         "[red]truncate[/red]"^8;
        #         fit=:fit,
        #         justify=justify
        #     ), 68, 4)

        # testpanel(
        #     TextBox(
        #         "[red]truncate[/red]test"^8;
        #         fit=:fit,
        #         justify=justify
        #     ), 100, 4)

        # testpanel(TextBox(
        #         "[red]tru\nncate[/red]test"^1;
        #         fit=:fit,
        #         justify=justify
        #     ), 13, 6)
    end
end

@testset "\e[34mPanel - padding" begin

    p = Panel("."^24; padding = (4, 4, 2, 2))
    testpanel(p, 88, 7)
    # @test string(p) == "\e[22m╭──────────────────────────────────────────────────────────────────────────────────────╮\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m│\e[22m    ........................                                                          \e[22m│\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m│\e[22m                                                                                \e[22m│\e[22m\n\e[22m╰──────────────────────────────────────────────────────────────────────────────────────╯\e[22m"

    p = Panel("."^24; padding = (4, 4, 2, 2), fit=true)
    testpanel(p, 34, 7)

end


text:
using UVH5Splicer
using Test

@testset "UVH5Splicer.jl" begin
    # Write your tests here.
end

ext: jl


text:
using OrdinaryDiffEq
using Trixi

###############################################################################
# semidiscretization of the (inviscid) Burgers' equation

equations = InviscidBurgersEquation1D()

function initial_condition_linear_stability(x, t, equation::InviscidBurgersEquation1D)
  k = 1
  2 + sinpi(k * (x[1] - 0.7)) |> SVector
end

polydeg = 3
surface_flux = flux_ec
volume_flux = flux_ec
solver = DGSEM(polydeg, surface_flux, VolumeIntegralFluxDifferencing(volume_flux))

coordinates_min = -1.0
coordinates_max =  1.0
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=4,
                n_cells_max=10_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_linear_stability, solver)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval,
                                     extra_analysis_errors=(:l2_error_primitive,
                                                            :linf_error_primitive))

alive_callback = AliveCallback(analysis_interval=analysis_interval)

stepsize_callback = StepsizeCallback(cfl=0.8)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary


text:
module GetSpectrum

using DSP
using FFTW

export getSpectrum
export getWaterfall
#greet() = print("Hello World!")

function getSpectrum(fs,sig;N=nothing)
	if isnothing(N)
		N = length(sig);
	end
    freqAx = collect(((0:N-1)./N .- 0.5)*fs);
	ss	   = @view sig[1:N];
	y	   = 10*log10.(abs2.(fftshift(fft(ss))));
	return (freqAx,y);
end
getSpectrum(sig) = getSpectrum(1,sig);

function getWaterfall(fe,sig;sizeFFT=1024)
    nbSeg   = length(sig) ÷ sizeFFT;
    ss      =  @views sig[1:nbSeg*sizeFFT];
    ss      = reshape(ss,sizeFFT,nbSeg)
    sMatrix = zeros(Float64,sizeFFT,nbSeg);
    for iN = 1 : 1 : nbSeg
        sMatrix[:,iN] = abs2.(fftshift(fft(ss[:,iN])));
    end
    fAx = collect(((0:1:sizeFFT-1)./sizeFFT .- 0.5) .* fe);
    tAx = (0:nbSeg-1) * (sizeFFT/fe);
    return tAx,fAx,sMatrix;
end
getWaterfall(sig;sizeFFT=1024) = getWaterfall(1,sig;sizeFFT=sizeFFT);
end # module


text:
#!/usr/bin/env julia --project
using Test

include("hydrostatic_spindown.jl")
ClimateMachine.init()

const FT = Float64

#################
# RUN THE TESTS #
#################
@testset "$(@__FILE__)" begin

    include("../refvals/hydrostatic_spindown_refvals.jl")

    # simulation time
    timeend = FT(15 * 24 * 3600) # s
    tout = FT(24 * 3600) # s
    timespan = (tout, timeend)

    # DG polynomial order
    N = Int(4)

    # Domain resolution
    Nˣ = Int(5)
    Nʸ = Int(5)
    Nᶻ = Int(8)
    resolution = (N, Nˣ, Nʸ, Nᶻ)

    # Domain size
    Lˣ = 1e6  # m
    Lʸ = 1e6  # m
    H = 400  # m
    dimensions = (Lˣ, Lʸ, H)

    BC = (
        OceanBC(Impenetrable(FreeSlip()), Insulating()),
        OceanBC(Penetrable(FreeSlip()), Insulating()),
    )
    config = SplitConfig(
        "rotating_bla",
        resolution,
        dimensions,
        Coupled(),
        Rotating();
        solver = SplitExplicitSolver,
        boundary_conditions = BC,
    )

    #=
    BC = (
        ClimateMachine.Ocean.SplitExplicit01.OceanFloorFreeSlip(),
        ClimateMachine.Ocean.SplitExplicit01.OceanSurfaceNoStressNoForcing(),
    )

    config = SplitConfig(
        "rotating_jmc",
        resolution,
        dimensions,
        Coupled(),
        Rotating();
        solver = SplitExplicitLSRK2nSolver,
        boundary_conditions = BC,
    )
    =#

    run_split_explicit(
        config,
        timespan,
        dt_fast = 300,
        dt_slow = 300, # 90 * 60,
        # refDat = refVals.ninety_minutes,
        analytic_solution = true,
    )
end


text:
export Dataset, DatasetVarInfo, DatasetParseParams, DatasetWriteParams
export DoubleArray
export Header
export IdArray
export Network
export Pattern
export Node
export NodeDefinition
export NodeValue
export SysCoordinates
export IntArray

type Dataset
	ptr::Ptr{Void}

	function Dataset()
		ptr = ccall( (:createDataset, LIB_SMILE), Ptr{Void}, ())
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeDataset, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
	function Dataset( dset::Dataset )
		ptr = ccall( (:copyDataset, LIB_SMILE), Ptr{Void}, (Ptr{Void},), dset.ptr)
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeDataset, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
end

immutable DMatrix
	ptr::Ptr{Void}
end

immutable Header
	ptr::Ptr{Void}
end

for str = ("DatasetParseParams", "DatasetWriteParams", "DoubleArray", "IdArray","Network","Pattern",)
	fname = symbol(str)
	f_sym

	op_cr = QuoteNode(symbol("create$str"))
	op_fr = QuoteNode(symbol("free$str"))

	@eval begin
		type $fname
			ptr::Ptr{Void}

			($fname)(ptr::Ptr{Void}) = new(ptr)
			function ($fname)()
				ptr = ccall( ($op_cr, LIB_SMILE), Ptr{Void}, ())
				smart_p = new(ptr)
				finalizer(smart_p, obj -> ccall( ($op_fr, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
				smart_p
			end
		end
	end
end

immutable Node
	ptr::Ptr{Void}
end

immutable NodeDefinition
	ptr::Ptr{Void}
end

immutable NodeValue
	ptr::Ptr{Void}
end

type SysCoordinates
	ptr::Ptr{Void}

	SysCoordinates( ptr::Ptr{Void} ) = new(ptr)
	function SysCoordinates( nodedef::NodeDefinition )
		ptr = ccall( (:createSysCoordinatesFromNodeDefinition, LIB_SMILE), Ptr{Void}, (Ptr{Void},), nodedef.ptr )
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeSysCoordinates, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
	function SysCoordinates( nodeval::NodeValue )
		ptr = ccall( (:createSysCoordinatesFromNodeValue, LIB_SMILE), Ptr{Void}, (Ptr{Void},), nodeval.ptr )
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeSysCoordinates, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
end

type IntArray
	ptr::Ptr{Void}

	function IntArray( ptr::Ptr{Void} )
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeIntArray, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
	function IntArray()
		ptr = ccall( (:createIntArray, LIB_SMILE), Ptr{Void}, () )
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeIntArray, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
	function IntArray(initialSize::Integer)
		ptr = ccall( (:createIntArray_InitialSize, LIB_SMILE), Ptr{Void}, (Cint,), initialSize )
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeIntArray, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
	function IntArray(likeThisOne::IntArray)
		ptr = ccall( (:createIntArray_Copy, LIB_SMILE), Ptr{Void}, (Ptr{Void},), likeThisOne.ptr )
		smart_p = new(ptr)
		finalizer(smart_p, obj -> ccall( (:freeIntArray, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))
		smart_p
	end
	function IntArray{I<:Integer}(arr::Vector{I})
	    n = length(arr)

	    ptr = ccall( (:createIntArray_InitialSize, LIB_SMILE), Ptr{Void}, (Cint,), n )
	    smart_p = new(ptr)
	    finalizer(smart_p, obj -> ccall( (:freeIntArray, LIB_SMILE), Void, (Ptr{Void},), obj.ptr ))

	    for i = 1 : n
	        ccall( (:intarray_SetIndex, LIB_SMILE), Void, (Ptr{Void},Int32,Int32), ptr, i-1, arr[i] )
	    end

	    smart_p
	end
end


text:
const HubbardConf = Array{Int8, 2} # conf === hsfield === discrete Hubbard Stratonovich field (Hirsch field)
const HubbardDistribution = Int8[-1,1]

"""
Famous attractive (negative U) Hubbard model on a cubic lattice.
Discrete Hubbard Stratonovich transformation (Hirsch transformation) in the density/charge channel,
such that HS-field is real.

    HubbardModelAttractive(; dims, L[, kwargs...])

Create an attractive Hubbard model on `dims`-dimensional cubic lattice
with linear system size `L`. Additional allowed `kwargs` are:

 * `mu::Float64=0.0`: chemical potential
 * `U::Float64=1.0`: onsite interaction strength, "Hubbard U"
 * `t::Float64=1.0`: hopping energy
"""
@with_kw_noshow struct HubbardModelAttractive{C<:AbstractCubicLattice} <: Model
    # user mandatory
    dims::Int
    L::Int

    # user optional
    mu::Float64 = 0.0
    U::Float64 = 1.0
    @assert U >= 0. "U must be positive."
    t::Float64 = 1.0

    # non-user fields
    l::C = choose_lattice(HubbardModelAttractive, dims, L)
    neighs::Matrix{Int} = neighbors_lookup_table(l)
    flv::Int = 1
end

function choose_lattice(::Type{HubbardModelAttractive}, dims::Int, L::Int)
    if dims == 1
        return Chain(L)
    elseif dims == 2
        return SquareLattice(L)
    else
        return CubicLattice(dims, L)
    end
end

"""
    HubbardModelAttractive(params::Dict)
    HubbardModelAttractive(params::NamedTuple)

Create an attractive Hubbard model with (keyword) parameters as specified in the
dictionary/named tuple `params`.
"""
HubbardModelAttractive(params::Dict{Symbol, T}) where T = HubbardModelAttractive(; params...)
HubbardModelAttractive(params::NamedTuple) = HubbardModelAttractive(; params...)

# cosmetics
import Base.summary
import Base.show
Base.summary(model::HubbardModelAttractive) = "$(model.dims)D attractive Hubbard model"
Base.show(io::IO, model::HubbardModelAttractive) = print(io, "$(model.dims)D attractive Hubbard model, L=$(model.L) ($(model.l.sites) sites)")
Base.show(io::IO, m::MIME"text/plain", model::HubbardModelAttractive) = print(io, model)

# implement `Model` interface
@inline nsites(m::HubbardModelAttractive) = nsites(m.l)

# implement `DQMC` interface: mandatory
@inline Base.rand(::Type{DQMC}, m::HubbardModelAttractive, nslices::Int) = rand(HubbardDistribution, nsites(m), nslices)

"""
Calculates the hopping matrix \$T_{i, j}\$ where \$i, j\$ are
site indices.

Note that since we have a time reversal symmetry relating spin-up
to spin-down we only consider one spin sector (one flavor) for the attractive
Hubbard model in the DQMC simulation.

This isn't a performance critical method as it is only used once before the
actual simulation.
"""
function hopping_matrix(mc::DQMC, m::HubbardModelAttractive)
    N = nsites(m)
    neighs = m.neighs # row = up, right, down, left; col = siteidx

    T = diagm(0 => fill(-m.mu, N))

    # Nearest neighbor hoppings
    @inbounds @views begin
      for src in 1:N
        for nb in 1:size(neighs,1)
          trg = neighs[nb,src]
          T[trg,src] += -m.t
        end
      end
    end

    return T
end

"""
Calculate the interaction matrix exponential `expV = exp(- power * delta_tau * V(slice))`
and store it in `result::Matrix`.

This is a performance critical method.
"""
@inline function interaction_matrix_exp!(mc::DQMC, m::HubbardModelAttractive,
            result::Matrix, conf::HubbardConf, slice::Int, power::Float64=1.)
    dtau = mc.p.delta_tau
    lambda = acosh(exp(m.U * dtau/2))
    result .= spdiagm(0 => exp.(sign(power) * lambda * conf[:,slice]))
    nothing
end

@inline function propose_local(mc::DQMC, m::HubbardModelAttractive, i::Int, slice::Int, conf::HubbardConf)
    # see for example dos Santos (2002)
    greens = mc.s.greens
    dtau = mc.p.delta_tau
    lambda = acosh(exp(m.U * dtau/2))

    ΔE_boson = -2. * lambda * conf[i, slice]
    γ = exp(ΔE_boson) - 1
    detratio = (1 + γ * (1 - greens[i,i]))^2 # squared because of two spin sectors.

    return detratio, ΔE_boson, γ
end

@inline function accept_local!(mc::DQMC, m::HubbardModelAttractive, i::Int, slice::Int, conf::HubbardConf, delta, detratio, ΔE_boson::Float64)
    greens = mc.s.greens
    γ = delta

    u = -greens[:, i]
    u[i] += 1.
    # TODO: OPT: speed check, maybe @views/@inbounds
    greens .-= kron(u * 1. /(1 + γ * u[i]), transpose(γ * greens[i, :]))
    conf[i, slice] *= -1
    nothing
end

# implement DQMC interface: optional
"""
Green's function is real for the attractive Hubbard model.
"""
@inline greenseltype(::Type{DQMC}, m::HubbardModelAttractive) = Float64

"""
Calculate energy contribution of the boson, i.e. Hubbard-Stratonovich/Hirsch field.
"""
@inline function energy_boson(m::HubbardModelAttractive, hsfield::HubbardConf)
  dtau = mc.p.delta_tau
    lambda = acosh(exp(m.U * dtau/2))
    return lambda * sum(hsfield)
end

include("observables.jl")

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
using Random,
    DataFrames, XLSX, LinearAlgebra, Statistics, StatsBase, Distributions, Plots, CSV, RCall, JLD;
include("fcns/mainsup.jl")

function phidraw(y,x,m,p,nsmpl,t)
    maxλ = 1
    ϕd = nothing
    while maxλ >= 1
        pos = rand(1:nsmpl-t)
        yn = y[pos:pos+t-1,:]
        xn = x[pos:pos+t-1,:]
        βd, ϕd = OLSbetas(yn, xn, m, p, Xblock = true, Nx = 3)
        maxλ = maximum(abs.(eigvals(ϕd)))
    end
    return ϕd
end

function boots(B,Y,E,p,nsmpl)
    t, m = size(Y)
    nx  = m*p+1
    ebo = E[rand(1:t,nsmpl),:]
    rst = rand(p+1:t)
    y = Array{Float32,2}(undef,m,nsmpl)
    x = Array{Float32,2}(undef,nx,nsmpl)
    x[:,1] = [vec(Y[rst:-1:rst-1,:]') ; 1]
    βm = B'
    for i in 1:nsmpl-1
        y[:,i] = βm * x[:,i] + ebo[i,:]
        x[:,i+1] = [y[:,i] ; x[1:m,i] ; 1]
    end
    y[:,nsmpl] = βm * x[:,nsmpl] + ebo[nsmpl,:]
    y = y'
    x = x'
    collΦ = Array{Any,1}(undef,nrep)
    for i in 1:nrep
        collΦ[i] = phidraw(y,x,m,p,nsmpl,t)
    end
    return collΦ
end

p = 2
h = 40
nrep = 1000
nsmpl = 100000
dd = load("GSdata.jld")
dataset = dd["dataset"]
nc = nfields(dataset);
Φdraw = Array{Any,1}(undef,nc)
for k in 1:nc
    Y, X = VARData(dataset[k],p)
    B = (X'*X) \ (X'*Y)
    E = Y - X*B
    Φdraw[k] = boots(B,Y,E,p,nsmpl)
end

max_stars_repo_name: ChrisVeJa/GlobalShock

max_issues_repo_name: ChrisVeJa/GlobalShock

max_forks_repo_name: ChrisVeJa/GlobalShock


text:
# Load the packages
using GraphNeuralNetworks, DiffEqFlux, DifferentialEquations
using Flux: onehotbatch, onecold
using Flux.Losses: logitcrossentropy
using Statistics: mean
using MLDatasets: Cora
using CUDA
# CUDA.allowscalar(false) # Some scalar indexing is still done by DiffEqFlux

# device = cpu # `gpu` not working yet
device = CUDA.functional() ? gpu : cpu

# LOAD DATA
data = Cora.dataset()
g = GNNGraph(data.adjacency_list) |> device
X = data.node_features |> device
y = onehotbatch(data.node_labels, 1:data.num_classes) |> device
train_ids = data.train_indices |> device
val_ids = data.val_indices |> device
test_ids = data.test_indices |> device
ytrain = y[:, train_ids]

# Model and Data Configuration
nin = size(X, 1)
nhidden = 16
nout = data.num_classes
epochs = 40

# Define the Neural GDE
diffeqsol_to_array(x) = reshape(device(x), size(x)[1:2])

# GCNConv(nhidden => nhidden, graph=g),

node_chain = GNNChain(GCNConv(nhidden => nhidden, relu),
                      GCNConv(nhidden => nhidden, relu)) |> device

node = NeuralODE(WithGraph(node_chain, g),
                (0.f0, 1.f0), Tsit5(), save_everystep = false,
                reltol = 1e-3, abstol = 1e-3, save_start = false) |> device

model = GNNChain(GCNConv(nin => nhidden, relu),
                 Dropout(0.5),
                 node,
                 diffeqsol_to_array,
                 Dense(nhidden, nout)) |> device

# Loss
loss(x, y) = logitcrossentropy(model(g, x), y)
accuracy(x, y) = mean(onecold(model(g, x)) .== onecold(y))

# # Training
# ## Model Parameters
ps = Flux.params(model);

# ## Optimizer
opt = ADAM(0.01)

# ## Training Loop
for epoch in 1:epochs
    gs = gradient(() -> loss(X, y), ps)
    Flux.Optimise.update!(opt, ps, gs)
    @show(accuracy(X, y))
end


text:
#Tests the jacobian
network1 = @reaction_network begin
    (2.0,1.0),   ∅ ↔ X
    (3.0,1.0),   ∅ ↔ Y
    (5.0,2.0),   X + Y ↔ XY
end
J = zeros(3,3)
network1.jac(J,[1.,1.,1.],[],0.)
@test J == [-6. -5. 2.; -5. -6. 2.; 5. 5. -2.]

network2 = @min_reaction_network begin
    (2.0,1.0),   ∅ ↔ X
    (3.0,1.0),   ∅ ↔ Y
    (5.0,2.0),   X + Y ↔ XY
end
addodes!(network2)
J = zeros(3,3)
network2.jac(J,[1.,1.,1.],[],0.)
@test J == [-6. -5. 2.; -5. -6. 2.; 5. 5. -2.]

network3 = @reaction_network begin
    (p1,1.0),   ∅ ↔ X
    (p2,1.0),   ∅ ↔ Y
    (p3*X,1.0), X + Y ↔ XY
end p1 p2 p3
J = zeros(3,3)
for i = 1:100
    u = 10*rand(3)
    p = 10*rand(3)
    network3.jac(J,u,p,0.)
    @test maximum(J .- [-1-2*p[3]*u[2]*u[1] -p[3]*u[1]*u[1] 1.; -2*p[3]*u[2]*u[1] -1-p[3]*u[1]*u[1] 1; 2*p[3]*u[2]*u[1] p[3]*u[1]*u[1] -1.])<0.0001
end

network4 = @reaction_network begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
    k6, 0 --> 2B
    hill(A,k7,k8,2), ∅ --> B
end k1 k2 k3 k4 k5 k6 k7 k8
p = rand(8)
u = rand(3)
t = 0.
function n4f(du, u, p, t)
    A = u[1]; B = u[2]; C = u[3];
    k1 = p[1]; k2 = p[2]; k3 = p[3]; k4 = p[4]; k5 = p[5]; k6 = p[6]; k7 = p[7]; k8 = p[8]
    du[1] = -k1*A^2 + 2*k2*B - k3*A*B + k4*C + k5*C^3/2
    du[2] = .5*k1*A^2 -k2*B -k3*A*B + k4*C + 2*k6 + k7*A^2/(k8^2+A^2)
    du[3] = k3*A*B-k4*C-k5*C^3/2
end
function n4Jac(J, u, p, t)
    A = u[1]; B = u[2]; C = u[3];
    k1 = p[1]; k2 = p[2]; k3 = p[3]; k4 = p[4]; k5 = p[5]; k6 = p[6]; k7 = p[7]; k8 = p[8]
    J[1,1] = -2*k1*A - k3*B
    J[1,2] = 2*k2 - k3*A
    J[1,3] = k4 + 3*k5*C^2 / 2
    J[2,1] = k1*A -k3*B + 2*k7*k8^2*A/(k8^2+A^2)^2
    J[2,2] = -k2-k3*A
    J[2,3] = k4
    J[3,1] = k3*B
    J[3,2] = k3*A
    J[3,3] = -k4 - 3*k5*C^2/2
end
du  = zeros(3); network4.f(du,u,p,t);
dua = zeros(3); n4f(dua,u,p,t);
J   = zeros(3,3); network4.jac(J,u,p,t);
Ja  = zeros(3,3); n4Jac(Ja,u,p,t);
@test maximum(abs.(du .- dua)) < 10*eps(Float64)
@test maximum(abs.(J .- Ja)) < 10*eps(Float64)

# Test the sparse jacobian
network5 = @min_reaction_network begin
    k1, 2A --> B
    k2, B --> 2A
    k3, A + B --> C
    k4, C --> A + B
    k5, 3C --> 3A
    k6, 0 --> 2B
    hill(A,k7,k8,2), ∅ --> B
end k1 k2 k3 k4 k5 k6 k7 k8
addodes!(network5, sparse_jac=true)
Js = sparse(ones(3,3))
network5.jac(Js, u, p, t)
@test maximum(abs.(Matrix(Js) .- Ja)) < 10*eps(Float64)

#Tests the parameter jacobian.
network1 = @reaction_network begin
    (p1,1.0),   ∅ ↔ X
    (1.0,p2),   X ↔ Y
    (p3,p4),    Y ↔ Z
end p1 p2 p3 p4
J = zeros(3,4)
network1.paramjac(J,[2.,3.,4.],[2.,3.,4.,5.],0.)
@test J == [1. 3. 0. 0.; -0. -3. -3. 4.; 0. 0. 3. -4.]

network2 = @reaction_network begin
    (1.,p1*p2),   ∅ ↔ X
    (p3^3+p1,p2),   2X ↔ Y
end p1 p2 p3
J = zeros(2,3)
for i = 1:100
    u = 10*rand(2)
    p = 10*rand(3)
    network2.paramjac(J,u,p,0.)
    @test maximum(J .- [-u[1]*p[2]-u[1]^2 -u[1]*p[1]+2*u[2] -6*p[3]^2*u[1]^2/2; u[1]^2/2 -u[2] 3*p[3]^2*u[1]^2/2])<0.0001
end


text:
using DiscreteValueIteration

pomdp = TigerPOMDP()

mdp = UnderlyingMDP(pomdp)

@test n_states(mdp) == n_states(pomdp)
@test states(mdp) == states(pomdp)
s_mdp = rand(MersenneTwister(1), initial_state_distribution(mdp))
s_pomdp = rand(MersenneTwister(1), initial_state_distribution(pomdp))

@test s_mdp == s_pomdp

solver = ValueIterationSolver(max_iterations = 100)
mdp_policy = solve(solver, mdp, verbose=false)
pomdp_policy = solve(solver, pomdp, verbose=false)
@test mdp_policy.util == pomdp_policy.util

action_index(mdp, 1)


text:
abstract type StaticLoad <: ElectricLoad end

struct PowerLoad <: StaticLoad
    name::String
    available::Bool
    bus::Bus
    maxactivepower::Float64 # [MW]
    maxreactivepower::Float64 # [MVAr]
    scalingfactor::TimeSeries.TimeArray
end

function PowerLoad(name::String, available::Bool, bus::Bus, maxactivepower::Float64, maxreactivepower::Float64, scalingfactor=nothing)
    scalingfactor=TimeSeries.TimeArray(Dates.today(),ones(1))
    return PowerLoad(name, available, bus, maxactivepower, maxreactivepower, scalingfactor)
end

function PowerLoadPF(name::String, available::Bool, bus::Bus, maxactivepower::Float64, power_factor::Float64, scalingfactor::TimeSeries.TimeArray)
    maxreactivepower = maxactivepower*sin(acos(power_factor))
    return PowerLoad(name, available, bus, maxactivepower, maxreactivepower, scalingfactor)
end

function PowerLoadPF(name::String, available::Bool, bus::Bus, maxactivepower::Float64, power_factor::Float64, scalingfactor=nothing)
    scalingfactor=TimeSeries.TimeArray(Dates.today(),ones(1))
    maxreactivepower = maxactivepower*sin(acos(power_factor))
    return PowerLoad(name, available, bus, maxactivepower, maxreactivepower, scalingfactor)
end

PowerLoadPF(; name = "init", available = true, bus = Bus(), maxactivepower = 0.0, power_factor=1.0,
            scalingfactor=TimeSeries.TimeArray(Dates.today(),ones(1))) = PowerLoadPF(name, available, bus, maxactivepower, power_factor, scalingfactor)

PowerLoad(; name = "init", available = true, bus = Bus(), maxactivepower = 0.0, maxreactivepower=0.0,
            scalingfactor=TimeSeries.TimeArray(Dates.today(),ones(1))) = PowerLoad(name, available, bus, maxactivepower, maxreactivepower, scalingfactor)


text:
using LinearAlgebra

@with_kw mutable struct AffineParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    A::Matrix{T}
    b::Vector{T}
end

struct Affine <: AbstractProx
    params::AffineParameters{Float64}
    H::Cholesky{Float64,Matrix{Float64}}

    function (::Type{Affine})(; kwargs...)
        params = AffineParameters(; kwargs...)
        @unpack A = params
        return new(params,cholfact(A*A'))
    end
end

function prox!(affine::Affine,step::Real,xprev::AbstractVector,gcurr::AbstractVector,xcurr::AbstractVector)
    @unpack A,b = affine.params
    x = xprev - step*gcurr
    xcurr .= x + A'*(affine.H\(b-A*x))
end


text:
#=
    Helper functions to store constraint mappings
=#
cmap(model::LinQuadOptimizer) = model.constraint_mapping
# dummy method for isvalid
constrdict(::LinQuadOptimizer, ::CI{F, S}) where {F, S} = Dict{CI{F, S}, Any}()

function Base.getindex(model::LinQuadOptimizer, index::CI)
    return constrdict(model, index)[index]
end

function delete_constraint_name(model::LinQuadOptimizer, index)
    if haskey(model.constraint_names, index)
        name = model.constraint_names[index]
        delete!(model.constraint_names_rev, name)
        delete!(model.constraint_names, index)
    end
end

"""
    has_integer(model::LinQuadOptimizer)::Bool

A helper function to determine if `model` has any integer components (i.e.
binary, integer, special ordered sets, semicontinuous, or semi-integer
variables).
"""
function has_integer(model::LinQuadOptimizer)
    constraint_map = cmap(model)
    return length(constraint_map.integer) > 0 ||
        length(constraint_map.binary) > 0 ||
        length(constraint_map.sos1) > 0 ||
        length(constraint_map.sos2) > 0 ||
        length(constraint_map.semicontinuous) > 0 ||
        length(constraint_map.semiinteger) > 0
end

#=
    MOI.isvalid
=#

function MOI.isvalid(model::LinQuadOptimizer, index::CI{F,S}) where F where S
    dict = constrdict(model, index)
    return haskey(dict, index)
end

"""
    __assert_valid__(model::LinQuadOptimizer, index::MOI.Index)

Throw an MOI.InvalidIndex error if `MOI.isvalid(model, index) == false`.
"""
function __assert_valid__(model::LinQuadOptimizer, index::MOI.Index)
    if !MOI.isvalid(model, index)
        throw(MOI.InvalidIndex(index))
    end
end

"""
    __assert_supported_constraint__(model::LinQuadOptimizer, ::Type{F}, ::Type{S})

Throw an `UnsupportedConstraint{F, S}` error if the model cannot add constraints
of type `F`-in-`S`.
"""
function __assert_supported_constraint__(model::LinQuadOptimizer, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet}
    if !((F,S) in supported_constraints(model))
        throw(MOI.UnsupportedConstraint{F, S}())
    end
end

#=
    Get number of constraints
=#
function MOI.canget(model::LinQuadOptimizer, ::MOI.NumberOfConstraints{F, S}) where F where S
    return (F,S) in supported_constraints(model)
end

function MOI.get(model::LinQuadOptimizer, attribute::MOI.NumberOfConstraints{F, S}) where F where S
    __assert_supported_constraint__(model, F, S)
    return length(constrdict(model, CI{F,S}(UInt(0))))
end

#=
    Get list of constraint references
=#
function MOI.canget(model::LinQuadOptimizer, ::MOI.ListOfConstraintIndices{F, S}) where F where S
    return (F,S) in supported_constraints(model)
end

function MOI.get(model::LinQuadOptimizer, attribute::MOI.ListOfConstraintIndices{F, S}) where F where S
    __assert_supported_constraint__(model, F, S)
    dict = constrdict(model, CI{F,S}(UInt(0)))
    indices = collect(keys(dict))
    return sort(indices, by=x->x.value)
end

#=
    Get list of constraint types in model
=#
function MOI.canget(::LinQuadOptimizer, ::MOI.ListOfConstraints)
    return true
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ListOfConstraints)
    return [(F, S) for (F, S) in supported_constraints(model) if
        MOI.get(model, MOI.NumberOfConstraints{F,S}()) > 0]
end

#=
    Get constraint names
=#
function MOI.canget(::LinQuadOptimizer, ::MOI.ConstraintName, ::Type{<:MOI.ConstraintIndex})
    return true
end
function MOI.get(model::LinQuadOptimizer, ::MOI.ConstraintName, index::MOI.ConstraintIndex)
    if haskey(model.constraint_names, index)
        return model.constraint_names[index]
    else
        return ""
    end
end
function MOI.supports(::LinQuadOptimizer, ::MOI.ConstraintName, ::Type{<:MOI.ConstraintIndex})
    return true
end
function MOI.set!(model::LinQuadOptimizer, ::MOI.ConstraintName,
                  index::MOI.ConstraintIndex, name::String)
    if haskey(model.constraint_names_rev, name)
        if model.constraint_names_rev[name] != index
            error("Duplicate constraint name: $(name)")
        end
    elseif name != ""
        if haskey(model.constraint_names, index)
            # we're renaming an existing constraint
            old_name = model.constraint_names[index]
            delete!(model.constraint_names_rev, old_name)
        end
        model.constraint_names[index] = name
        model.constraint_names_rev[name] = index
    end
    return nothing
end

#=
    Get constraint by name
=#

# this covers the non-type-stable get(m, ConstraintIndex) case
function MOI.canget(model::LinQuadOptimizer, ::Type{MOI.ConstraintIndex}, name::String)
    return haskey(model.constraint_names_rev, name)
end
function MOI.get(model::LinQuadOptimizer, ::Type{<:MOI.ConstraintIndex}, name::String)
    return model.constraint_names_rev[name]
end

# this covers the type-stable get(m, ConstraintIndex{F,S}, name)::CI{F,S} case
function MOI.canget(model::LinQuadOptimizer, ::Type{FS}, name::String) where FS <: MOI.ConstraintIndex
    return haskey(model.constraint_names_rev, name) && typeof(model.constraint_names_rev[name]) == FS
end
function MOI.get(model::LinQuadOptimizer, ::Type{MOI.ConstraintIndex{F,S}}, name::String) where F where S
    model.constraint_names_rev[name]::MOI.ConstraintIndex{F,S}
end

"""
    CSRMatrix{T}

Matrix given in compressed sparse row (CSR) format.

`CSRMatrix` is analgous to the structure in Julia's `SparseMatrixCSC` but with
the rows and columns flipped. It contains three vectors:
 - `row_pointers` is a vector pointing to the start of each row in
    `columns` and `coefficients`;
 - `columns` is a vector of column numbers; and
 - `coefficients` is a vector of corresponding nonzero values.

The length of `row_pointers` is the number of rows in the matrix.

This struct is not a subtype of `AbstractSparseMatrix` as it is intended to be a
collection of the three vectors as they are required by solvers such as Gurobi.
It is not intended to be used for general computation.
"""
struct CSRMatrix{T}
    row_pointers::Vector{Int}
    columns::Vector{Int}
    coefficients::Vector{T}
    function CSRMatrix{T}(row_pointers, columns, coefficients) where T
        @assert length(columns) == length(coefficients)
        new(row_pointers, columns, coefficients)
    end
end

#=
    Below we add constraints.
=#

include("constraints/singlevariable.jl")
include("constraints/vectorofvariables.jl")
include("constraints/scalaraffine.jl")
include("constraints/vectoraffine.jl")
include("constraints/scalarquadratic.jl")


text:
#----------------------------------------------------------------------
# Boids algorithm visualization
# Based on pseudocode at http://www.kfish.org/boids/pseudocode.html
# Simulate flocking birds using simple rules.
# Initial version developed and contributed by Iain Dunning (@IainNZ)
# as part of the JuliaCon 2015 workshop
# TODO:
# - Add interactivity for parameters that govern behavior.
# - Add a force that attracts to mouse pointer.
#----------------------------------------------------------------------
using Reactive, GeometryTypes, GLAbstraction, GLPlot, GLVisualize, Colors, ColorTypes

typealias Position Point{2, Float32}
typealias Velocity Vec{2, Float32}

immutable Boids
    position::Vector{Position}
    velocity::Vector{Velocity}
end

# Create new boid at random location and velocity
Boids(n = 200) = Boids([rand(Position) for i=1:n], [rand(Velocity)/1000 for i=1:n])

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

function simulate!(t, boids)
    len = Float32(length(boids.position))
    boid_center = sum(boids.position)/len
    @inbounds for (i, vel) in enumerate(boids.velocity)
        boids.velocity[i] = vel + Vec(boid_center - boids.position[i])/600
    end

    # Force 2: Avoid others
    @inbounds for (i, boidpos) in enumerate(boids.position)
        avoidance = zero(Velocity)
        for (j, other_boidpos) in enumerate(boids.position)
            i == j && continue
            if norm(boidpos - other_boidpos) <= 0.2
                avoidance -= Vec(other_boidpos - boidpos)
            end
        end
        bvel = boids.velocity[i]
        boids.velocity[i] = bvel + (avoidance/1000f0)
    end
    @inbounds for (i, bvel) in enumerate(boids.velocity)
        perceived_vel = zero(Velocity)
        for (j, other_boidvel) in enumerate(boids.velocity)
            i == j && continue
            perceived_vel += other_boidvel
        end
        perceived_vel /= len-1f0
        boids.velocity[i] = bvel + (perceived_vel - bvel)/900f0
    end
    # Limit max velocity
    MAXVEL = 0.5f0
    @inbounds for (i, vel) in boids.velocity
        absvel = norm(vel)
        if absvel >= MAXVEL
            boids.velocity[i] = vel / absvel*MAXVEL
        end
    end
    # Update positions
    @inbounds for (i, vel) in enumerate(boids.velocity)
        boids.position[i] = boids.position[i] + Point(vel)
    end
    boids.position, boids.velocity
end
Base.clamp{T}(x::T) = clamp(x, zero(T), one(T))

to_color(velocities, mul) = RGBA{N0f8}[RGBA{N0f8}(clamp(velocity[1]*mul), clamp(velocity[2]*mul), 1.0, 0.8) for velocity in velocities]
function main()
    # Create population of boids
    boids = Boids()
    pv = map(simulate!, bounce(0:1000), Signal(boids))
    positions = map(first, pv)
    velocity = map(pv) do pv
        map(norm, pv[2])
    end
    glplot(
       (Circle, positions), scale = Vec2f0(0.05),
       intensity = velocity,
       color_map = GLVisualize.default(Vector{RGBA}),
       color_norm = Vec2f0(0, 0.02)
    )

end

x = main()


text:
import Base: show

function show(io::IO, s::Spectrum)
    println(io, "==== $(s.name) ($(s.expno)) ====")
    println(io, "FID size: $(length(s.fid))")
    println(io, "# processings: $(length(s.procs))")
    println(io, "Default proc. no.: $(s.default_proc)")
    out = PipeBuffer()
    for k in sort(collect(keys(s.procs)))
        show(out, s.procs[k])
        println(out, "")
    end
    for l in eachline(out)
        println(io, "    $l")
    end
end

function show(io::IO, p::ProcessedSpectrum)
    println(io, "---- Proc. no. $(p.procno) ----")
    println(io, strip(p.title))
    println(io, "Processed size: $(p["SI"]) points")
end


text:
#= test/cur.jl
=#

  @testset "cur" begin
    m = 128
    n =  64
    M = matrixlib(:fourier, rand(m), rand(n))
    Mh = M[1:n,1:n]; Mh += Mh';
    Ms = M[1:n,1:n]; Ms += transpose(Ms);
    opts = LRAOptions(maxdet_tol=0., sketch_randn_niter=1)

    @time for (t, s) in ((:none,                 :none ),
                         (:RandomGaussian,       :randn),
                         (:RandomSubset,         :sub  ),
                         (:SRFT,                 :srft ),
                         (:SparseRandomGaussian, :sprn ))
      opts.sketch = s
      for T in (Float32, Float64, ComplexF32, ComplexF64)
        let A
            println("  $t/$T")

            rtol = 5*eps(real(T))
            approx_rtol = 1000*rtol
            opts.rtol = rtol

            for (N, ρ, q) in ((M, m, n), (Mh, n, n), (Ms, n, n))
              A = convert(Array{T}, T <: Real ? real(N) : N)

              U  = curfact(A, opts)
              F  = CUR(A, U)
              @test norm(A - Matrix(F)) < approx_rtol*norm(A)

              xp = rand(T, ρ)
              xq = rand(T, q)
              y = A  *xq; @test norm(y - F  *xq) < approx_rtol*norm(y)
              y = A' *xp; @test norm(y - F' *xp) < approx_rtol*norm(y)
              y = transpose(A)*xp; @test norm(y - transpose(F)*xp) < approx_rtol*norm(y)

              Bp = rand(T, ρ, ρ)
              Bq = rand(T, q, q)
              C = A   *Bq  ; @test norm(C - F   *Bq  ) < approx_rtol*norm(C)
              C = A   *Bq' ; @test norm(C - F   *Bq' ) < approx_rtol*norm(C)
              C = A   *transpose(Bq); @test norm(C - F   *transpose(Bq)) < approx_rtol*norm(C)
              C = A'  *Bp  ; @test norm(C - F'  *Bp  ) < approx_rtol*norm(C)
              C = A'  *Bp' ; @test norm(C - F'  *Bp' ) < approx_rtol*norm(C)
              C = transpose(A) *Bp  ; @test norm(C - transpose(F) *Bp  ) < approx_rtol*norm(C)
              C = transpose(A) *transpose(Bp); @test norm(C - transpose(F) *transpose(Bp)) < approx_rtol*norm(C)
              C = Bp  *A   ; @test norm(C - Bp  *F   ) < approx_rtol*norm(C)
              C = Bq  *A'  ; @test norm(C - Bq  *F'  ) < approx_rtol*norm(C)
              C = Bq  *transpose(A) ; @test norm(C - Bq  *transpose(F) ) < approx_rtol*norm(C)
              C = Bp' *A   ; @test norm(C - Bp' *F   ) < approx_rtol*norm(C)
              C = Bq' *A'  ; @test norm(C - Bq' *F'  ) < approx_rtol*norm(C)
              C = transpose(Bp)*A   ; @test norm(C - transpose(Bp)*F   ) < approx_rtol*norm(C)
              C = transpose(Bq)*transpose(A) ; @test norm(C - transpose(Bq)*transpose(F) ) < approx_rtol*norm(C)
          end
        end
      end
    end
  end


text:
"""
    AbstractSpan

Section of a parent String with a specific meaning for Franklin. All subtypes of
`AbstractBlock` must have an `ss` field corresponding to the substring associated to the
block. This field is necessarily of type `SubString{String}`.
"""
abstract type AbstractSpan end

from(s::AbstractSpan)::Int             = from(s.ss::SS)
to(s::AbstractSpan)::Int               = to(s.ss::SS)
parent_string(s::AbstractSpan)::String = parent_string(s.ss::SS)
content(s::AbstractSpan)::SS           = s.ss
content(s::SS)::SS                     = s
content(s::String)::SS                 = subs(s)
Base.isempty(o::AbstractSpan)::Bool    = isempty(strip(o.ss))

"""
    Token <: AbstractSpan

A token is a subtype of `AbstractSpan` which typically determines the start or end of
a block. It can also be used for special characters.
"""
struct Token <: AbstractSpan
    name::Symbol
    ss::SS
end

is_eos(t::Token) = t.name == :EOS
to(t::Token)     = ifelse(is_eos(t), from(t.ss), to(t.ss))

const EMPTY_TOKEN      = Token(:NONE, subs(""))
const EMPTY_TOKEN_SVEC = @view (Token[])[1:0]

"""
    Block <: AbstractSpan

Blocks are defined by an opening and a closing `Token`, they may be nested. For instance
braces block are formed of an opening `{` and a closing `}`.
"""
struct Block <: AbstractSpan
    name::Symbol
    open::Token
    close::Token
    ss::SS
    inner_tokens::SubVector{Token}
end

function Block(n::Symbol, p::Pair{Token, Token}, it=EMPTY_TOKEN_SVEC)
    o, c = p.first, p.second
    ss = subs(parent_string(o), from(o), to(c))
    return Block(n, o, c, ss, it)
end

function Block(n::Symbol, s::SubVector{Token})
    it = @view s[2:end-1]
    return Block(n, s[1] => s[end], it)
end

function Block(n::Symbol, ss::SS, it=EMPTY_TOKEN_SVEC)
    return Block(n, EMPTY_TOKEN, EMPTY_TOKEN, ss, it)
end

TokenBlock(t::Token) = Block(t.name, t, EMPTY_TOKEN, t.ss, EMPTY_TOKEN_SVEC)

"""
    TextBlock

Spans of text which should be left to the fallback engine (such as CommonMark for
instance). Text blocks can also have inner tokens that are non-block delimiters such as
emojis or html entities.
"""
function TextBlock(ss::SS, it=EMPTY_TOKEN_SVEC)::Block
    isempty(it) && return Block(:TEXT, ss)
    fss = from(ss)
    tss = to(ss)
    i = findfirst(t -> fss <= from(t), it)
    j = findlast(t -> to(t) <= tss && !is_eos(t), it)
    any(isnothing, (i, j)) && return Block(:TEXT, ss)
    inner_tokens = @view it[i:j]
    return Block(:TEXT, ss, inner_tokens)
end

"""
    content(block)

Return the content of a `Block`, for instance the content of a `{...}` block would be
`...`. Note EOS is a special '0 length' case to  deal with the fact that a text can end
with a token (which would then be an overlapping token and an EOS).
"""
function content(b::Block)::SS
    b.name == :TEXT && return b.ss
    # find the relevant range of the parent string
    s = parent_string(b.ss)
    t = from(b.close)
    idxo = nextind(s, to(b.open))
    idxc = ifelse(is_eos(b.close), t, prevind(s, t))
    # return the substring corresponding to the range
    return subs(s, idxo, idxc)
end

"""
    BlockTemplate

Template for a block to find. A block goes from a token with a given `opening` name to
one of several possible `closing` names. Blocks can allow or disallow nesting. For
instance brace blocks can be nested `{.{.}.}` but not comments.
When nesting is enabled, Franklin will try to find the closing token taking into account
the balance in opening-closing tokens.
"""
struct BlockTemplate
    name::Symbol
    opening::Symbol
    closing::NTuple{N, Symbol} where N
    nesting::Bool
end

BlockTemplate(n, o, c::Symbol, ne) = BlockTemplate(n, o, (c,), ne)
BlockTemplate(a...; nesting=false) = BlockTemplate(a..., nesting)

const NO_CLOSING = (:NONE,)

SingleTokenBlockTemplate(name::Symbol) = BlockTemplate(name, name, NO_CLOSING, false)

"""
    Group <: AbstractSpan

A Group contains 1 or more more Blocks and will map to either a Paragraph or
something else like a code block.
"""
struct Group <: AbstractSpan
    role::Symbol
    blocks::Vector{Block}
    ss::SS
end

function Group(blocks::Vector{Block}; role::Symbol=:none)
    ps = parent_string(first(blocks))
    return Group(
        role,
        blocks,
        subs(ps, from(first(blocks)), to(last(blocks)))
    )
end

Group(block::Block; kw...) = Group([block]; kw...)


text:
@testset "PointCloud Model tests" begin
    x_test = rand(Float32, 256, 3, 1)   #Reduced npoints to speed up tests
    npoints = size(x_test,1)

    @testset "PointNet Model tests" begin
        MODEL = PointNet
        @info "Testing $(MODEL)..."
        for num_classes in [10, 40]
            @testset "PointNet num_classes: $(num_classes)" begin
                model = MODEL(num_classes)

                y_test = model(x_test)

                # Testing Forward Propagation
                @test y_test isa AbstractArray
                @test size(y_test) == (num_classes,1)

                # Testing Backward Propagation
                @test gradtest(x_test -> model(x_test), x_test)
            end
        end
    end

    @testset "DGCNN Model tests" begin
        MODEL = DGCNN
        @info "Testing $(MODEL)..."
        for num_classes in [10, 40]
            @testset "DGCNN num_classes: $(num_classes)" begin
                model = MODEL(num_classes,10,npoints)

                y_test = model(x_test)

                # Testing Forward Propagation
                @test y_test isa AbstractArray
                @test size(y_test) == (num_classes,1)

                # Testing Backward Propagation
                @test gradtest(x_test -> model(x_test), x_test)
            end
        end
    end

end


text:
# Undergraduate Student: Arturo Burgos
# Professor: João Rodrigo Andrade
# Federal University of Uberlândia - UFU, Fluid Mechanics Laboratory - MFLab, Block 5P, Uberlândia, MG, Brazil

using Printf # Just the same as import a module
using Statistics
using CPUTime

A = rand(1:50,3,3,3)

println("\n\nThe A[1,1,1] is: ",A[1,1,1])
println("The A[1,1,2] is: ",A[1,1,2])
println("The A[1,1,3] is: ",A[1,1,3])
println("\n\nThe A[3,3,1] is: ",A[3,3,1])
println("The A[3,3,2] is: ",A[3,3,2])
println("The A[3,3,3] is: ",A[3,3,3])

function loop(mtx)

    for i = 1:3
        for j=1:3
            mtx[i, j, 1] = mtx[i, j, 2]
            mtx[i, j, 3] = mtx[i, j, 1]
        end
    end
    return mtx # retornar algo à variavel alocada, no caso "A"
end

A = loop(A)

println("\n\n\nThe new Matrix is: \n\n")

println("\n\nThe A[1,1,1] is: ",A[1,1,1])
println("The A[1,1,2] is: ",A[1,1,2])
println("The A[1,1,3] is: ",A[1,1,3])

println("\n\nThe A[3,3,1] is: ",A[3,3,1])
println("The A[3,3,2] is: ",A[3,3,2])
println("The A[3,3,3] is: ",A[3,3,3])

@time @CPUtime loop(A)

max_stars_repo_name: bangyen/Comparing-Languages

max_issues_repo_name: bangyen/Comparing-Languages

max_forks_repo_name: bangyen/Comparing-Languages


text:
using MarkovTools
using Base.Test

tests = [
    "test_permutations",
    "prng/test_input"
]

for t in tests
    test_fn = "$t.jl"
    include(test_fn)
end

max_stars_repo_name: PawelPamula/MarkovTools

max_issues_repo_name: PawelPamula/MarkovTools

max_forks_repo_name: PawelPamula/MarkovTools


text:
__precompile__()

module LightGBM
using Libdl
using Dates
import StatsBase

# if LIGHTGBM_PATH is difined, automatical define LIGHTGBM_PATH
function lgbm_librarysetup()
    try
        println("ENV[\"LIGHTGBM_PATH\"] is ",ENV["LIGHTGBM_PATH"])
        return ENV["LIGHTGBM_PATH"]
    catch
        println("Don't find ENV[\"LIGHTGBM_PATH\"].LightGBM library setup start...")
        global ligthgbmpath=joinpath(abspath(joinpath(dirname(Base.find_package("LightGBM")), "..")),"deps/usr/lib")
        if Sys.islinux()
            global prefix=joinpath(ligthgbmpath,"lib_lightgbm.so")
        elseif Sys.iswindows()
            global prefix=joinpath(ligthgbmpath,"lib_lightgbm.dll")
        elseif Sys.isapple()
            global prefix=joinpath(ligthgbmpath,"lib_lightgbm.dylib")
        else
            global prefix=""
        end

        if !isfile(prefix)
            println("Instaling library")
            include(joinpath(abspath(joinpath(dirname(Base.find_package("LightGBM")), "..")),"deps/build.jl"))
        end

        if isfile(prefix)
            println("Set library path:",ligthgbmpath)
            return ligthgbmpath
        else
            println("Not find LightGBM library")
        end
    end
end

# pre initialization
ENV["LIGHTGBM_PATH"] = lgbm_librarysetup();

function __init__()
    println("Start __init__")

    try
        println("LIGHTGBM_PATH:",ENV["LIGHTGBM_PATH"]);
    catch
        println("LIGHTGBM_PATH setup1...")
        ENV["LIGHTGBM_PATH"] = lgbm_librarysetup()
        global LGBM_library = Libdl.find_library(["lib_lightgbm.so", "lib_lightgbm.dll",
            "lib_lightgbm.dylib"], [ENV["LIGHTGBM_PATH"]])
    end

    if !haskey(ENV, "LIGHTGBM_PATH")
        error("Environment variable LIGHTGBM_PATH not found. ",
            "Set this variable to point to the LightGBM directory prior to loading LightGBM.jl ",
            "(e.g. `ENV[\"LIGHTGBM_PATH\"] = \"../LightGBM\"`).")
    end

    println("Finished __init__()")
end

const LGBM_library = find_library(["lib_lightgbm.so", "lib_lightgbm.dll", "lib_lightgbm.dylib"], [ENV["LIGHTGBM_PATH"]])

include("wrapper.jl")
include("estimators.jl")
include("utils.jl")
include("fit.jl")
include("predict.jl")
include("cv.jl")
include("search_cv.jl")
include("LightGBM-util2.jl")

export fit, predict, cv, search_cv, savemodel, loadmodel
export LGBMEstimator, LGBMRegression, LGBMBinary, LGBMMulticlass
export metaformattedclassresult, metaformattedclassresult, formattedclassfit, predict2

end # module LightGBM


text:
NetCDF.open(fn1, "v1", mode = NC_WRITE) do v1
    @test isa(v1, NetCDF.NcVar{Float64,3,6})
    @test v1[1, 1, 1] == x1[1, 1, 1]

    @test dropdims(v1[2, :, 2], dims = (1, 3)) == x1[2, :, 2]

    @test dropdims(v1[:, 3, 1], dims = (2, 3)) == x1[:, 3, 1]
    inds = bitrand(size(x1))
    @test v1[inds] == x1[inds]

    xnew = rand(Float64, size(x1))
    v1[:, :, :] = xnew
    @test v1 == xnew
end


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.11.4/agnostic/x64/nccl_2.11.4-1+cuda11.4_x86_64.txz",
                      "639a86d7f90a0b9f8143ef2044b9cdccd7de51b28cccce2f9852d3bdd0e4d114")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.11.4/agnostic/ppc64le/nccl_2.11.4-1+cuda11.4_ppc64le.txz",
                      "6347d56c6e13fc321858fd8779741c62ac0bf55a13ab57d39e2780d62181a6b4")],
    Platform("aarch64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.11.4/agnostic/aarch64sbsa/nccl_2.11.4-1+cuda11.4_aarch64.txz",
                      "4bef02db97ef45146b80e3d20218198f57cd24f115624f50c5eca06fff3f216f")],
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function read_field_datatypes(f::JLDFile, attrs::Vector{ReadAttribute})
    for attr in attrs
        if attr.name == :field_datatypes
            return read_attr_data(f, attr, ReferenceDatatype(),
                                  ReadRepresentation{RelOffset,RelOffset}())
        end
    end
    RelOffset[]
end

function check_empty(attrs::Vector{ReadAttribute})
    for attr in attrs
        if attr.name == :empty
            return true
        end
    end
    false
end

# jltype is the inverse of h5type, providing a ReadRepresentation for an
# H5Datatype. We handle committed datatypes here, and other datatypes below.
function jltype(f::JLDFile, cdt::CommittedDatatype)
    haskey(f.h5jltype, cdt) && return f.h5jltype[cdt]::ReadRepresentation
    dt, attrs = read_committed_datatype(f, cdt)

    julia_type_attr = nothing
    written_type_attr = nothing
    for attr in attrs
        if attr.name == :julia_type
            julia_type_attr = attr
        elseif attr.name == :written_type
            written_type_attr = attr
        end
    end

    if isa(julia_type_attr, Nothing)
        throw(InvalidDataException())
    end
    julia_type_attr = julia_type_attr::ReadAttribute

    # If type of datatype is this datatype, then this is the committed
    # datatype that describes a datatype
    if h5offset(f, julia_type_attr.datatype_offset) == cdt.header_offset
        # Verify that the datatype matches our expectations
        if dt != H5TYPE_DATATYPE
            error("""The HDF5 datatype representing a Julia datatype does not match
                     the expectations of this version of JLD.

                     You may need to update JLD to read this file.""")
        end
        f.jlh5type[DataType] = cdt
        f.datatypes[cdt.index] = dt
        return (f.h5jltype[cdt] = ReadRepresentation{DataType, DataTypeODR()}())
    end

    datatype = read_attr_data(f, julia_type_attr)
    if written_type_attr !== nothing
        # Custom serialization
        readas = datatype
        datatype = read_attr_data(f, written_type_attr)
        if isa(readas, UnknownType)
            @warn("custom serialization of $(typestring(readas))" *
                  " encountered, but the type does not exist in the workspace; the data will be read unconverted")
            rr = (constructrr(f, datatype, dt, attrs)::Tuple{ReadRepresentation,Bool})[1]
            canonical = false
        else
            rr, canonical = constructrr(f, datatype, dt, attrs)::Tuple{ReadRepresentation,Bool}
            rrty = typeof(rr)
            rr = ReadRepresentation{readas, CustomSerialization{rrty.parameters[1], rrty.parameters[2]}}()
            canonical = canonical && writeas(readas) === datatype
        end
    else
        rr, canonical = constructrr(f, datatype, dt, attrs)::Tuple{ReadRepresentation,Bool}
    end

    canonical && (f.jlh5type[datatype] = cdt)
    f.datatypes[cdt.index] = dt
    f.h5jltype[cdt] = rr
end

# Constructs a ReadRepresentation for a given opaque (bitstype) type
function constructrr(::JLDFile, T::DataType, dt::BasicDatatype, attrs::Vector{ReadAttribute})
    dt.class == DT_OPAQUE || throw(UnsupportedFeatureException())
    if T.size == dt.size && isempty(T.types)
        (ReadRepresentation{T,T}(), true)
    else
        empty = check_empty(attrs)
        if empty
            if !hasdata(T)
                (ReadRepresentation{T,nothing}(), true)
            else
                @warn("$T has $(T.size*8) bytes, but written type was empty; reconstructing")
                reconstruct_bitstype(T.name.name, dt.size, empty)
            end
        else
            if isempty(T.types)
                @warn("primitive type $T has $(T.size*8) bits, but written type has $(dt.size*8) bits; reconstructing")
            else
                @warn("$T is a non-primitive type, but written type is a primitive type with $(dt.size*8) bits; reconstructing")
            end
            reconstruct_bitstype(T.name.name, dt.size, empty)
        end
    end
end

"""
constructrr(f::JLDFile, T::DataType, dt::CompoundType, attrs::Vector{ReadAttribute},
            hard_failure::Bool=false)

Constructs a ReadRepresentation for a given type. This is the generic method for all
types not specially handled below.

If hard_failure is true, then throw a TypeMappingException instead of attempting
reconstruction. This helps in cases where we can't know if reconstructed parametric types
will have a matching memory layout without first inspecting the memory layout.
"""
struct TypeMappingException <: Exception; end
function constructrr(f::JLDFile, T::DataType, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute},
                     hard_failure::Bool=false)
    field_datatypes = read_field_datatypes(f, attrs)

    # If read type is not a leaf type, reconstruct
    if !isconcretetype(T)
        @warn("read type $T is not a leaf type in workspace; reconstructing")
        return reconstruct_compound(f, string(T), dt, field_datatypes)
    end

    # Map names in dt to their indices
    dtnames = Dict{Symbol,Int}()
    for i = 1:length(dt.names)
        dtnames[dt.names[i]] = i
    end
    mapped = falses(length(dt.names))

    offsets = Vector{Int}(undef, length(T.types))
    types = Vector{Any}(undef, length(T.types))
    odrs = Vector{Any}(undef, length(T.types))
    fn = fieldnames(T)
    samelayout = isbitstype(T) && T.size == dt.size
    dtindex = 0
    for i = 1:length(T.types)
        wstype = T.types[i]
        writtenas = writeas(wstype)
        if !hasfielddata(writtenas)
            types[i] = wstype
            odrs[i] = wstype === writtenas ? nothing : CustomSerialization{writtenas,nothing}
            offsets[i] = dtindex == 0 ? 0 : (dt.offsets[dtindex] + dt.members[dtindex].size::UInt32)
        else
            if !haskey(dtnames, fn[i])
                hard_failure && throw(TypeMappingException())
                @warn("saved type $T is missing field $(fn[i]) in workspace type; reconstructing")
                return reconstruct_compound(f, string(T), dt, field_datatypes)
            end

            dtindex = dtnames[fn[i]]
            if !isempty(field_datatypes) && (ref = field_datatypes[dtindex]) != NULL_REFERENCE
                dtrr = jltype(f, f.datatype_locations[ref])
            else
                dtrr = jltype(f, dt.members[dtindex])
            end

            readtype, odrtype = typeof(dtrr).parameters

            if typeintersect(readtype, wstype) === Union{} &&
               !hasmethod(convert, Tuple{Type{wstype}, readtype})
                # Saved type does not match type in workspace and no
                # convert method exists, so we definitely need to reconstruct.
                hard_failure && throw(TypeMappingException())
                @warn("saved type $T has field $(fn[i])::$(readtype)" *
                      ", but workspace type has field $(fn[i])::$(wstype)" *
                      ", and no applicable convert method exists; reconstructing")
                return reconstruct_compound(f, string(T), dt, field_datatypes)
            end

            types[i] = readtype
            odrs[i] = odrtype
            offsets[i] = dt.offsets[dtindex]

            # The on disk representation of T can only be the same as in memory
            # if the offsets are the same, field type on disk (readtype) and in memory (wstype)
            # are the same and if no CustomSerialization is involved
            samelayout = samelayout &&
                offsets[i] == fieldoffset(T, i) &&
                types[i] === wstype &&
                # An OnDiskRepresentation as odr means that something "fixable" went wrong
                # for this field
                !(odrs[i] isa OnDiskRepresentation) &&
                !(odrs[i] <: CustomSerialization)

            mapped[dtindex] = true
        end
    end

    if !all(mapped)
        @warn("the following fields are present in type $T" *
              " saved in the file but not present in the type the workspace:\n\n" *
              "$(join(dt.names[.!mapped], "\n"))," *
              "\n\nData in these fields will not be accessible")
    end

    if samelayout
        (ReadRepresentation{T,T}(), true)
    else
        wodr = odr(T)
        # This should theoretically be moved inside the if statement, but then it returns
        # the wrong result due to a bug in type inference on 0.6
        typeof_wodr = typeof(wodr)
        offsets = (offsets...,)
        if wodr isa OnDiskRepresentation
            odr_offsets = typeof_wodr.parameters[1]
            odr_types = typeof_wodr.parameters[2].parameters
            odr_h5types = typeof_wodr.parameters[3].parameters
            tequal = length(odr_types) == length(types)
            if tequal
                for i = 1:length(types)
                    if !(odr_types[i] <: types[i]) || odr_h5types[i] != odrs[i]
                        tequal = false
                        break
                    end
                end
                if tequal && odr_offsets == offsets
                    # This should not be necessary, but type inference mistakenly changes
                    # the value of wodr here
                    wodr = typeof_wodr()
                    return (ReadRepresentation{T,wodr}(), true)
                end
            end
        end
        return (ReadRepresentation{T,OnDiskRepresentation{offsets, Tuple{types...}, Tuple{odrs...}}()}(), false)
    end
end

function constructrr(f::JLDFile, T::UnionAll, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute},
                     hard_failure::Bool=false)
    @warn("read type $T is not a leaf type in workspace; reconstructing")
    return reconstruct_compound(f, string(T), dt, read_field_datatypes(f, attrs))
end

# Find types in modules
# returns the result of searching for the type in the specified module m
function _resolve_type_singlemodule(::ReadRepresentation{T,DataTypeODR()},
                                    m,
                                    parts,
                                    mypath,
                                    hasparams::Bool,
                                    params) where T
    for part in parts
        sym = Symbol(part)
        if !isa(m, Module) || !isdefined(m, sym)
            return hasparams ? UnknownType(mypath, params) : UnknownType(mypath)
        end
        m = getfield(m, sym)
    end
    if !isa(m, DataType) && !isa(m, UnionAll)
        return hasparams ? UnknownType(mypath, params) : UnknownType(mypath)
    end
    return m
end

_is_not_unknown_type(x::UnknownType) = false
_is_not_unknown_type(x) = true

function _resolve_type(rr::ReadRepresentation{T,DataTypeODR()},
                       f::JLDFile,
                       ptr::Ptr,
                       header_offset::RelOffset,
                       mypath,
                       hasparams::Bool,
                       params) where T
    parts = split(mypath, '.')
    modules = vcat([Main], collect(keys(Base.module_keys)), stdlibmodules(Main))
    unique!(modules)
    for mod in modules
        resolution_attempt = _resolve_type_singlemodule(rr,
                                                        mod,
                                                        parts,
                                                        mypath,
                                                        hasparams,
                                                        params)
        if _is_not_unknown_type(resolution_attempt)
            return resolution_attempt
        end
    end
    return hasparams ? UnknownType(mypath, params) : UnknownType(mypath)
end

function types_from_refs(f::JLDFile, ptr::Ptr)
    # Test for a potential null pointer indicating an empty array
    isinit = jlunsafe_load(convert(Ptr{UInt32}, ptr)) != 0
    unknown_params = false
    if isinit
        refs = jlconvert(ReadRepresentation{RelOffset, Vlen{RelOffset}}(), f, ptr, NULL_REFERENCE)
        params =  Any[let
            # If the reference is to a committed datatype, read the datatype
            nulldt = CommittedDatatype(UNDEFINED_ADDRESS, 0)
            cdt = get(f.datatype_locations, ref, nulldt)
            res = cdt !== nulldt ? (typeof(jltype(f, cdt)::ReadRepresentation)::DataType).parameters[1] : load_dataset(f, ref)
            unknown_params = unknown_params || isa(res, UnknownType)
            res
        end for ref in refs]
        return params, unknown_params
    end
    return [], unknown_params
end

# Read a type. Returns an instance of UnknownType if the type or parameters
# could not be resolved.
function jlconvert(rr::ReadRepresentation{T,DataTypeODR()},
                   f::JLDFile,
                   ptr::Ptr,
                   header_offset::RelOffset) where T

    params, unknown_params = types_from_refs(f, ptr+odr_sizeof(Vlen{UInt8}))
    # For cross-platform compatibility convert integer type parameters to system precision
    params = [p isa Union{Int64,Int32} ? Int(p) : p for p in params]
    hasparams = !isempty(params)
    mypath = String(jlconvert(ReadRepresentation{UInt8,Vlen{UInt8}}(), f, ptr, NULL_REFERENCE))

    if mypath in keys(f.typemap)
        m = f.typemap[mypath]
    else
        m = _resolve_type(rr, f, ptr, header_offset, mypath, hasparams, hasparams ? params : nothing)
        m isa UnknownType && return m
    end

    if hasparams
        try
            m = m{params...}
        catch e
            return UnknownType(m, params)
        end
    elseif m === Tuple
        # Need to instantiate with no parameters, since Tuple is really
        # Tuple{Vararg{Any}}
        m = Tuple{}
    end
    track_weakref!(f, header_offset, m)
    return m
end

constructrr(::JLDFile, ::Type{T}, dt::CompoundDatatype, ::Vector{ReadAttribute}) where {T<:DataType} =
    dt == H5TYPE_DATATYPE ? (ReadRepresentation{DataType,DataTypeODR()}(), true) :
                            throw(UnsupportedFeatureException())

## Type reconstruction

module ReconstructedTypes end

isreconstructed(x) = isreconstructed(typeof(x))

function isreconstructed(::Type{T}) where T
    return supertype(T) == Any && parentmodule(T) == ReconstructedTypes
end

function reconstruct_bitstype(name::Union{Symbol,String}, size::Integer, empty::Bool)
    sym = gensym(name)
    Core.eval(ReconstructedTypes, empty ? :(struct $(sym) end) : :(primitive type $(sym) $(Int(size)*8) end))
    T = getfield(ReconstructedTypes, sym)
    (ReadRepresentation{T, empty ? nothing : T}(), false)
end

function constructrr(f::JLDFile, unk::UnknownType, dt::BasicDatatype,
                     attrs::Vector{ReadAttribute})
    @warn("type $(typestring(unk)) does not exist in workspace; reconstructing")
    reconstruct_bitstype(typestring(unk), dt.size, check_empty(attrs))
end

# Convert an ordinary type or an UnknownType to a corresponding string. This is
# only used to create gensymmed names for reconstructed types.
function typestring(T::UnknownType)
    tn = IOBuffer()
    print(tn, T.name)
    if isdefined(T, :parameters)
        write(tn, '{')
        for i = 1:length(T.parameters)
            x = T.parameters[i]
            i != 1 && write(tn, ',')
            if isa(x, UnknownType)
                write(tn, typestring(x))
            else
                print(tn, x)
            end
        end
        write(tn, '}')
    end
    String(take!(tn))
end

function constructrr(f::JLDFile, unk::UnknownType{DataType}, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute})
    field_datatypes = read_field_datatypes(f, attrs)
    if unk.name == Tuple
        # For a tuple with unknown fields, we should reconstruct the fields
        rodr = reconstruct_odr(f, dt, field_datatypes)
        # This is a "pseudo-RR" since the tuple is not fully parametrized, but
        # the parameters must depend on the types actually encoded in the file
        (ReadRepresentation{Tuple,rodr}(), false)
    else
        @warn("read type $(typestring(unk)) was parametrized, but type " *
              "$(unk.name) in workspace is not; reconstructing")
        reconstruct_compound(f, typestring(unk), dt, field_datatypes)
    end
end

"""
    behead(T)

Given a UnionAll type, recursively eliminates the `where` clauses
"""
behead(T::UnionAll) = behead(T.body)
behead(@nospecialize T) = T

function constructrr(f::JLDFile, unk::UnknownType{UnionAll}, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute})
    field_datatypes = read_field_datatypes(f, attrs)
    body = behead(unk.name)
    if length(body.parameters) != length(unk.parameters)
        @warn("read type $(typestring(unk)) has a different number of parameters from type " *
              "$(unk.name) in workspace; reconstructing")
        reconstruct_compound(f, typestring(unk), dt, field_datatypes)
    else
        params = copy(unk.parameters)
        for i = 1:length(params)
            if isa(params[i], UnknownType)
                param = body.parameters[i]::TypeVar
                params[i] = param.ub
            end
        end

        # Try to construct the rr for the relaxed type. On failure, fall back to
        # reconstruct_compound.
        local T
        try
            T = unk.name{params...}
        catch err
            @warn("type parameters for $(typestring(unk)) do not match type $(unk.name) in workspace; reconstructing")
            return reconstruct_compound(f, typestring(unk), dt, field_datatypes)
        end

        try
            (rr,) = constructrr(f, T, dt, attrs, true)
            @warn("some parameters could not be resolved for type $(typestring(unk)); reading as $T")
            return (rr, false)
        catch err
            !isa(err, TypeMappingException) && rethrow(err)
            @warn("some parameters could not be resolved for type $(typestring(unk)); reconstructing")
            return reconstruct_compound(f, typestring(unk), dt, field_datatypes)
        end
    end
end

# An UnknownType for which we were not able to resolve the name.
function constructrr(f::JLDFile, unk::UnknownType{String}, dt::CompoundDatatype,
                     attrs::Vector{ReadAttribute})
    ts = typestring(unk)
    @warn("type $ts does not exist in workspace; reconstructing")
    reconstruct_compound(f, ts, dt, read_field_datatypes(f, attrs))
end

# Reconstruct the ODR of a type from the CompoundDatatype and field_datatypes
# attribute
function reconstruct_odr(f::JLDFile, dt::CompoundDatatype,
                         field_datatypes::Vector{RelOffset})
    # Get the type and ODR information for each field
    types = Vector{Any}(undef, length(dt.names))
    h5types = Vector{Any}(undef, length(dt.names))
    for i = 1:length(dt.names)
        if !isempty(field_datatypes) && (ref = field_datatypes[i]) != NULL_REFERENCE
            dtrr = jltype(f, f.datatype_locations[ref])
        else
            dtrr = jltype(f, dt.members[i])
        end
        types[i], h5types[i] = typeof(dtrr).parameters
    end
    return OnDiskRepresentation{(dt.offsets...,), Tuple{types...}, Tuple{h5types...}}()
end

# Reconstruct type that is a "lost cause": either we were not able to resolve
# the name, or the workspace type has additional fields, or cannot convert
# fields to workspace types
function reconstruct_compound(f::JLDFile, T::String, dt::H5Datatype,
                              field_datatypes::Union{Vector{RelOffset},Nothing})
    rodr = reconstruct_odr(f, dt, field_datatypes)
    types = typeof(rodr).parameters[2].parameters

    # Now reconstruct the type
    reconname = gensym(T)
    Core.eval(ReconstructedTypes,
              Expr(:struct, false, reconname,
                   Expr(:block, Any[ Expr(Symbol("::"), dt.names[i], types[i]) for i = 1:length(types) ]...,
                        # suppress default constructors, plus a bogus `new()` call to make sure
                        # ninitialized is zero.
                        Expr(:if, false, Expr(:call, :new)))))
    T = getfield(ReconstructedTypes, reconname)

    (ReadRepresentation{T,rodr}(), false)
end

# At present, we write Union{} as an object of Core.TypeofBottom. The method above
# basically works, but `Expr(:new, Type{Union{}})` is a bit weird and causes problems for
# inference. Better to define a separate method.
jlconvert(::ReadRepresentation{Core.TypeofBottom,nothing}, f::JLDFile, ptr::Ptr,
          header_offset::RelOffset) = Union{}

# This jlconvert method handles compound types with padding or references
@generated function jlconvert(::ReadRepresentation{T,S}, f::JLDFile, ptr::Ptr,
                              header_offset::RelOffset) where {T,S}
    isa(S, DataType) && return :(convert(T, jlunsafe_load(convert(Ptr{S}, ptr))))
    @assert isa(S, OnDiskRepresentation)

    offsets = typeof(S).parameters[1]
    types = typeof(S).parameters[2].parameters
    odrs = typeof(S).parameters[3].parameters

    blk = Expr(:block)
    args = blk.args
    # Separate treatment for mutable structs for better readability
    if ismutabletype(T)
        push!(args, quote
            obj = $(Expr(:new, T))
            track_weakref!(f, header_offset, obj)
        end)
        fn = fieldnames(T)
        for i = 1:length(types)
            offset = offsets[i]
            rtype = types[i]
            odr = odrs[i]

            rr = ReadRepresentation{rtype,odr}()

            fni = QuoteNode(fn[i])
            ttype = T.types[i]
            if odr === nothing
                # Type is not stored or single instance
                newi = Expr(:new, ttype)
                push!(args, :(setfield!(obj, $fni, $newi)))
            else
                loadfield = :(setfield!(obj, $fni, rconvert($ttype, jlconvert($rr, f, ptr+$offset, NULL_REFERENCE)::$rtype)::$ttype))
                if jlconvert_canbeuninitialized(rr)
                    push!(args, :(jlconvert_isinitialized($rr, ptr+$offset) && $(loadfield)))
                else
                    push!(args, loadfield)
                end
            end
        end

        push!(args, (:obj))
        return blk
    end
    if isbitstype(T)
        # For bits types, we should always inline, because otherwise we'll just
        # pass a lot of crap around in registers
        push!(args, Expr(:meta, :inline))
    end
    fsyms = []
    fn = T === Tuple ? [Symbol(i) for i = 1:length(types)] : fieldnames(T)
    for i = 1:length(types)
        offset = offsets[i]
        rtype = types[i]
        odr = odrs[i]

        fsym = Symbol("field_", fn[i])
        push!(fsyms, fsym)

        rr = ReadRepresentation{rtype,odr}()

        if odr === nothing
            # Type is not stored or single instance
            if T.types[i] == Union{}
                # This cannot be defined
                @assert !ismutabletype(T)
                push!(args, Expr(:return, Expr(:new, T, fsyms[1:i-1]...)))
                return blk
            else
                newi = Expr(:new, T.types[i])
                push!(args, :($fsym = $newi))
            end
        else
            if jlconvert_canbeuninitialized(rr)
                push!(args, quote
                    if !jlconvert_isinitialized($rr, ptr+$offset)
                        $(if T <: Tuple || i <= ninitialized(T)
                            # Reference must always be initialized
                            :(throw(UndefinedFieldException(T,$(QuoteNode(fn[i])))))
                        else
                            Expr(:return, Expr(:new, T, fsyms[1:i-1]...))
                        end)
                    end
                end)
            end
            if T === Tuple
                # Special case for reconstructed tuples, where we don't know the
                # field types in advance
                push!(args, :($fsym = jlconvert($rr, f, ptr+$offset, NULL_REFERENCE)::$rtype))
            else
                ttype = T.types[i]
                push!(args, :($fsym = rconvert($ttype, jlconvert($rr, f, ptr+$offset, NULL_REFERENCE)::$rtype)::$ttype))
            end
        end
    end

    push!(args, T <: Tuple ? Expr(:tuple, fsyms...) : Expr(:new, T, fsyms...))

    blk
end


text:
using RecurrenceAnalysis
using Test
using Downloads

# Download some test timeseries
tsfolder = joinpath(@__DIR__, "timeseries")
todownload1 = ["$n.csv" for n in 1:4]
todownload = ["test_time_series_lorenz_standard_N_10000_multivariate.csv", "test_time_series_roessler_N_10000_multivariate.csv"]
append!(todownload, todownload1)
repo = "https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/timeseries"
mkpath(tsfolder)
for a in todownload
    Downloads.download(repo*"/"*a, joinpath(tsfolder, a))
end

ti = time()

@testset "RecurrenceAnalysis tests" begin
    include("dynamicalsystems.jl")
    include("smallmatrix.jl")
    include("skeletontest.jl")
    include("deprecations.jl")
end

ti = time() - ti
println("\nTest took total time of:")
println(round(ti, digits=3), " seconds or ", round(ti/60, digits=3), " minutes")


text:
struct Tabs
    options::SimpleList
    activetab::Observable{Int}
end

Tabs(options::SimpleList) = Tabs(options, Observable(1))
function Tabs(nt::Union{AbstractDict, NamedTuple}, activetab=Observable(1))
    options = Iterators.map(pairs(nt)) do (k, v)
        return SimpleDict("key" => k, "value" => v)
    end
    return Tabs(collect(Any, options), activetab)
end

function jsrender(session::Session, tabs::Tabs)
    options, activetab = tabs.options, tabs.activetab
    activeClasses = ["shadow", "bg-white"]
    inactiveClasses = String[]

    nodes = [DOM.li(
        class="text-blue-800 text-2xl font-semibold rounded mr-4 px-4 py-2 cursor-pointer hover:bg-gray-200",
        onclick=js"JSServe.update_obs($activetab, $i)",
        getkey(option)
    ) for (i, option) in enumerate(options)]
    headers = DOM.ul(class="flex mb-12", nodes)

    onjs(session, activetab, js"""
        function (idx) {
            $(UtilitiesJS).styleSelected($(nodes), idx - 1, $activeClasses, $inactiveClasses);
        }
    """)
    activetab[] = activetab[]

    contents = map(enumerate(options)) do (i, option)
        display = activetab[] == i ? "block" : "none"
        content = DOM.div(style="display: $display;", class="flex-auto", getvalue(option))
        onjs(
            session,
            activetab,
            js"""
                function (val) {
                    $(content).style.display = (val == $i) ? "block" : "none";
                }
            """
        )
        return jsrender(session, content)
    end
    return DOM.div(
        class="flex flex-col h-screen py-8",
        DOM.div(class="flex-initial", headers),
        contents...
    )
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule beefalo_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("beefalo")
JLLWrappers.@generate_main_file("beefalo", UUID("a30d6a8e-903d-54eb-aa61-387f835467ce"))
end  # module beefalo_jll


text:
type SAGsd <: StepData
	d
	y
	m::Int64
	numChunks::Int64
	batchSize::Int64
	getStep::Function
	stepString::AbstractString
	function SAGsd(numVars::Int64,y, stepString::AbstractString,numDp::Int64,numChunks::Int64, batchSize::Int64)
		d= sum(y)
		new(d, y, numDp, numChunks, batchSize, SAGComputation, stepString)
	end
end

function SAG(numVars::Int64, y, numDP::Int64,numChunks::Int64,batchSize::Int64)
	SAGsd(numVars,y, "SAG.$numChunks.$batchSize", numDP,numChunks,batchSize)
end

function SAGComputation(x, k, gnum, sd::SAGsd, problem::Problem; outputLevel =0)

	outputLevel >0 && println("  Starting SAG step computation")
	indices = sample(1:sd.numChunks,sd.batchSize;replace=false)
	outputLevel >1 && println("   indices = $indices")

	for i in indices

		outputLevel >1 && println("    Working on index $i")

		# step 1
		sd.d -= sd.y[i]
		sd.y[i] = zeros(sd.y[i] )

		# step 2
		for j in (i-1)*round(Int,floor(sd.m/sd.numChunks))+1:i*round(Int,floor(sd.m/sd.numChunks))
			outputLevel >1 && println("     Working on sample point $j")
			(f,sampleG) = (problem.getSampleFunctionAt(j))(x)
			sd.y[i]+=sampleG
			gnum += 1
		end

		#step 3
		sd.d += sd.y[i]

	end

	outputLevel >0 && println("  Finished SAG step computation")

	(sd.d/sd.m, gnum)
end


text:
# Example of vOptGeneric, for a 2IP (the 2ap03.mop is a 2LAP)
# the numerical instance is read on a file, according to the MOP format
# July 2017

using vOptGeneric
using GLPK ; using GLPKMathProgInterface

m = parseMOP( "2ap03.mop" , solver=GLPKSolverMIP() ) # comment: set up the path to the file if required
@time solve( m , method = :epsilon , step = 0.5 )

print_X_E(m)
getY_N(m)

max_stars_repo_name: faithyyyy/vOptSolver

max_issues_repo_name: faithyyyy/vOptSolver

max_forks_repo_name: faithyyyy/vOptSolver


text:
#=
   Logistic regression - Challenger

   See ~/jags/logistic_regression_challenger.R
   """
   model {
     for (i in 1:N){
       y[i] ~ dbern(p[i])
       p[i] <- 1 / (1 + exp(-z[i]))
       z[i] <- w0 + w1 * x[i]
     }
     w0 ~ dnorm(0, .001)
     w1 ~ dnorm(0, .0001)

   Output:
      Mean     SD  Naive SE Time-series SE
w0 16.8900 7.7265 0.0122167       0.290180
w1 -0.2602 0.1135 0.0001794       0.004232

2. Quantiles for each variable:

     2.5%     25%     50%     75%    97.5%
w0  4.131 11.3621 16.1240 21.5029 34.27954
w1 -0.516 -0.3279 -0.2488 -0.1789 -0.07332
   """

  From https://www.zinkov.com/posts/2012-06-27-why-prob-programming-matters/
  "Logistic Regression"
  """
  Logistic Regression can be seen as a generalization of Linear Regression where the output is
  transformed
  to lie between 0 and 1. This model only differs from the previous one by a single line, illustrating that
  adding this complexity does not require starting from scratch. The point with probabilistic programming
  is you are able to explore slightly more complex models very easily.
  """

  From https://www.stat.ubc.ca/~bouchard/courses/stat520-sp2014-15/lecture/2015/02/27/notes-lecture3.html
  x = 66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81,76,79,75,76,58
  y = 1,0,1,1,1,1,1,1,0,0,0,1,1,0,1,1,1,1,1,1,0,1,0

  Compared to logistic_regression_challenger.jl here we use the logicstic function instead.

    Summary Statistics
    parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
        Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

            w0   -9.1270    5.3364     0.1688    1.1745    16.5012    1.0096        2.1562
            w1    0.1893    0.4241     0.0134    0.0245   167.3714    1.0025       21.8700
    y_post[1]    0.7570    0.4291     0.0136    0.0251   578.5337    1.0030       75.5957
    y_post[2]    0.8380    0.3686     0.0117    0.0158   653.2435    0.9990       85.3578
    y_post[3]    0.8240    0.3810     0.0120    0.0177   548.3477    0.9990       71.6513
    y_post[4]    0.8140    0.3893     0.0123    0.0171   519.4787    1.0042       67.8791
    y_post[5]    0.7880    0.4089     0.0129    0.0174   827.3196    1.0005      108.1040
    y_post[6]    0.8680    0.3387     0.0107    0.0117   499.8656    0.9991       65.3163
    y_post[7]    0.9100    0.2863     0.0091    0.0125   418.0910    1.0039       54.6310
    y_post[8]    0.8440    0.3630     0.0115    0.0203   380.3478    0.9995       49.6992
        ⋮           ⋮         ⋮         ⋮          ⋮         ⋮          ⋮           ⋮
                                                                            15 rows omitted

    Posterior predictive:
    1: true: 1 pred: 1 (0.757)
    2: true: 0 pred: 1 (0.838)!
    3: true: 1 pred: 1 (0.824)
    4: true: 1 pred: 1 (0.814)
    5: true: 1 pred: 1 (0.788)
    6: true: 1 pred: 1 (0.868)
    7: true: 1 pred: 1 (0.91)
    8: true: 1 pred: 1 (0.844)
    9: true: 0 pred: 0 (0.415)
    10: true: 0 pred: 1 (0.673)!
    11: true: 0 pred: 1 (0.83)!
    12: true: 1 pred: 1 (0.935)
    13: true: 1 pred: 1 (0.77)
    14: true: 0 pred: 0 (0.297)
    15: true: 1 pred: 1 (0.767)
    16: true: 1 pred: 1 (0.911)
    17: true: 1 pred: 1 (0.853)
    18: true: 1 pred: 1 (0.952)
    19: true: 1 pred: 1 (0.922)
    20: true: 1 pred: 1 (0.949)
    21: true: 0 pred: 1 (0.92)!
    22: true: 1 pred: 1 (0.923)
    23: true: 0 pred: 0 (0.437)
    num_correct: 19 / 23 (82.6086956521739%)

  Cf logistic_regression_challenger.jl
      ~/blog/logistic_regression-challenger.blog
     ~/webppl/logistic_regression-challenger.wppl

=#

using Turing, StatsPlots

include("jl_utils.jl")

@model function logistic_regression_challenger2(xs,ys)
    n = length(xs)
    w0 ~ Normal(0, 10)
    w1 ~ Normal(0, 10)

    p = tzeros(n)
    for i in 1:n
        p[i] = logistic(w0 + w1 * xs[i])
        ys[i] ~ Bernoulli(p[i])
    end

    y_post = tzeros(n)
    for i in 1:n
        y_post[i] ~ Bernoulli(p[i])
    end
end

xs = [66,70,69,68,67,72,73,70,57,63,70,78,67,53,67,75,70,81,76,79,75,76,58];
# ys = [true,false,true,true,true,true,true,true,false,false,false,true,true,false,true,true,true,true,true,true,false,true,false];
ys = [1,0,1,1,1,1,1,1,0,0,0,1,1,0,1,1,1,1,1,1,0,1,0];
#     1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3

model = logistic_regression_challenger2(xs,ys)

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 1_000)
# chns = sample(model, SMC(), 10_000)
# chns = sample(model, IS(), 10_000)

display(chns)
# display(plot(chns))

# display(plot(chns[[:w0,:w1]]) # quite bad!

println("Posterior predictive:")
num_correct = 0
n = length(ys)
for i in 1:n
    pred1 = mean(chns["y_post[$i]"])
    pred2 = round(Int64,pred1)
    c = "!"
    if pred2 == ys[i]
        global num_correct += 1
        c = ""
    end
    println("$i: true: $(ys[i]) pred: $pred2 ($pred1)$c" )
end
println("num_correct: $num_correct / $n ($(100*num_correct/n)%)")

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
export ProdAnalytical

import Base: prod

"""
    ProdAnalytical

`ProdAnalytical` is one of the strategies for `prod` function. This strategy uses analytical prod methods but does not constraint a prod to be in any specific form.
It fails if no analytical rules is available, use `ProdGeneric` prod strategy to fallback to approximation methods.

See also: [`prod`](@ref), [`ProdPreserveType`](@ref), [`ProdGeneric`](@ref)
"""
struct ProdAnalytical end

"""
    prod(strategy, left, right)

`prod` function is used to find a product of two probability distrubution over same variable (e.g. 𝓝(x|μ_1, σ_1) × 𝓝(x|μ_2, σ_2)).
There are multiple strategies for prod function, e.g. `ProdAnalytical`, `ProdGeneric` or `ProdPreserveType`.

# Examples:
```jldoctest
using ReactiveMP

product = prod(ProdAnalytical(), NormalMeanVariance(-1.0, 1.0), NormalMeanVariance(1.0, 1.0))

mean(product), var(product)

# output
(0.0, 0.5)
```

See also: [`prod_analytical_rule`](@ref), [`ProdAnalytical`](@ref), [`ProdGeneric`](@ref)
"""
prod(::ProdAnalytical, left, right) = error("No analytical rule available to compute a product of distributions $(left) and $(right).")

prod(::ProdAnalytical, ::Missing, right)     = right
prod(::ProdAnalytical, left, ::Missing)      = left
prod(::ProdAnalytical, ::Missing, ::Missing) = missing