text:
using SignatureGB
using Singular
P, (x1, x2, x3, x4, x5, x6, x7, x8, x9)  = Singular.PolynomialRing(Fp(65521), ["x1", "x2", "x3", "x4", "x5", "x6", "x7", "x8", "x9"]);
gb = SignatureGB.f5([x1, x2, x3, x4, x5, x6, x7, x8, x9]);F = [3844*x1^4+841*x1^2*x2^2+992*x1^2*x2*x6+2916*x1^2*x5^2+10044*x1^2*x5*x9+289*x1^2*x6^2+3968*x1^2*x6*x7-722*x1*x2^2*x6+3886*x1*x2*x3*x5+3596*x1*x2*x4*x8-4104*x1*x2*x5^2-4640*x1*x2*x5*x8+490*x1*x2*x5*x9-2924*x1*x2*x6*x7-6696*x1*x5*x6*x7+3564*x1*x5*x7^2+2268*x1*x5*x7*x8+782*x1*x6^2*x9+1836*x1*x6*x7*x8+324*x2^4-3096*x2^3*x7+1768*x2^2*x5^2-2508*x2^2*x5*x9+593*x2^2*x6^2+828*x2^2*x6*x9+7396*x2^2*x7^2+1944*x2^2*x7*x8+1089*x2^2*x9^2-3082*x2*x3*x5*x6-3528*x2*x3*x5*x8-2852*x2*x4*x6*x8+4712*x2*x5*x6*x7+3680*x2*x5*x6*x8+3734*x2*x5*x6*x9-2436*x2*x5*x7^2-1596*x2*x5*x7*x8-3492*x2*x5*x8^2+936*x2*x5*x8*x9+512*x2*x6^2*x7-8048*x2*x6*x7*x9-9288*x2*x7^2*x8+2178*x2*x7^2*x9+1386*x2*x7*x8*x9+4489*x3^2*x5^2+6561*x3^2*x7^2-4050*x3^2*x7*x9+9604*x3^2*x8^2+629*x3^2*x9^2+8308*x3*x4*x5*x8+88*x3*x4*x6*x9-10720*x3*x5^2*x8-7102*x3*x5^2*x9-316*x3*x6^2*x9+136*x3*x6*x7*x9-296*x3*x6*x8*x9-12312*x3*x7^3-392*x3*x7^2*x8+3800*x3*x7^2*x9+19012*x3*x8^3-5112*x3*x8^2*x9+484*x4^2*x6^2+3844*x4^2*x8^2-9920*x4*x5*x8^2-6572*x4*x5*x8*x9-3476*x4*x6^3+1496*x4*x6^2*x7-3256*x4*x6^2*x8-176*x4*x6*x8^2+6400*x5^2*x8^2+8480*x5^2*x8*x9+9370*x5^2*x9^2+5184*x5*x6*x7*x9+6241*x6^4-5372*x6^3*x7+11692*x6^3*x8+6024*x6^2*x7^2-5032*x6^2*x7*x8+6108*x6^2*x8^2+529*x6^2*x9^2-4092*x6*x7^3-2604*x6*x7^2*x8-272*x6*x7*x8^2+2484*x6*x7*x8*x9+592*x6*x8^3+6869*x7^4+1386*x7^3*x8+2969*x7^2*x8^2+104*x7^2*x8*x9+9425*x8^4-5044*x8^3*x9+676*x8^2*x9^2+2108*x1^2*x8-2520*x1*x2*x5-2040*x1*x2*x6+13720*x1*x3*x8-280*x1*x7^2+13580*x1*x8^2-3640*x1*x8*x9-2160*x2^3+10320*x2^2*x7+272*x2*x6*x8-2760*x2*x6*x9-6480*x2*x7*x8+9072*x3*x4*x7-2800*x3*x4*x9-9234*x3*x7^2-4536*x3*x7*x8+2850*x3*x7*x9+1400*x3*x8*x9-8512*x4*x7^2+2754*x5*x8*x9+1088*x6*x7*x8+8664*x7^3+4256*x7^2*x8+16308*x1^2+3600*x2^2+1472*x2*x6+3136*x4^2-6384*x4*x7-3136*x4*x8+14904*x5*x9+5888*x6*x7+3249*x7^2+3192*x7*x8+1073*x8^2+3128*x8+8464, 1682*x1^2*x2+992*x1^2*x6-1444*x1*x2*x6+3886*x1*x3*x5+3596*x1*x4*x8-4104*x1*x5^2-4640*x1*x5*x8+490*x1*x5*x9-2924*x1*x6*x7+1296*x2^3-9288*x2^2*x7+3536*x2*x5^2-5016*x2*x5*x9+1186*x2*x6^2+1656*x2*x6*x9+14792*x2*x7^2+3888*x2*x7*x8+2178*x2*x9^2-3082*x3*x5*x6-3528*x3*x5*x8-2852*x4*x6*x8+4712*x5*x6*x7+3680*x5*x6*x8+3734*x5*x6*x9-2436*x5*x7^2-1596*x5*x7*x8-3492*x5*x8^2+936*x5*x8*x9+512*x6^2*x7-8048*x6*x7*x9-9288*x7^2*x8+2178*x7^2*x9+1386*x7*x8*x9-2520*x1*x5-2040*x1*x6-6480*x2^2+20640*x2*x7+272*x6*x8-2760*x6*x9-6480*x7*x8+7200*x2+1472*x6, 3886*x1*x2*x5-3082*x2*x5*x6-3528*x2*x5*x8+8978*x3*x5^2+13122*x3*x7^2-8100*x3*x7*x9+19208*x3*x8^2+1258*x3*x9^2+8308*x4*x5*x8+88*x4*x6*x9-10720*x5^2*x8-7102*x5^2*x9-316*x6^2*x9+136*x6*x7*x9-296*x6*x8*x9-12312*x7^3-392*x7^2*x8+3800*x7^2*x9+19012*x8^3-5112*x8^2*x9+13720*x1*x8+9072*x4*x7-2800*x4*x9-9234*x7^2-4536*x7*x8+2850*x7*x9+1400*x8*x9, 3596*x1*x2*x8-2852*x2*x6*x8+8308*x3*x5*x8+88*x3*x6*x9+968*x4*x6^2+7688*x4*x8^2-9920*x5*x8^2-6572*x5*x8*x9-3476*x6^3+1496*x6^2*x7-3256*x6^2*x8-176*x6*x8^2+9072*x3*x7-2800*x3*x9-8512*x7^2+6272*x4-6384*x7-3136*x8, 5832*x1^2*x5+10044*x1^2*x9+3886*x1*x2*x3-8208*x1*x2*x5-4640*x1*x2*x8+490*x1*x2*x9-6696*x1*x6*x7+3564*x1*x7^2+2268*x1*x7*x8+3536*x2^2*x5-2508*x2^2*x9-3082*x2*x3*x6-3528*x2*x3*x8+4712*x2*x6*x7+3680*x2*x6*x8+3734*x2*x6*x9-2436*x2*x7^2-1596*x2*x7*x8-3492*x2*x8^2+936*x2*x8*x9+8978*x3^2*x5+8308*x3*x4*x8-21440*x3*x5*x8-14204*x3*x5*x9-9920*x4*x8^2-6572*x4*x8*x9+12800*x5*x8^2+16960*x5*x8*x9+18740*x5*x9^2+5184*x6*x7*x9-2520*x1*x2+2754*x8*x9+14904*x9, 992*x1^2*x2+578*x1^2*x6+3968*x1^2*x7-722*x1*x2^2-2924*x1*x2*x7-6696*x1*x5*x7+1564*x1*x6*x9+1836*x1*x7*x8+1186*x2^2*x6+828*x2^2*x9-3082*x2*x3*x5-2852*x2*x4*x8+4712*x2*x5*x7+3680*x2*x5*x8+3734*x2*x5*x9+1024*x2*x6*x7-8048*x2*x7*x9+88*x3*x4*x9-632*x3*x6*x9+136*x3*x7*x9-296*x3*x8*x9+968*x4^2*x6-10428*x4*x6^2+2992*x4*x6*x7-6512*x4*x6*x8-176*x4*x8^2+5184*x5*x7*x9+24964*x6^3-16116*x6^2*x7+35076*x6^2*x8+12048*x6*x7^2-10064*x6*x7*x8+12216*x6*x8^2+1058*x6*x9^2-4092*x7^3-2604*x7^2*x8-272*x7*x8^2+2484*x7*x8*x9+592*x8^3-2040*x1*x2+272*x2*x8-2760*x2*x9+1088*x7*x8+1472*x2+5888*x7, 3968*x1^2*x6-2924*x1*x2*x6-6696*x1*x5*x6+7128*x1*x5*x7+2268*x1*x5*x8+1836*x1*x6*x8-3096*x2^3+14792*x2^2*x7+1944*x2^2*x8+4712*x2*x5*x6-4872*x2*x5*x7-1596*x2*x5*x8+512*x2*x6^2-8048*x2*x6*x9-18576*x2*x7*x8+4356*x2*x7*x9+1386*x2*x8*x9+13122*x3^2*x7-4050*x3^2*x9+136*x3*x6*x9-36936*x3*x7^2-784*x3*x7*x8+7600*x3*x7*x9+1496*x4*x6^2+5184*x5*x6*x9-5372*x6^3+12048*x6^2*x7-5032*x6^2*x8-12276*x6*x7^2-5208*x6*x7*x8-272*x6*x8^2+2484*x6*x8*x9+27476*x7^3+4158*x7^2*x8+5938*x7*x8^2+208*x7*x8*x9-560*x1*x7+10320*x2^2-6480*x2*x8+9072*x3*x4-18468*x3*x7-4536*x3*x8+2850*x3*x9-17024*x4*x7+1088*x6*x8+25992*x7^2+8512*x7*x8-6384*x4+5888*x6+6498*x7+3192*x8, 3596*x1*x2*x4-4640*x1*x2*x5+2268*x1*x5*x7+1836*x1*x6*x7+1944*x2^2*x7-3528*x2*x3*x5-2852*x2*x4*x6+3680*x2*x5*x6-1596*x2*x5*x7-6984*x2*x5*x8+936*x2*x5*x9-9288*x2*x7^2+1386*x2*x7*x9+19208*x3^2*x8+8308*x3*x4*x5-10720*x3*x5^2-296*x3*x6*x9-392*x3*x7^2+57036*x3*x8^2-10224*x3*x8*x9+7688*x4^2*x8-19840*x4*x5*x8-6572*x4*x5*x9-3256*x4*x6^2-352*x4*x6*x8+12800*x5^2*x8+8480*x5^2*x9+11692*x6^3-5032*x6^2*x7+12216*x6^2*x8-2604*x6*x7^2-544*x6*x7*x8+2484*x6*x7*x9+1776*x6*x8^2+1386*x7^3+5938*x7^2*x8+104*x7^2*x9+37700*x8^3-15132*x8^2*x9+1352*x8*x9^2+2108*x1^2+13720*x1*x3+27160*x1*x8-3640*x1*x9+272*x2*x6-6480*x2*x7-4536*x3*x7+1400*x3*x9+2754*x5*x9+1088*x6*x7+4256*x7^2-3136*x4+3192*x7+2146*x8+3128, 10044*x1^2*x5+490*x1*x2*x5+782*x1*x6^2-2508*x2^2*x5+828*x2^2*x6+2178*x2^2*x9+3734*x2*x5*x6+936*x2*x5*x8-8048*x2*x6*x7+2178*x2*x7^2+1386*x2*x7*x8-4050*x3^2*x7+1258*x3^2*x9+88*x3*x4*x6-7102*x3*x5^2-316*x3*x6^2+136*x3*x6*x7-296*x3*x6*x8+3800*x3*x7^2-5112*x3*x8^2-6572*x4*x5*x8+8480*x5^2*x8+18740*x5^2*x9+5184*x5*x6*x7+1058*x6^2*x9+2484*x6*x7*x8+104*x7^2*x8-5044*x8^3+1352*x8^2*x9-3640*x1*x8-2760*x2*x6-2800*x3*x4+2850*x3*x7+1400*x3*x8+2754*x5*x8+14904*x5];
sort!(F, by = p -> leading_monomial(p))
println(@elapsed SignatureGB.f5(F, verbose=true));exit()


text:
module CSV2Parquet

using DataConvenience

Base.@ccallable function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    for file in ARGS
        if !isfile(file)
            error("could not find file $file")
        end

        chunks = DataConvenience.CsvChunkIterator(file)
        for (n, chunk) in enumerate(chunks)
            println("chunk #: ", n)
            println("chunk size: ", size(chunk))
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end # module

ext: jl

max_stars_repo_name: xiaodaigh/csv-to-parquet

max_issues_repo_name: xiaodaigh/csv-to-parquet

max_forks_repo_name: xiaodaigh/csv-to-parquet


text:
using LinearAlgebra: I

@testset "R" begin
    @testset "radrec" begin
        @test [1.0, 0.0, 0.0] ≈ radrec(1.0, 0.0, 0.0)
        @test [0.0, 1.0, 0.0] ≈ radrec(1.0, π/2, 0.0)
        @test [0.0, 0.0, 1.0] ≈ radrec(1.0, 0.0, π/2)
    end
    @testset "rav2xf" begin
        e = [1.0, 0.0, 0.0]
        rz = [0.0 1.0 0.0; -1.0 0.0 0.0; 0.0 0.0 1.0]
        exp = [0.0  1.0  0.0  0.0  0.0  0.0;
               -1.0  0.0  0.0  0.0  0.0  0.0;
               0.0  0.0  1.0  0.0  0.0  0.0;
               0.0  0.0  1.0  0.0  1.0  0.0;
               0.0  0.0  0.0 -1.0  0.0  0.0;
               0.0 -1.0  0.0  0.0  0.0  1.0]
        act = rav2xf(rz, e)
        @testset for i in eachindex(exp, act)
            @test exp[i] ≈ act[i]
        end
    end
    @testset "raxisa" begin
        axis = [1.0, 2.0, 3.0]
        angle = 0.2π
        rotate_matrix = axisar(axis, angle)
        axout, angout = raxisa(rotate_matrix)
        expected_angout = [0.26726124, 0.53452248, 0.80178373]
        @test angout ≈ 0.62831853
        @test axout ≈ expected_angout
        @test_throws SpiceError raxisa(randn(3, 3))
        @test_throws ArgumentError raxisa(randn(2, 2))
    end
    @testset "reccyl" begin
        expected1 = (0.0, 0.0, 0.0)
        expected2 = (1.0, π/2, 0.0)
        expected3 = (1.0, 3/2 * π, 0.0)
        @test all(reccyl([0.0, 0.0, 0.0]) .≈ expected1)
        @test all(reccyl([0.0, 1.0, 0.0]) .≈ expected2)
        @test all(reccyl([0.0, -1.0, 0.0]) .≈ expected3)
        @test_throws ArgumentError reccyl([0.0, -1.0])
    end
    @testset "recgeo" begin
        try
            furnsh(path(CORE, :pck))
            radii = bodvrd("EARTH", "RADII", 3)
            flat = (radii[1] - radii[3]) / radii[1]
            x = [-2541.748162, 4780.333036, 3360.428190]
            lon, lat, alt = recgeo(x, radii[1], flat)
            actual = [rad2deg(lon), rad2deg(lat), alt]
            expected = [118.0, 32.0, 0.001915518]
            @testset for i in eachindex(actual, expected)
                @test actual[i] ≈ expected[i] rtol=1e-4
            end
            @test_throws ArgumentError recgeo(x[1:2], radii[1], flat)
            @test_throws SpiceError recgeo(x, -radii[1], flat)
            @test_throws SpiceError recgeo(x, radii[1], 1.0)
        finally
            kclear()
        end
    end
    @testset "reclat" begin
        act1 = reclat([1.0, 0.0, 0.0])
        act2 = reclat([0.0, 1.0, 0.0])
        act3 = reclat((-1.0, 0.0, 0.0))
        @test [act1[1], act1[2], act1[3]] ≈ [1.0, 0.0, 0.0]
        @test [act2[1], act2[2], act2[3]] ≈ [1.0, deg2rad(90.0), 0.0]
        @test [act3[1], act3[2], act3[3]] ≈ [1.0, deg2rad(180.0), 0.0]
        @test_throws ArgumentError reclat([1.0, 2.0])
    end
    @testset "recpgr" begin
        try
            furnsh(path(CORE, :lsk), path(CORE, :pck), path(CORE, :spk))
            radii = bodvrd("MARS", "RADII", 3)
            flat = (radii[1] - radii[3]) / radii[1]
            x = [0.0, -2620.678914818178, 2592.408908856967]
            lon, lat, alt = recpgr("MARS", x, radii[1], flat)
            actual = [rad2deg(lon), rad2deg(lat), alt]
            expected = [90., 45, 300]
            @test actual ≈ expected
            @test_throws ArgumentError recpgr("MARS", x[1:2], radii[1], flat)
            @test_throws SpiceError recpgr("MARS", x, -radii[1], flat)
            @test_throws SpiceError recpgr("MARS", x, radii[1], 1.0)
        finally
            kclear()
        end
    end
    @testset "recrad" begin
        act1 = collect(recrad([1.0, 0.0, 0.0]))
        act2 = collect(recrad([0.0, 1.0, 0.0]))
        act3 = collect(recrad([0.0, 0.0, 1.0]))
        exp1 = [1.0, 0.0, 0.0]
        @testset for i in eachindex(act1, exp1)
            @test act1[i] ≈ exp1[i]
        end
        exp2 = [1.0, deg2rad(90), 0.0]
        @testset for i in eachindex(act2, exp2)
            @test act2[i] ≈ exp2[i]
        end
        exp3 = [1.0, 0.0, deg2rad(90)]
        @testset for i in eachindex(act3, exp3)
            @test act3[i] ≈ exp3[i]
        end
    end
    @testset "recsph" begin
        v1 = [-1.0, 0.0, 0.0]
        @test all(recsph(v1) .≈ (1.0, π/2, π))
    end
    @testset "removc" begin
        cell = SpiceCharCell(10, 10)
        items = ["one", "two", "three", "four"]
        for i in items
            insrtc!(cell, i)
        end
        remove_items = ["three", "four"]
        for r in remove_items
            removc!(cell, r)
        end
        expected = ["one", "two"]
        @test cell == expected
    end
    @testset "removd" begin
        cell = SpiceDoubleCell(10)
        items = [0.0, 1.0, 1.0, 2.0, 3.0, 5.0, 8.0, 13.0, 21.0]
        for i in items
            insrtd!(cell, i)
        end
        remove_items = [0.0, 2.0, 4.0, 6.0, 8.0, 12.0]
        for r in remove_items
            removd!(cell, r)
        end
        expected = [1.0, 3.0, 5.0, 13.0, 21.0]
        @testset for i in eachindex(cell, expected)
            @test cell[i] == expected[i]
        end
    end
    @testset "removi" begin
        cell = SpiceIntCell(10)
        items = [0, 1, 1, 2, 3, 5, 8, 13, 21]
        for i in items
            insrti!(cell, i)
        end
        remove_items = [0, 2, 4, 6, 8, 12]
        for r in remove_items
            removi!(cell, r)
        end
        expected = [1, 3, 5, 13, 21]
        @testset for i in eachindex(cell, expected)
            @test cell[i] == expected[i]
        end
    end
    @testset "reordc" begin
        array = ["one", "three", "two", "zero"]
        iorder = [4, 1, 3, 2]
        outarray = SPICE._reordc(iorder, copy(array))
        @test_skip outarray == array[iorder]
    end
    @testset "reordd" begin
        array = [1.0, 3.0, 2.0]
        iorder = [1, 3, 2]
        outarray = SPICE._reordd(iorder, copy(array))
        @test outarray == array[iorder]
    end
    @testset "reordi" begin
        array = [1, 3, 2]
        iorder = [1, 3, 2]
        outarray = SPICE._reordi(iorder, copy(array))
        @test outarray == array[iorder]
    end
    @testset "reordl" begin
        array = [true, true, false]
        iorder = [1, 3, 2]
        outarray = SPICE._reordl(iorder, copy(array))
        @test outarray == array[iorder]
    end
    @testset "repmc" begin
        stringtestone = "The truth is #"
        outstringone = SPICE._repmc(stringtestone, "#", "SPICE")
        @test replace(stringtestone, "#" => "SPICE") == outstringone
    end
    @testset "rotate" begin
        act = rotate(π / 4, 3)
        exp = [sqrt(2) / 2.0 sqrt(2) / 2.0 0.0;
              -sqrt(2) / 2.0 sqrt(2) / 2.0 0.0;
               0.0 0.0 1.0]
        @testset for i in eachindex(act, exp)
             @test act[i] ≈ exp[i]
        end
    end
    @testset "rotmat" begin
        ident = Array{Float64}(I, 3, 3)
        expected_r = [0.0 0.0 -1.0;
                      0.0 1.0 0.0;
                      1.0 0.0 0.0]
        r_out = rotmat(ident, π/2, 2)
        @test r_out ≈ expected_r
        @test_throws ArgumentError rotmat(ident[1:2, 1:2], 0.0, 3)
    end
    @testset "rotvec" begin
        vin = [sqrt(2), 0.0, 0.0]
        angle = π / 4
        iaxis = 3
        v_expected = [1.0, -1.0, 0.0]
        vout = rotvec(vin, angle, iaxis)
        @test vout ≈ v_expected
        @test_throws ArgumentError rotvec(vin[1:2], angle, iaxis)
    end
    @testset "rpd" begin
        @test deg2rad(1.0) == SPICE._rpd()
    end
    @testset "rquad" begin
        # solve x^2 + 2x + 3 = 0
        root1, root2 = rquad(1.0, 2.0, 3.0)
        expected_root1 = [-1.0, sqrt(2.0)]
        expected_root2 = [-1.0, -sqrt(2.0)]
        @test root1 == expected_root1
        @test root2 == expected_root2
        # This test prints a funny CSPICE error message
        # @test_throws SpiceError rquad(0.0, 0.0, 3.0)
    end
end


text:
using Parameters
using Blink
using Interact
using Plots
gr()

include("lv_model.jl")

function iso1(u2, p)
    @unpack α12, K1 = p
    return K1 - α12 * u2 / K1
end

function iso2(u2, p)
    @unpack α21, K2 = p
    return (K2 - u2) / α21
end

let
    u2s = range(0, 5, length = 100)
    ui = @manipulate for α12 in slider(0.1:0.01:1.0; label = "α12", value = 0.8),
                         α21 in slider(0.1:0.01:1.0; label = "α21", value = 0.3),
                         K1 in slider(0.1:0.01:2.0; label = "K1", value = 1.5),
                         K2 in slider(0.1:0.01:2.0; label = "K2", value = 1.0)
        p = LVPar(α12 = α12, α21 = α21, K1 = K1, K2 = K2)

        plot(u2s, [iso1(u2, p) for u2 in u2s], label = "u1 = 0")
        plot!(u2s, [iso2(u2, p) for u2 in u2s], label = "u2 = 0")
        xlims!(0, 3)
        ylims!(0, 5)
        xlabel!("u2")
        ylabel!("u1")
    end

    w = Blink.AtomShell.Window()
    body!(w, ui)
end

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
using Pkg
Pkg.status()
using Plots
using DifferentialEquations

#----------------------------------------------
# Set Parameters:
# Membrane capacitance and potential
C_m = 1.0
V_m = -65.0

# Max sodium conductance and membrane potential
g_Na = 120.0
V_Na = 56.0

# Max potassium conductance and membrane potential
g_K = 36.0
V_K = -77.0

# Max leak conductance and membrane potential
g_l = 0.3
V_l = -60.0

# Time to integrate over
time = 0.0:0.1:1000

# ion variables
m = 0.05
h = 0.6
n = 0.32

#----------------------------------------------
# Potassium Rate Functoins
a_n(V_m) =  0.01 * (V_m + 55) / (1 - exp(-(V_m + 55.0) / 10))

B_n(V_m) = 0.125 * exp(-(V_m + 65.0) / 80)

#----------------------------------------------
# Sodium Rate Functions
a_m(V_m) = 0.1 * (40 + V_m) / (1.0 - exp(-(V_m + 40.0) / 10.0))

a_h(V_m) = 0.07 * exp(-(V_m + 65.0) / 20.0)

B_m(V_m) = 4 * exp(-(V_m + 65.0) / 18.0)

B_h(V_m) = 1 / (exp(-(35 + V_m) / 10) + 1)

#----------------------------------------------
# Ion Currents
IK(g_K, n, V_m, V_K) = g_K * n^(4) * (V_m - V_K)

INa(g_Na, m, h, V_m, V_Na) =  g_Na * m^(3) * h * (V_m - V_Na)

LI(g_l, V_m, V_l) = g_l * (V_m - V_l)

#----------------------------------------------
# Input current
I_Inj = 20

input(t::AbstractFloat) = ifelse(0 < t < 200, zero(t), ifelse(200<t<300, one(t)*I_Inj, zero(t)))

#----------------------------------------------
# Hodking Huxley model
function hodgkin_huxley!(du, u, p, t)
    # dVdt, dndt, dmdt, dhdt = u
    #
    #
    V_m, n, m, h = u
    g_K, V_K, g_Na, V_Na, g_l, C_m = p

    # Total current through the membrane
    du[1] = input(t) - IK(g_K,n,V_m,V_K) - INa(g_Na, m, h, V_m, V_Na) - LI(g_l, V_m, V_l) / C_m

    # Derivative of n, potassium channel activation, w.r.t. time
    du[2] = a_n(V_m) * (1 - n) - B_n(V_m) * n

    # Derivative of m, sodium channel activion, w.r.t. time
    du[3] = a_m(V_m) * (1 - m) - B_m(V_m) * m

    # Derivative of h, leaky channel in-activion, w.r.t. time

    du[4] = a_h(V_m) * (1 - h) - B_h(V_m) * h

end

#----------------------------------------------
# Define and solve ODE problem
p = [g_K, V_K, g_Na, V_Na, g_l, C_m]
u0 = [V_m, n, m, h]
tspan = (0.0,500.0)

prob = ODEProblem(hodgkin_huxley!,u0,tspan,p)
sol = solve(prob)

#----------------------------------------------
# Plotting Solutions
plot(sol, vars=1)
#animate(sol, vars=(1), every=20)

plot(sol, vars=(2,3))

plot(sol, vars=(2,3,4))


text:
"""
Tools for the exploration of chaos and nonlinear dynamics
"""
module ChaosTools

using Reexport
@reexport using DynamicalSystemsBase

using DynamicalSystemsBase: DS, DDS, CDS
using DynamicalSystemsBase: MDI, TDI
using DynamicalSystemsBase: stateeltype
using OrdinaryDiffEq
import OrdinaryDiffEq: ODEIntegrator

# Lyapunovs:
include("lyapunovs.jl")

# Phase space related:
include("orbitdiagrams/discrete_diagram.jl")
include("orbitdiagrams/poincare.jl")

# Entropies and Dimension Estimation:
include("dimensions/entropies.jl")
include("dimensions/dims.jl")

# Nonlinear Timeseries Analysis:
include("nlts.jl")
include("estimate_reconstruction.jl")

# Periodicity:
include("periodic.jl")

# Chaos Detection:
include("chaos_detection.jl")

# Ugly methods that shouldn't exist:
include("ugliness.jl")

end # module


text:
using Test

using iLQGames:
    LinearSystem,
    discretize,
    linearize,
    linearize_discrete,
    LinearizationStyle,
    TrivialLinearization

using StaticArrays

A = @SMatrix [1 2; 3 4]
B = @SMatrix [1; 3]
ls = LinearSystem{0}(A, B)
ΔT = 0.1241241

# some test point
x0 = @SVector [1., 2.]
u0 = @SVector [2.]
t0 = .11

d1 = discretize(ls, Val{ΔT}())

@testset "linearization style" begin
    @test LinearizationStyle(ls) isa TrivialLinearization
    @test LinearizationStyle(d1) isa TrivialLinearization
end

@testset "linearize" begin
    # linearizing a linear system should not make a difference
    ls_lin = linearize(ls, x0, u0, t0)
    @test ls.A == ls_lin.A
    @test ls.B == ls_lin.B
end;


text:
# To test implementation of parser, and understand the working structure
using JLD, Knet
include("../parser/types.jl")
include("../parser/preprocess.jl")
include("../parser/parser.jl")
include("../parser/helper.jl")
include("../parser/modelutils.jl")
include("set_testenv.jl")

function test_parser(;mode=:random)
    p1, s1 = get_parser_sentence()

end

function compare_moves(p1, s1)
    p_r = deepcopy(p1);
    old_acs = nothing
    for i in 1:500
        p_r = deepcopy(p1)

        acs = done_randmoves(p_r, s1)
        (acs == old_acs) && prinlnt("Previous and currrent actions are the same")
        old_acs = acs
        if p_r.head == s1.head
            println("Correct sequence is catch")
        end
    end
end

# extract the gold moves from given parser
function done_goldmoves(porig, s1)
    actions = []
    p1 = deepcopy(porig)
    totmoves = 0
    parserdone = false
    while !parserdone
        mcosts = movecosts(p1, p1.sentence.head, p1.sentence.deprel)
        goldmove = indmin(mcosts)
        if mcosts[goldmove] == typemax(Cost)
            parserdone = true
        else
            push!(actions, goldmove)
            move!(p1, goldmove)
            totmoves += 1;
        end
    end
    println("GoldMoves: Total moves $totmoves total words $(length(s1.word))")
    return p1, actions
end

# Done random moves
function done_randmoves(pin, s1)
    p1 = deepcopy(pin)
    ac_seq = []
    totmoves = 0
    mcosts = movecosts(p1, p1.sentence.head, p1.sentence.deprel)
    n = length(mcosts)
    while anyvalidmoves(p1)
        mcosts = movecosts(p1, p1.sentence.head, p1.sentence.deprel)
        validmoves = find(mcosts .< typemax(Cost))
        m = validmoves[rand(1:length(validmoves))]
        push!(ac_seq, m)
        totmoves += 1;
        move!(p1, m)
    end
    println("RandomMoves: Total moves $totmoves total words $(length(s1.word))")
    return ac_seq
end

function get_parser_sentence()
    d = load(chmodel); vocab = create_vocab(d)
    corpora = []
    for f in [real_tdata3, real_ddata3]
        c = loadcorpus(f, vocab)
        push!(corpora, c)
    end
    sentence = corpora[1][4]
    p1 = ArcHybridR1(sentence)
    return p1, sentence
end

# LEFT-ARC : (σ|wi,wj|β,A) ⇒ (σ,wj|β,A∪{(wj,r,wi)})

# To perform Left-Arc operation you need two preconditions:
# - Stack and buffer are non-empty and wi ≠ ROOT

# RIGHT-ARC : (σ|wi , wj|β, A) ⇒ (σ, wi|β, A∪{(wi , r, wj)})
# To perform right-arc operation the only precondition is that the stack and buffer are non-empty.

# we have total of p.moves that are represented as :
# shiftmove has a value of 1
# leftmoves are the even numbers from 2 to p.nmove-2
# rightmoves are the odd numbers from 3 to p.nmove-1

# Here are some other notes on that:
# nword stands for number of words in that sentence, ndeps stands for number of dependency labels in that sentence

# init! function modifies the number of possible moves in that parser, and then applies the shift operation to the parser

# In shift move:
# - increase p.sptr by 1,
# - make p.stack[p.sptr] = p.wptr -> that puts ith word on buffer to stack.
# - increase the p.wptr by 1

# SHIFT MOVE implementation
# shift(p::Parser)=(p.stack[p.sptr+=1]=p.wptr; p.wptr+=1)

# Parser fields:
# p.wptr -> indicates the word pointer in Buffer, therefore if the shift move happens, then -> p.wptr += 1
# p.sptr -> indicates the word pointer in stack,

# movecosts take parser, parser.sentence.head, parser.sentence.deprel as an input

# Here are the steps to do:
# You need to find a correct move set of a given sentence
# You need to make 3 experiments:
# - Use random moves
# - Use correct moves
# - Use predicted moves


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using Pkg
using BinaryBuilder

name = "LibtaskDylib"
version_str = Pkg.TOML.parsefile(joinpath(@__DIR__, "../Project.toml"))["version"] |> strip |> (x) -> lstrip(x, ['v'])
version = VersionNumber(version_str)
# see https://github.com/JuliaPackaging/BinaryBuilder.jl/issues/336
ENV["CI_COMMIT_TAG"] = ENV["TRAVIS_TAG"] = "v" * version_str

event_file = get(ENV, "GITHUB_EVENT_PATH", "")
# run(`cat $event_file`)

# Collection of sources required to build Libtask
function get_commit_id()
    is_pr = get(ENV, "TRAVIS_PULL_REQUEST", "false")
    if is_pr != "false"
        return get(ENV, "TRAVIS_PULL_REQUEST_SHA", "")
    end

    ref = "HEAD"

    gaction = get(ENV, "GITHUB_ACTIONS", "")
    if !isempty(gaction)
        # .pull_request.head.sha, .release.tag_name,
        ref = readlines(`jq --raw-output '.pull_request.head.sha' $event_file`)[1]
        if ref == "null"
            ref = readlines(`jq --raw-output '.release.tag_name' $event_file`)[1]
        end
    end

    if ref == "null"
        ref = "HEAD"
    end
    return readlines(`git rev-parse $ref`)[1]
end

sources = [
    "https://github.com/TuringLang/Libtask.jl.git" => get_commit_id(),
]

# Bash recipe for building across all platforms
script = read(joinpath(dirname(@__FILE__), "build_dylib.sh"), String)

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, :glibc),
    Linux(:x86_64, :glibc),
    Linux(:armv7l, :glibc, :eabihf),
    Linux(:aarch64),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libtask_v1_0", :libtask_v1_0)
    LibraryProduct(prefix, "libtask_v1_1", :libtask_v1_1)
    LibraryProduct(prefix, "libtask_v1_2", :libtask_v1_2)
    LibraryProduct(prefix, "libtask_v1_3", :libtask_v1_3)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
# build_file = "products/build_$(name).v$(version_str).jl"
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
#=##############################################################################################
Filename: wing_cl_alpha_sweep.jl
Author: Ryan Anderson
Contact: rymanderson@gmail.com
README: this is a template file. Convenience methods are provided to prepare a single rotor and
        sweep over advance ratios.
=###############################################################################################

# initialize parameters
"""
    CLAlphaSweep{TF,TS,FUN} <: Parameters

# Fields:

* `CFs::Array{Float64,2}` : [:,j]th element is the [CD,CY,CL] force coefficients of the aircraft at the jth step
* `CMs::Array{Float64,2}` : [:,j]th element is the [CMx,CMy,CMz] moment coefficients of the aircraft at the jth step
* `wake_function::Vector{FUN}`
* `plot_directory::TS`
* `plot_base_name::TS`
* `plot_extension::TS`

"""
struct CLAlphaSweep{TF,TS,FUN} <: Parameters
    CFs::Array{TF,2}
    CMs::Array{TF,2}
    wake_function::Vector{FUN}
    plot_directory::TS
    plot_base_name::TS
    plot_extension::TS
end

# function (rs::RotorSweepParameters{V1, V2, V3, V4, V5})(alphas, omega...)

# end

# rs  = RotorSweepParameter(args...)

# rs(alphas, omega, nblades...)

"""
    cl_alpha_sweep_template(alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip;
        plot_directory = joinpath(topdirectory, "data","plots",TODAY),
        plot_base_name = "default",
        plot_extension = ".pdf",
        step_symbol = L"/alpha[^/circ]", # these forward slashes actually represent backslashes, but those cause an error in the docstrings
        kwargs...)

# Arguments:

* `alphas`
* `wing_b`
* `wing_TR`
* `wing_AR`
* `wing_θroot`
* `wing_θtip`

# Keyword Arguments:

* `plot_directory`
* `plot_base_name`
* `plot_extension`
* `step_symbol`

"""
function cl_alpha_sweep_template(alphas, wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_le_sweep, wing_ϕ;
            plot_directory=joinpath(topdirectory,"data","plots",TODAY),
            plot_base_name="default",
            plot_extension=".pdf",
            step_symbol=L"\alpha[^\circ]",
            kwargs...)

    # prepare subsystems
    wings = simplewing_system(wing_b, wing_TR, wing_AR, wing_θroot, wing_θtip, wing_le_sweep, wing_ϕ; kwargs...)
    rotors = CCBladeSystem(
        Vector{CC.Rotor}(undef,0),
        Vector{Vector{CC.Section{Float64, Float64, Float64, nothing}}}(undef,0),
        Vector{Vector{Float64}}(undef,0),
        Vector{Int64}(undef,0),
        Vector{Vector{Float64}}(undef,0),
        Vector{Vector{Float64}}(undef,0),
        Vector{Bool}(undef,0)
    )

    # rotors::Vector{CC.Rotor}#{TF, Int64, Bool, T1, T2, T3, T4}}
    # sectionlists::Vector{Vector{CCBlade.Section{TF1, TF2, TF3, TAF}}}
    # rlists::Vector{Vector{TF}}
    # index::Vector{Int64}
    # positions::Vector{Vector{TF}}
    # orientations::Vector{Vector{TF}}
    # spin_directions::Vector{Bool}
    inertia = nothing
    nonliftingbodies = nothing
    structures = nothing
    motors = nothing
    batteries = nothing

    # build system struct
    aircraft = Aircraft(wings, rotors, inertia, nonliftingbodies, structures, motors, batteries)

    # compile actions
    actions = [solve_wing_CF_CM]

    # initialize parameters
    wake_function, CFs, CMs = solve_wing_CF_CM(aircraft, alphas) # wake_function, CFs, CMs

    # check sizes and instantiate struct
    @assert size(CFs) == size(CMs) "size of CFs and CMs inconsistent"
    @assert size(CFs)[2] == length(alphas) "length of parameter CFs and alphas inconsistent"
    @assert size(CFs)[1] == 3 "first dimension of CFs should be of length 1"

    parameters = CLAlphaSweep(CFs, CMs, wake_function, plot_directory, plot_base_name, plot_extension)

    # build freestream_function
    function freestream_function(aircraft, parameters, environment, alphas, stepi)

        # calculate freestream
        Vinf = 1.0 # arbitrary for CL-alpha sweep
        alpha = alphas[stepi]
        beta = 0.0
        Omega = StaticArrays.@SVector zeros(3)
        freestream = Freestream(Vinf, alpha, beta, Omega)

        return freestream
    end

    # build environment_function
    function environment_function(aircraft, parameters, alphas, stepi)
        Environment() # arbitrary
    end

    # compile postactions
    postactions = [post_plot_cf_cm_alpha_sweep]

    # build objective_function
    objective_function(aircraft, parameters, freestream, environment, alphas) = 0.0

    return aircraft, parameters, actions, freestream_function, environment_function, postactions, objective_function, alphas, step_symbol
end

max_stars_repo_name: byuflowlab/AircraftSystems

max_issues_repo_name: byuflowlab/AircraftSystems

max_forks_repo_name: byuflowlab/AircraftSystems


text:
export
    solve, evaluate

# http://www.alpertron.com.ar/METHODS.HTM

include("diophantine_solution.jl")
include("diophantine_linear.jl")
include("diophantine_quadratic.jl")


text:
# Rose Teague 06/07/2018
# AdiabaticPropagationTest.jl
#
# Runs the Adiabatic propagation of an eigenstate
push!(LOAD_PATH,"../src/")

using TheDancer
using Base.Test

N, Edisorder, Jdisorder, modelJ, B, dipolestrength, r = init!(10, 0.0, 0.0, 0.001, 298, 20)

"""
Test full time-dependent propagation for a simple 3-site system. With small transfer
integral (J=0.001eV).

Unitary propagation matrix should be unitary - test U'U = I

Prepare system in an eigenstate, unitary propagation should leave psi.^2 unchanged.
"""
function main()
    S,E,H,psi,dipoles=prepare_model()
    density = abs.(psi.^2)
    psi_1 = psi

    for i in 1:200
        psi,U=TimeDependentPropagation(psi_1,H,1,decompose=false,verbose=false,test=true)
        density = abs.(psi.^2)
        unit = sum(U'U)/N
        e = 1E-16
        @test 1 ≈ unit atol = e
    end

    diff = sum(density - abs.(psi_1.^2))

    ε = 1E-7
    @test  0 ≈ diff atol = ε
end

main()


text:
"""A package for graph computations related to graph Laplacians

Graphs are represented by sparse adjacency matrices, etc.
"""
module Laplacians

  using LinearAlgebra, SparseArrays, Statistics

  include("approxCholTypes.jl")
  include("approxChol.jl")
  include("graphOps.jl")
  include("graphUtils.jl")

  export adj, lap, extendMatrix, approxchol_lap_pc, LDLsolver!, LDLinv

end # module Laplacians.jl


text:
# ---
# title: 854. K-Similar Strings
# id: problem854
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Breadth-first Search, Graph
# link: <https://leetcode.com/problems/k-similar-strings/description/>
# hidden: true
# ---
#
# Strings `A` and `B` are `K`-similar (for some non-negative integer `K`) if we
# can swap the positions of two letters in `A` exactly `K` times so that the
# resulting string equals `B`.
#
# Given two anagrams `A` and `B`, return the smallest `K` for which `A` and `B`
# are `K`-similar.
#
# **Example 1:**
#
#
#
#     Input: A = "ab", B = "ba"
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: A = "abc", B = "bca"
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: A = "abac", B = "baca"
#     Output: 2
#
#
# **Example 4:**
#
#
#
#     Input: A = "aabc", B = "abca"
#     Output: 2
#
# **Note:**
#
#   1. `1 <= A.length == B.length <= 20`
#   2. `A` and `B` contain only lowercase letters from the set `{'a', 'b', 'c', 'd', 'e', 'f'}`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Convert vector values to quantile scores
function quantiles(vec::Array; ignorenothing::Bool=true)
    if ignorenothing == true
        qfinder = ecdf(filter(!isnothing, vec))
        qvec = replace(x -> !isnothing(x) ? qfinder(x) : x, vec)
    else
        qfinder = ecdf(vec)
        qvec = qfinder(vec)
    end
    return qvec
end

# Convert & mutate vector values to quantile scores
function quantiles!(vec::Array; ignorenothing::Bool=true)
    if ignorenothing == true
        qfinder = ecdf(filter(!isnothing, vec))
        replace!(x -> !isnothing(x) ? qfinder(x) : x, vec)
    else
        qfinder = ecdf(vec)
        replace!(x -> qfinder(x), vec)
    end
    return vec
end

# Convert SimpleSDMLayer values to quantile scores
function quantiles(layer::SimpleSDMLayer)
    qfinder = ecdf(Float32.(filter(!isnothing, layer.grid)))
    qgrid = replace(x -> isnothing(x) ? x : qfinder(x), layer.grid)
    qlayer = SimpleSDMResponse(qgrid, layer)
    return qlayer
end

# Convert & mutate SimpleSDMLayer values to quantile scores
function quantiles!(layer::SimpleSDMLayer)
    qfinder = ecdf(Float32.(filter(!isnothing, layer.grid)))
    replace!(x -> isnothing(x) ? x : qfinder(x), layer.grid)
    return layer
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
### expm

function codeblock(stmts...)
    stmts_ = Any[]
    for s in stmts
        if !(s == nothing)
            push!(stmts_, s)
        end
    end
    isempty(stmts_) ? nothing : Expr(:block, stmts_...)
end

macro rotate3(A, c, s, i,j)
  l = setdiff(1:3, [i,j])[1]
  codeblock([
        quote
             $(symbol("a",l,k)) = $(esc(A))[$l,$k]
             $(symbol("a",i,k))  = $(esc(A))[$i,$k] * $(esc(c)) + $(esc(A))[$j,$k] * $(esc(s))
             $(symbol("a",j,k))  = $(esc(A))[$j,$k] * $(esc(c)) - $(esc(A))[$i,$k] * $(esc(s))
        end
        for k in 1:3]...,        quote Mat(((a11,a21,a31),(a12,a22,a32),(a13,a23,a33))) end
    )
end

macro rotate3t(A, c, s, i,j)
  l = setdiff(1:3, [i,j])[1]
  codeblock([
        quote
             $(symbol("a",k,l)) = $(esc(A))[$k,$l]
             $(symbol("a",k,i)) = $(esc(A))[$k,$i] * $(esc(c)) + $(esc(A))[$k,$j] * $(esc(s))
             $(symbol("a",k,j)) = $(esc(A))[$k,$j] * $(esc(c)) - $(esc(A))[$k,$i] * $(esc(s))
        end
        for k in 1:3]...,        quote Mat(((a11,a21,a31),(a12,a22,a32),(a13,a23,a33))) end
    )
end

macro dblrotate3(A, c, s, i,j)
  k = setdiff(1:3, [i,j])[1]
  quote
             $(symbol("a",i,k)) = $(esc(A))[$i,$k] * $(esc(c)) + $(esc(A))[$j,$k] * $(esc(s))
             $(symbol("a",j,k)) = $(esc(A))[$j,$k] * $(esc(c)) - $(esc(A))[$i,$k] * $(esc(s))
             $(symbol("a",k,i)) = $(esc(A))[$k,$i] * $(esc(c)) + $(esc(A))[$k,$j] * $(esc(s))
             $(symbol("a",k,j)) = $(esc(A))[$k,$j] * $(esc(c)) - $(esc(A))[$k,$i] * $(esc(s))

             $(symbol("a",j,j)) = $(esc(A))[$j,$j] * abs2($(esc(c))) + $(esc(A))[$i,$i] * abs2($(esc(s))) - ($(esc(A))[$i,$j] + $(esc(A))[$j,$i]) * $(esc(s)) * $(esc(c))
             $(symbol("a",j,i)) = $(esc(A))[$j,$i] * abs2($(esc(c))) - $(esc(A))[$i,$j] * abs2($(esc(s))) + ($(esc(A))[$j,$j] - $(esc(A))[$i,$i]) * $(esc(c)) * $(esc(s))
             $(symbol("a",i,j)) = $(esc(A))[$i,$j] * abs2($(esc(c))) - $(esc(A))[$j,$i] * abs2($(esc(s))) + ($(esc(A))[$j,$j] - $(esc(A))[$i,$i]) * $(esc(c)) * $(esc(s))
             $(symbol("a",i,i)) = $(esc(A))[$i,$i] * abs2($(esc(c))) + $(esc(A))[$j,$j] * abs2($(esc(s))) + ($(esc(A))[$i,$j] + $(esc(A))[$j,$i]) * $(esc(s)) * $(esc(c))

             $(symbol("a",k,k)) = $(esc(A))[$k,$k]

            Mat(((a11,a21,a31),(a12,a22,a32),(a13,a23,a33)))
   end
end
function hessenberg3(A)
    c, s, _ = LinAlg.givensAlgorithm(A[2,1],A[3,1])
    @dblrotate3(A, c, s, 2, 3)
end

function reflect3{T}(x1,x2,x3, τ::Number, A::Mat{3, 3, T}) # apply reflector from left

            vA1 = τ'*(A[1, 1] +  x2'*A[2, 1] +  x3'*A[3, 1])
            a11 = A[1,1] - vA1
            a21 = A[2,1] - x2*vA1
            a31 = A[3,1] - x3*vA1

            vA2 = τ'*(A[1, 2] +  x2'*A[2, 2] +  x3'*A[3, 2])
            a12 = A[1,2] - vA2
            a22 = A[2,2] - x2*vA2
            a32 = A[3,2] - x3*vA2

            vA3 = τ'*(A[1, 3] +  x2'*A[2, 3] +  x3'*A[3, 3])
            a13 = A[1,3] - vA3
            a23 = A[2,3] - x2*vA3
            a33 = A[3,3] - x3*vA3
            Mat(((a11,a21,a31),(a12,a22,a32),(a13,a23,a33)))::Mat{3, 3, T}

end

function reflect3t{T}(x1,x2,x3, τ::Number, A::Mat{3, 3, T}) # apply reflector from left

            vA1 = τ'*(A[1, 1] +  x2'*A[1, 2] +  x3'*A[1, 3])
            a11 = A[1,1] - vA1
            a21 = A[1,2] - x2*vA1
            a31 = A[1,3] - x3*vA1

            vA2 = τ'*(A[2, 1] +  x2'*A[2, 2] +  x3'*A[2, 3])
            a12 = A[2,1] - vA2
            a22 = A[2,2] - x2*vA2
            a32 = A[2,3] - x3*vA2

            vA3 = τ'*(A[3, 1] +  x2'*A[3, 2] +  x3'*A[3, 3])
            a13 = A[3,1] - vA3
            a23 = A[3,2] - x2*vA3
            a33 = A[3,3] - x3*vA3
            Mat(((a11,a12,a13),(a21,a22,a23),(a31,a32,a33)))::Mat{3, 3, T}

end

#francis QR algorithm with single or double shift

function francisdbl{T<:Real}(AA::Mat{3, 3, T})
    A = hessenberg3(AA)
    EPS = eps(T)
    i = 0
    loc = 0 #location of zero
    la1::T = la2::T = la3::T = 0.
    for i in 1:200
        if abs(A[3,2]) <= EPS * (abs(A[2,2])+abs(A[3,3]))
            la1 = A[3,3]
            loc = 1
            break
        elseif abs(A[2,1]) <= EPS * (abs(A[1,1])+abs(A[2,2]))
            la1 = A[1,1]
            loc = 3
            break
        end
        l = 3
        # if real ev in sub 2x2 matrix do twice a raleigh shift
        if 0 < abs2(A[l,l] - A[2,2]) + 4*A[l,2]*A[2,l]
            mu = A[l,l]
            A = A - mu*(eye(typeof(A))::Mat{3, 3, T})
            c, s, _ = LinAlg.givensAlgorithm(A[1,1],A[2,1])::Tuple{T,T,T}
            A2 = @rotate3(A, c, s, 1, 2)
            c2, s2, _ = LinAlg.givensAlgorithm(A2[2,2],A2[3,2])::Tuple{T,T,T}
            A = @rotate3(A2, c2, s2, 2, 3)

            A2 = @rotate3t(A, c, s, 1, 2)
            A = @rotate3t(A2, c2, s2, 2, 3)

            c, s, _ = LinAlg.givensAlgorithm(A[1,1],A[2,1])::Tuple{T,T,T}
            A2 = @rotate3(A, c, s, 1, 2)
            c2, s2, _ = LinAlg.givensAlgorithm(A2[2,2],A2[3,2])::Tuple{T,T,T}
            A = @rotate3(A2, c2, s2, 2, 3)

            A2 = @rotate3t(A, c, s, 1, 2)
            A = @rotate3t(A2, c2, s2, 2, 3)
            A = A + mu*(eye(typeof(A))::Mat{3, 3, T})
        else #do an implicit francis double shift

            x = (A[1,1]*A[1,1] + A[1,2]*A[2,1] + A[3,1]*A[1,3]) -  A[1,2]*A[2,1] -
                (A[3,3] + A[2,2])*A[1,1]  + A[3,3] * A[2,2] - A[2,3] * A[3,2]
            y = A[2,1]*(A[1,1] - A[3,3])
            z = A[2,1]*A[3,2]

            ξ1 = x
            normu = abs2(x) + abs2(y) + abs2(z)
            if normu == zero(normu)
               g = zero(ξ1/normu)
            else
                normu = sqrt(normu)
                ν = copysign(normu, real(ξ1))
                ξ1 += ν
                x = -ν
                y /= ξ1
                z /= ξ1
                g = ξ1/ν
            end
            A = reflect3(x,y,z, g, A)
            A = reflect3t(x,y,z, g, A)
            A = hessenberg3(A)
        end
    end

    # check convergence
    if loc == 0
        converged = false
        loc = 1
    else
        converged = true
    end

    # compute eigenvalues of 2x2 submatrix
    spur = A[loc,loc] + A[2,2]
    det2 = A[loc,loc] * A[2,2] - A[2,loc] * A[loc,2]
    dis = abs2(A[loc,loc] - A[2,2]) + 4*A[loc,2]*A[2,loc]

    # return eigenvalue, discriminant of 2x2 matrix of end result and convergence status
    la1, 0.5*spur, dis/4, converged

end

# default version using conversion to and from Matrix type
expm{N,T}(m::Mat{N,N,T}) = Mat{N,N,T}(expm(convert(Matrix{T}, m)))

expm{T}(A::Mat{1, 1, T}) = Mat{1, 1, T}(((expm(A[1,1]),),))
function expm{T<:Complex}(A::Mat{2, 2, T})
 	a = A[1,1]
	b = A[1,2]
	c = A[2,1]
	d = A[2,2]

	z = sqrt((a-d)*(a-d) + 4.0*b*c )
	e = exp(a/2.0 + d/2.0 - z/2.0)
	f = exp(a/2.0 + d/2.0 + z/2.0)

    Mat{2, 2, T}(
        ( -(e*(a - d - z))/(2.0* z) + (f*(a - d + z))/(2.0* z), -((e * c)/z) + (f * c)/z),
        ( -((e * b)/z) + (f * b)/z,	-(e*(-a + d - z))/(2.0* z) + (f*(-a + d + z))/(2.0* z))
    )
end

# presumably better without resorting to complex numbers, but for now...
function expm{T<:Real}(A::Mat{2, 2, T})
 	a = A[1,1]
	b = A[1,2]
	c = A[2,1]
	d = A[2,2]

	z = sqrt(Complex((a-d)*(a-d) + 4.0*b*c))
	e = exp(a/2.0 + d/2.0 - z/2.0)
	f = exp(a/2.0 + d/2.0 + z/2.0)

    Mat{2, 2, T}(
        ( real(-(e*(a - d - z))/(2.0* z) + (f*(a - d + z))/(2.0* z)), real(-((e * c)/z) + (f * c)/z)),
        ( real(-((e * b)/z) + (f * b)/z), real(-(e*(-a + d - z))/(2.0* z) + (f*(-a + d + z))/(2.0* z)))
    )
end

function expm{T<:Real}(AA::Mat{3, 3, T})
    t = sqrt(sum(AA.^2))

    A = AA/t #normalize

    la1, x, dis, conv = francisdbl(A)
    if !conv # if convergence takes to long, call lapack
       return Mat{3,3,T}(expm(convert(Matrix{T}, AA)))
    end
    putzer(t, A, la1, x, dis)
end

function dexp(t, x, y)
    if abs2(x-y) > eps()
        2exp(t*(x+y)/2)*sinh(t*(x-y)/2)/((x-y))
    else
        exp(t*(x+y)/2)*(t+(t^3*(x-y)^2)/24)
    end
end

function ddexp(t, x, y, z)
    d, i = findmax([abs2(y-x), abs2(z-y),abs2(z-x)])
    if d < eps()
        exp(t*(x+y+z)/3)*t^2/2
    elseif i == 1
        (dexp(t, y, z) - dexp(t, x,z))/(y-x)
    elseif i==2
        (dexp(t, z, x) - dexp(t, y,x))/(z-y)
    elseif i==3
        (dexp(t, z, y) - dexp(t, x,y))/(z-x)
    end
end

function putzer{T}(t, A::Mat{3,3,T}, x, y, d)
    deltaabs = sqrt(abs(d))
    delta = sqrt(complex(d))
    r1 = real(exp(t * (y+delta)))
    la1 = (y+delta)
    la2 = (y-delta)
    la3 = x
    r2 = real(dexp(t, la1, la2))
    r3 = real(ddexp(t, la1, la2, la3))
    R = r2 - la2*r3
    r3 * A*A +  real(R - la1*r3)*A +  real(r1 - la1 * R)*eye(Mat{3,3,T})
end


text:
module TestExecutorTypes

include("preamble.jl")
using Referenceables: referenceable
using Transducers: Executor, PreferParallel, executor_type, maybe_set_simd

struct Opinionated{T<:Executor} <: AbstractVector{Int} end
Opinionated(::Type{T}) where {T} = Opinionated{T}()
Transducers.executor_type(::Opinionated{T}) where {T} = T

# Minimal interface to make printing happy:
Base.size(::Opinionated) = (0,)

@testset "executor_type (promote)" begin
    @testset for (xs, ex) in [
        (1:2, PreferParallel),
        (Channel(0), SequentialEx),
        (pairs(1:2), PreferParallel),
        (pairs(Opinionated(DistributedEx)), DistributedEx),
        (zip(1:2, 3:4), PreferParallel),
        (zip(1:2, Channel(0)), SequentialEx),
        (Iterators.product(1:2, 3:4), PreferParallel),
        (Iterators.product(1:2, Channel(0)), SequentialEx),
        (Opinionated(DistributedEx), DistributedEx),
        (zip(Opinionated(DistributedEx), 1:2), DistributedEx),
        (zip(1:2, Opinionated(DistributedEx)), DistributedEx),
        (zip(Opinionated(DistributedEx), 1:2, 3:4), DistributedEx),
        (zip(1:2, Opinionated(DistributedEx), 3:4), DistributedEx),
        (zip(1:2, 3:4, Opinionated(DistributedEx)), DistributedEx),
        (referenceable([]), PreferParallel),
        (referenceable(Opinionated(SequentialEx)), SequentialEx),
    ]
        @test @inferred(executor_type(xs)) === ex
    end
end

@testset "executor_type (error)" begin
    err =
        @test_error executor_type(zip(Opinionated(ThreadedEx), Opinionated(DistributedEx)))
    msg = sprint(showerror, err)
    @test occursin("failed to promote executors", msg)
end

@testset "show" begin
    @test occursin("ThreadedEx()", sprint(show, ThreadedEx()))
    @test occursin("ThreadedEx(simd = true)", sprint(show, ThreadedEx(simd = true)))
    @test occursin(
        "ThreadedEx(simd = true, basesize = 2)",
        sprint(show, ThreadedEx(simd = true, basesize = 2)),
    )
end

@testset "repr" begin
    @testset for x in [
        ThreadedEx(),
        ThreadedEx(simd = true),
        ThreadedEx(simd = true, basesize = 2),
        ThreadedEx{Any}(),
        ThreadedEx{Any}(simd = true),
        ThreadedEx{Any}(simd = true, basesize = 2),
    ]
        str = repr(x)
        reconstructed = include_string(@__MODULE__, str)
        @test reconstructed == x
    end
end

@testset "maybe_set_simd" begin
    @test maybe_set_simd(SequentialEx(), true).kwargs.simd === true
    @test maybe_set_simd(SequentialEx(), false).kwargs.simd === false
    @test maybe_set_simd(SequentialEx(simd = false), true).kwargs.simd === false
    @test maybe_set_simd(SequentialEx(simd = true), false).kwargs.simd === true
end

end  # module


text:
using MonteCarlo, MonteCarloObservable, StableDQMC
using Test
using LinearAlgebra, Random, Dates
using MonteCarlo: @bm, TimerOutputs

# In case some test failed and left behind a .jld file
for f in readdir()
    if endswith(f, ".jld") || endswith(f, "jld2")
        @warn "Removing $f"
        rm(f)
    end
end

macro benchmark_test(name, code)
    TimerOutputs.timer_expr(MonteCarlo, true, name, code)
end

@testset "All Tests" begin
    @testset "Utilities" begin
        @bm function test1(x, y)
            sleep(x+y)
        end
        @bm test2(x, y) = sleep(x+y)
        # @eval MonteCarlo ... will put the body of test3, test4 in the wrong scope
        # eval(timer_expr(...)) doesn't work
        function test3(x, y)
            @benchmark_test "test1(::Any, ::Any)" begin sleep(x+y) end
        end
        test4(x, y) = @benchmark_test "test2(::Any, ::Any)" begin sleep(x+y) end

        x = code_lowered(test1, Tuple{Float64, Float64})[1]
        y = code_lowered(test3, Tuple{Float64, Float64})[1]
        @test x.code == y.code

        x = code_lowered(test2, Tuple{Float64, Float64})[1]
        y = code_lowered(test4, Tuple{Float64, Float64})[1]
        @test x.code == y.code

        @test !MonteCarlo.timeit_debug_enabled()
        enable_benchmarks()
        @test MonteCarlo.timeit_debug_enabled()
        disable_benchmarks()
        @test !MonteCarlo.timeit_debug_enabled()
    end

    @testset "Lattices" begin
        include("lattices.jl")
    end

    @testset "Models" begin
        include("modeltests_IsingModel.jl")
        include("modeltests_HubbardModel.jl")
    end

    @testset "Flavors" begin
        # include("flavortests_MC.jl")
        include("flavortests_DQMC.jl")
    end

    @testset "Measurements" begin
        include("measurements.jl")
    end

    @testset "Intergration tests" begin
        include("integration_tests.jl")
    end

    @testset "Exact Diagonalization" begin
        include("ED/ED_tests.jl")
    end

    @testset "File IO" begin
        include("FileIO.jl")
    end
end


text:
using GeometricFlux
using GraphSignals
using Flux
using Flux: throttle
using Flux.Losses: logitbinarycrossentropy
using Flux: @epochs
using JLD2
using Statistics
using SparseArrays
using LightGraphs.SimpleGraphs
using LightGraphs: adjacency_matrix

@load "data/cora_features.jld2" features
@load "data/cora_graph.jld2" g

num_nodes = 2708
num_features = 1433
h_dim = 32
z_dim = 16
target_catg = 7
epochs = 200

## Preprocessing data
masks = [rand(Float32, num_nodes, num_nodes).>0.1 for i in 1:10]
adj_mat = Matrix{Float32}(adjacency_matrix(g))
train_data = [(FeaturedGraph(adj_mat.*M, Matrix{Float32}(features)), adj_mat) for M in masks]

## Model
model = VGAE(GCNConv(num_features=>h_dim, relu; cache=false), h_dim, z_dim, σ)
encoder = model.encoder
decoder = model.decoder
ps = Flux.params(model)

l2_norm(p) = sum(abs2, p)

function loss(fg, Y, X=node_feature(fg), T=eltype(X), β=one(T), λ=T(0.01); debug=false)
    μ̂, logσ̂ = summarize(encoder, fg)
    Z = node_feature(encoder(fg))
    kl_q_p = -T(0.5) * sum(one(T) .+ T(2).*logσ̂ .- μ̂.^2 .- exp.(T(2).*logσ̂))
    logp_y_z = -sum(logitbinarycrossentropy(decoder(Z), Y, agg=identity)) / size(Y,2)
    l2reg = sum(l2_norm, ps)
    debug && begin
        @show kl_q_p
        @show logp_y_z
        @show l2reg
    end
    -logp_y_z + β*kl_q_p + λ*l2reg
end

average_loss(data) = mean(map(x -> loss(x...), data))

## Training
opt = ADAM(0.01)
evalcb() = @show(average_loss(train_data))

@epochs epochs Flux.train!(loss, ps, train_data, opt, cb=throttle(evalcb, 10))


text:
using Test, HMCResearchRandomWalks.PPD

"testing functions for locations"
function rand_x() return PPD.ESCAPE_X * rand() * 2 - PPD.ESCAPE_X end
function rand_y() return PPD.ESCAPE_Y * rand() end # greater than 0.0

@testset "new spawnrandompoint" begin
    for _ in 1:10
        initxy = PPD.spawnrandompoint()
        @test -12.0 <= initxy[1] <= 12.0
        @test 0.2 <= initxy[2] <= 1.0
    end
end

@testset "locations" begin
    @test PPD.insensor([rand_x(), 0.0])
    @test !PPD.insensor([rand_x(), 0.1])

    @test PPD.inwalls(13.0, 1.499)
    @test !PPD.inwalls(13.0, 1.501)
    @test !PPD.inwalls(12.0, 1.5)

    @test PPD.inppd([0.0, 0.15])
    @test !PPD.inppd([0.0, 0.16])
    @test !PPD.inppd([12.6, 0.15])

    @test PPD.inwater([0.0, 0.16])
    @test !PPD.inwater([13.0, 1.499])
    @test !PPD.inwater([0.0, 0.15])

    randxy = [rand_x(), rand_y()]
    @test !PPD.inenz(randxy)
    @test PPD.whereoutsidespawn(randxy) == "N"
    @test PPD.sensorcases(rand_x()) == "sensor"
end

@testset "boundarycheck cases" begin
    initxy = [15.0, 1.56]
    dx = 0.0
    dy = -1.0
    ss = 0.1
    output = PPD.wallcases(initxy, dx, dy, ss)
    @test output[1] == [15.0, 1.51]
    @test output[2] == "no collision"
end

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
using PicoQuant

# Create a tensor network circuit to find contraction plans for.
qasm_str = """OPENQASM 2.0;
              include "qelib1.inc";
              qreg q[3];
              h q[0];
              cx q[0],q[1];
              cx q[1],q[2];
              h q[2];
              cx q[0],q[1];
              cx q[0],q[2];"""

circ = load_qasm_as_circuit(qasm_str)
tng = convert_qiskit_circ_to_network(circ)
add_input!(tng, "000")
add_output!(tng, "000")

# %% Plot the optimal contraction tree
optimal_plan = netcon(tng)
optimal_contraction_tree = create_contraction_tree(tng, optimal_plan)

display(plot_contraction_tree(optimal_contraction_tree))

# %% Plot the random contration tree found by bgreedy
bgreedy_plan, time_cost = bgreedy(tng, 0.5, 2.5, 10)
bgreedy_contraction_tree = create_contraction_tree(tng, bgreedy_plan)

display(plot_contraction_tree(bgreedy_contraction_tree))


text:
module TestErrors

using Test
using MathOptInterface
const MOI = MathOptInterface

include("dummy.jl")

function test_errors_fallback_AddVariableNotAllowed()
    model = DummyModel()
    @test_throws MOI.AddVariableNotAllowed MOI.add_variable(model)
    try
        MOI.add_variable(model)
    catch err
        @test sprint(showerror, err) ==
              "MathOptInterface.AddVariableNotAllowed:" *
              " Adding variables cannot be performed. You may want to use a" *
              " `CachingOptimizer` in `AUTOMATIC` mode or you may need to call" *
              " `reset_optimizer` before doing this operation if the" *
              " `CachingOptimizer` is in `MANUAL` mode."
    end
    @test_throws MOI.AddVariableNotAllowed MOI.add_variables(model, 2)
    return
end

function test_errors_inconsistent_vectorscalar()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    @test_throws(
        MOI.ErrorException,
        MOI.add_constraint(
            model,
            MOI.VectorOfVariables([vi, vi]),
            MOI.EqualTo(0),
        )
    )
    @test_throws(
        MOI.ErrorException,
        MOI.add_constraint(model, vi, MOI.Nonnegatives(2))
    )
    return
end

function _test_errors_UnsupportedConstraint(f)
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    @test_throws(MOI.UnsupportedConstraint, f(model, vi, MOI.EqualTo(0)),)
    try
        f(model, vi, MOI.EqualTo(0))
    catch err
        @test sprint(showerror, err) ==
              "$(MOI.UnsupportedConstraint{MOI.VariableIndex,MOI.EqualTo{Int}}):" *
              " `$MOI.VariableIndex`-in-`$MOI.EqualTo{$Int}` constraint is" *
              " not supported by the model."
    end
    return
end

function test_errors_UnsupportedConstraint()
    _test_errors_UnsupportedConstraint(MOI.add_constraint)
    return
end

function test_errors_UnsupportedConstraint_shortcut()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    @test_throws MOI.UnsupportedConstraint begin
        MOI.add_constraint(model, vi, MOI.EqualTo(0))
    end
    @test_throws MOI.UnsupportedConstraint begin
        MOI.add_constraint(model, [vi, vi], MOI.Nonnegatives(2))
    end
    @test_throws MOI.UnsupportedConstraint begin
        MOI.add_constraints(model, [vi, vi], [MOI.EqualTo(0), MOI.EqualTo(0)])
    end
end

function test_errors_add_constraint()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    @test_throws(
        MOI.AddConstraintNotAllowed,
        MOI.add_constraint(model, vi, MOI.EqualTo(0.0)),
    )
    try
        MOI.add_constraint(model, vi, MOI.EqualTo(0.0))
    catch err
        @test sprint(showerror, err) ==
              "$(MOI.AddConstraintNotAllowed{MOI.VariableIndex,MOI.EqualTo{Float64}}):" *
              " Adding `$MOI.VariableIndex`-in-`$MOI.EqualTo{Float64}`" *
              " constraints cannot be performed. You may want to use a" *
              " `CachingOptimizer` in `AUTOMATIC` mode or you may need to call" *
              " `reset_optimizer` before doing this operation if the" *
              " `CachingOptimizer` is in `MANUAL` mode."
    end
    @test_throws(
        MOI.AddConstraintNotAllowed,
        MOI.add_constraint(model, vi, MOI.EqualTo(0.0)),
    )
    @test_throws(
        MOI.AddConstraintNotAllowed,
        MOI.add_constraint(model, [vi, vi], MOI.Zeros(2)),
    )
    @test_throws(
        MOI.AddConstraintNotAllowed,
        MOI.add_constraints(
            model,
            [vi, vi],
            [MOI.EqualTo(0.0), MOI.EqualTo(0.0)],
        ),
    )
    return
end

function test_errors_DeleteNotAllowed()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    @test_throws MOI.DeleteNotAllowed{typeof(vi)} MOI.delete(model, vi)
    try
        MOI.delete(model, vi)
    catch err
        @test sprint(showerror, err) ==
              "MathOptInterface.DeleteNotAllowed{MathOptInterface.VariableIndex}:" *
              " Deleting the index MathOptInterface.VariableIndex(1) cannot be" *
              " performed. You may want to use a `CachingOptimizer` in" *
              " `AUTOMATIC` mode or you may need to call `reset_optimizer`" *
              " before doing this operation if the `CachingOptimizer` is in" *
              " `MANUAL` mode."
    end
    @test_throws MOI.DeleteNotAllowed{typeof(ci)} MOI.delete(model, ci)
    try
        MOI.delete(model, ci)
    catch err
        @test sprint(showerror, err) ==
              "$(MathOptInterface.DeleteNotAllowed{MathOptInterface.ConstraintIndex{MathOptInterface.VariableIndex,MathOptInterface.EqualTo{Float64}}}):" *
              " Deleting the index $(MathOptInterface.ConstraintIndex{MathOptInterface.VariableIndex,MathOptInterface.EqualTo{Float64}}(1))" *
              " cannot be performed. You may want to use a `CachingOptimizer`" *
              " in `AUTOMATIC` mode or you may need to call `reset_optimizer`" *
              " before doing this operation if the `CachingOptimizer` is in" *
              " `MANUAL` mode."
    end
    return
end

function _test_unsupported_attribute(f)
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    @test_throws(
        MOI.UnsupportedAttribute,
        f(model, MOI.ObjectiveFunction{MOI.VariableIndex}(), vi),
    )
    @test_throws(
        MOI.UnsupportedAttribute,
        f(model, MOI.ConstraintDualStart(), ci, 0.0),
    )
    return
end

function test_errors_unsupported_attribute()
    _test_unsupported_attribute(MOI.set)
    return
end

function test_errors_SetAttributeNotAllowed()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    model = DummyModel()
    @test_throws(
        MOI.SetAttributeNotAllowed,
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE),
    )
    @test_throws(
        MOI.SetAttributeNotAllowed,
        MOI.set(model, MOI.ConstraintPrimalStart(), ci, 0.0),
    )
    return
end

function test_errors_ConstraintFunction_NotAllowed()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    @test_throws(
        MOI.SetAttributeNotAllowed,
        MOI.set(model, MOI.ConstraintFunction(), ci, vi)
    )
    @test_throws(
        ArgumentError,
        MOI.set(
            model,
            MOI.ConstraintFunction(),
            ci,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, vi)], 0.0),
        ),
    )
    return
end

function test_errors_ConstraintSet_NotAllowed()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    @test_throws(
        MOI.SetAttributeNotAllowed,
        MOI.set(model, MOI.ConstraintSet(), ci, MOI.EqualTo(1.0))
    )
    @test_throws(
        ArgumentError,
        MOI.set(model, MOI.ConstraintSet(), ci, MOI.EqualTo(1))
    )
    @test_throws(
        ArgumentError,
        MOI.set(model, MOI.ConstraintSet(), ci, MOI.GreaterThan(1.0))
    )
end

function test_errors_ModifyNotAllowed_constraint()
    model = DummyModel()
    vi = MOI.VariableIndex(1)
    ci = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Float64}}(1)
    change = MOI.ScalarConstantChange(1.0)
    err = MOI.ModifyConstraintNotAllowed(ci, change)
    @test_throws err MOI.modify(model, ci, change)
    @test sprint(showerror, err) ==
          "$(MathOptInterface.ModifyConstraintNotAllowed{MathOptInterface.VariableIndex,MathOptInterface.EqualTo{Float64},MathOptInterface.ScalarConstantChange{Float64}}):" *
          " Modifying the constraints $(MathOptInterface.ConstraintIndex{MathOptInterface.VariableIndex,MathOptInterface.EqualTo{Float64}}(1))" *
          " with MathOptInterface.ScalarConstantChange{Float64}(1.0) cannot" *
          " be performed. You may want to use a `CachingOptimizer` in" *
          " `AUTOMATIC` mode or you may need to call `reset_optimizer`" *
          " before doing this operation if the `CachingOptimizer` is in" *
          " `MANUAL` mode."
end

function test_errors_ModifyNotAllowed_objective()
    model = DummyModel()
    change = MOI.ScalarConstantChange(1.0)
    attr = MOI.ObjectiveFunction{MOI.VariableIndex}()
    err = MOI.ModifyObjectiveNotAllowed(change)
    @test_throws err MOI.modify(model, attr, change)
    @test sprint(showerror, err) ==
          "$(MathOptInterface.ModifyObjectiveNotAllowed{MathOptInterface.ScalarConstantChange{Float64}}):" *
          " Modifying the objective function with $(MathOptInterface.ScalarConstantChange{Float64}(1.0))" *
          " cannot be performed. You may want to use a `CachingOptimizer`" *
          " in `AUTOMATIC` mode or you may need to call `reset_optimizer`" *
          " before doing this operation if the `CachingOptimizer` is in" *
          " `MANUAL` mode."
end

function test_errors_show_SetAttributeNotAllowed()
    @test sprint(showerror, MOI.UnsupportedAttribute(MOI.Name())) ==
          "$MOI.UnsupportedAttribute{$MOI.Name}:" *
          " Attribute $MOI.Name() is not supported by the model."
    @test sprint(showerror, MOI.UnsupportedAttribute(MOI.Name(), "Message")) ==
          "$MOI.UnsupportedAttribute{$MOI.Name}:" *
          " Attribute $MOI.Name() is not supported by the model: Message"
    @test sprint(showerror, MOI.SetAttributeNotAllowed(MOI.Name())) ==
          "$MOI.SetAttributeNotAllowed{$MOI.Name}:" *
          " Setting attribute $MOI.Name() cannot be performed. You may want to use" *
          " a `CachingOptimizer` in `AUTOMATIC` mode or you may need to call" *
          " `reset_optimizer` before doing this operation if the" *
          " `CachingOptimizer` is in `MANUAL` mode."
    @test sprint(
              showerror,
              MOI.SetAttributeNotAllowed(MOI.Name(), "Message"),
          ) ==
          "$MOI.SetAttributeNotAllowed{$MOI.Name}:" *
          " Setting attribute $MOI.Name() cannot be performed: Message You may want" *
          " to use a `CachingOptimizer` in `AUTOMATIC` mode or you may need to call" *
          " `reset_optimizer` before doing this operation if the `CachingOptimizer`" *
          " is in `MANUAL` mode." ==
          "$MOI.SetAttributeNotAllowed{$MOI.Name}:" *
          " Setting attribute $MOI.Name() cannot be performed: Message You may want" *
          " to use a `CachingOptimizer` in `AUTOMATIC` mode or you may need to call" *
          " `reset_optimizer` before doing this operation if the `CachingOptimizer`" *
          " is in `MANUAL` mode."
    return
end

function test_errors_ResultIndexBoundsError()
    @test sprint(
        showerror,
        MOI.ResultIndexBoundsError(MOI.VariablePrimal(1), 0),
    ) ==
          "Result index of attribute MathOptInterface.VariablePrimal(1) out of" *
          " bounds. There are currently 0 solution(s) in the model."
end

function test_errors_InvalidCalbackUsage()
    @test sprint(
        showerror,
        MOI.InvalidCallbackUsage(MOI.LazyConstraintCallback(), MOI.UserCut(1)),
    ) ==
          "InvalidCallbackUsage: Cannot submit $(MOI.UserCut(1)) inside a MathOptInterface.LazyConstraintCallback()."
end

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$name", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

end

TestErrors.runtests()


text:
# ---
# title: 1015. Smallest Integer Divisible by K
# id: problem1015
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math
# link: <https://leetcode.com/problems/smallest-integer-divisible-by-k/description/>
# hidden: true
# ---
#
# Given a positive integer `K`, you need to find the **length** of the
# **smallest** positive integer `N` such that `N` is divisible by `K`, and `N`
# only contains the digit `1`.
#
# Return _the **length** of _`N`. If there is no such `N`, return -1.
#
# **Note:** `N` may not fit in a 64-bit signed integer.
#
#
#
# **Example 1:**
#
#
#
#     Input: K = 1
#     Output: 1
#     Explanation: The smallest answer is N = 1, which has length 1.
#
#
# **Example 2:**
#
#
#
#     Input: K = 2
#     Output: -1
#     Explanation: There is no such positive integer N divisible by 2.
#
#
# **Example 3:**
#
#
#
#     Input: K = 3
#     Output: 3
#     Explanation: The smallest answer is N = 111, which has length 3.
#
#
#
#
# **Constraints:**
#
#   * `1 <= K <= 105`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ArgParse

settings = ArgParseSettings()
@add_arg_table! settings begin
  "--N", "-N"
  help = "Number of sites"
  arg_type = Int
  default = 20
  "--Npart"
  help = "Number of particles"
  arg_type = Int
  default = 10
  "--t1"
  help = "Nearest neighbor hopping"
  arg_type = Float64
  default = 1.0
  "--t2"
  help = "Next-nearest neighbor hopping"
  arg_type = Float64
  default = 0.2
  "--U", "-U"
  help = "On-site potential"
  arg_type = Float64
  default = 1.0
  "--V1", "-V"
  help = "Nearest neighbor potential"
  arg_type = Float64
  default = 0.5
  "--nsweep", "-n"
  help = "Number of sweeps of DMRG"
  arg_type = Int64
  default = 6
  "--maxdim"
  help = "Maximum bond dimension in DMRG"
  nargs = '+'
  arg_type = Int64
  default = [50, 100, 200, 400, 800]
  "--mindim"
  help = "Minimum bond dimension in DMRG"
  nargs = '+'
  arg_type = Int64
  default = [10, 20]
  "--cutoff", "-c"
  help = "Truncation cutoff in DMRG"
  nargs = '+'
  arg_type = Float64
  default = [1e-12]
  "--noise"
  help = "Magnitude of noise to add in DMRG"
  nargs = '+'
  arg_type = Float64
  default = [1e-7, 1e-8, 1e-10, 0, 1e-11, 0]
  "input_file"
  help = "Input file"
  arg_type = String
end


text:
abstract type SKCE <: CalibrationErrorEstimator end

"""
    unsafe_skce_eval(k, p, y, p̃, ỹ)

Evaluate
```math
k((p, y), (p̃, ỹ)) - E_{z ∼ p}[k((p, z), (p̃, ỹ))] - E_{z̃ ∼ p̃}[k((p, y), (p̃, z̃))] + E_{z ∼ p, z̃ ∼ p̃}[k((p, z), (p̃, z̃))]
```
for kernel `k` and predictions `p` and `p̃` with corresponding targets `y` and `ỹ`.

This method assumes that `p`, `p̃`, `y`, and `ỹ` are valid and specified correctly, and
does not perform any checks.
"""
function unsafe_skce_eval end

# default implementation for classification
# we do not use the symmetry of `kernel` since it seems unlikely that `(p, y) == (p̃, ỹ)`
function unsafe_skce_eval(
    kernel::Kernel,
    p::AbstractVector{<:Real},
    y::Integer,
    p̃::AbstractVector{<:Real},
    ỹ::Integer,
)
    # precomputations
    n = length(p)

    @inbounds py = p[y]
    @inbounds p̃ỹ = p̃[ỹ]
    pym1 = py - 1
    p̃ỹm1 = p̃ỹ - 1

    tuple_p_y = (p, y)
    tuple_p̃_ỹ = (p̃, ỹ)

    # i = y, j = ỹ
    result = kernel((p, y), (p̃, ỹ)) * (1 - py - p̃ỹ + py * p̃ỹ)

    # i < y
    for i in 1:(y - 1)
        @inbounds pi = p[i]
        tuple_p_i = (p, i)

        # j < ỹ
        @inbounds for j in 1:(ỹ - 1)
            result += kernel(tuple_p_i, (p̃, j)) * pi * p̃[j]
        end

        # j = ỹ
        result += kernel(tuple_p_i, tuple_p̃_ỹ) * pi * p̃ỹm1

        # j > ỹ
        @inbounds for j in (ỹ + 1):n
            result += kernel(tuple_p_i, (p̃, j)) * pi * p̃[j]
        end
    end

    # i = y, j < ỹ
    @inbounds for j in 1:(ỹ - 1)
        result += kernel(tuple_p_y, (p̃, j)) * pym1 * p̃[j]
    end

    # i = y, j > ỹ
    @inbounds for j in (ỹ + 1):n
        result += kernel(tuple_p_y, (p̃, j)) * pym1 * p̃[j]
    end

    # i > y
    for i in (y + 1):n
        @inbounds pi = p[i]
        tuple_p_i = (p, i)

        # j < ỹ
        @inbounds for j in 1:(ỹ - 1)
            result += kernel(tuple_p_i, (p̃, j)) * pi * p̃[j]
        end

        # j = ỹ
        result += kernel(tuple_p_i, tuple_p̃_ỹ) * pi * p̃ỹm1

        # j > ỹ
        @inbounds for j in (ỹ + 1):n
            result += kernel(tuple_p_i, (p̃, j)) * pi * p̃[j]
        end
    end

    return result
end

# for binary classification with probabilities (corresponding to parameters of Bernoulli
# distributions) and boolean targets the expression simplifies to
# ```math
# k((p, y), (p̃, ỹ)) = (y(1-p) + (1-y)p)(ỹ(1-p̃) + (1-ỹ)p̃)(k((p, y), (p̃, ỹ)) - k((p, 1-y), (p̃, ỹ)) - k((p, y), (p̃, 1-ỹ)) + k((p, 1-y), (p̃, 1-ỹ)))
# ```
function unsafe_skce_eval(kernel::Kernel, p::Real, y::Bool, p̃::Real, ỹ::Bool)
    noty = !y
    notỹ = !ỹ
    z =
        kernel((p, y), (p̃, ỹ)) - kernel((p, noty), (p̃, ỹ)) -
        kernel((p, y), (p̃, notỹ)) + kernel((p, noty), (p̃, notỹ))
    return (y ? 1 - p : p) * (ỹ ? 1 - p̃ : p̃) * z
end

# evaluation for tensor product kernels
function unsafe_skce_eval(kernel::KernelTensorProduct, p, y, p̃, ỹ)
    κpredictions, κtargets = kernel.kernels
    return κpredictions(p, p̃) * unsafe_skce_eval_targets(κtargets, p, y, p̃, ỹ)
end

# resolve method ambiguity
function unsafe_skce_eval(
    kernel::KernelTensorProduct,
    p::AbstractVector{<:Real},
    y::Integer,
    p̃::AbstractVector{<:Real},
    ỹ::Integer,
)
    κpredictions, κtargets = kernel.kernels
    return κpredictions(p, p̃) * unsafe_skce_eval_targets(κtargets, p, y, p̃, ỹ)
end
function unsafe_skce_eval(kernel::KernelTensorProduct, p::Real, y::Bool, p̃::Real, ỹ::Bool)
    κpredictions, κtargets = kernel.kernels
    return κpredictions(p, p̃) * unsafe_skce_eval_targets(κtargets, p, y, p̃, ỹ)
end

function unsafe_skce_eval_targets(
    κtargets::Kernel,
    p::AbstractVector{<:Real},
    y::Integer,
    p̃::AbstractVector{<:Real},
    ỹ::Integer,
)
    # ensure that y ≤ ỹ (simplifies the implementation)
    y > ỹ && return unsafe_skce_eval_targets(κtargets, p̃, ỹ, p, y)

    # precomputations
    n = length(p)

    @inbounds begin
        py = p[y]
        pỹ = p[ỹ]
        p̃y = p̃[y]
        p̃ỹ = p̃[ỹ]
    end
    pym1 = py - 1
    pỹm1 = pỹ - 1
    p̃ym1 = p̃y - 1
    p̃ỹm1 = p̃ỹ - 1

    # i = y, j = ỹ
    result = κtargets(y, ỹ) * (1 - py - p̃ỹ + py * p̃ỹ)

    # i < y
    for i in 1:(y - 1)
        @inbounds pi = p[i]
        @inbounds p̃i = p̃[i]

        # i = j < y ≤ ỹ
        result += κtargets(i, i) * pi * p̃i

        # i < j < y ≤ ỹ
        @inbounds for j in (i + 1):(y - 1)
            result += κtargets(i, j) * (pi * p̃[j] + p[j] * p̃i)
        end

        # i < y < j < ỹ
        @inbounds for j in (y + 1):(ỹ - 1)
            result += κtargets(i, j) * (pi * p̃[j] + p[j] * p̃i)
        end

        # i < y ≤ ỹ < j
        @inbounds for j in (ỹ + 1):n
            result += κtargets(i, j) * (pi * p̃[j] + p[j] * p̃i)
        end
    end

    # y < i < ỹ
    for i in (y + 1):(ỹ - 1)
        @inbounds pi = p[i]
        @inbounds p̃i = p̃[i]

        # y < i = j < ỹ
        result += κtargets(i, i) * pi * p̃i

        # y < i < j < ỹ
        @inbounds for j in (i + 1):(ỹ - 1)
            result += κtargets(i, j) * (pi * p̃[j] + p[j] * p̃i)
        end

        # y < i < ỹ < j
        @inbounds for j in (ỹ + 1):n
            result += κtargets(i, j) * (pi * p̃[j] + p[j] * p̃i)
        end
    end

    # ỹ < i
    for i in (ỹ + 1):n
        @inbounds pi = p[i]
        @inbounds p̃i = p̃[i]

        # ỹ < i = j
        result += κtargets(i, i) * pi * p̃i

        # ỹ < i < j
        @inbounds for j in (i + 1):n
            result += κtargets(i, j) * (pi * p̃[j] + p[j] * p̃i)
        end
    end

    # handle special case y = ỹ
    if y == ỹ
        # i < y = ỹ, j = y = ỹ
        @inbounds for i in 1:(y - 1)
            result += κtargets(i, y) * (p[i] * p̃ym1 + pym1 * p̃[i])
        end

        # i = y = ỹ, j > y = ỹ
        @inbounds for j in (y + 1):n
            result += κtargets(y, j) * (pym1 * p̃[j] + p[j] * p̃ym1)
        end
    else
        # i < y
        for i in 1:(y - 1)
            @inbounds pi = p[i]
            @inbounds p̃i = p̃[i]

            # j = y < ỹ
            result += κtargets(i, y) * (pi * p̃y + pym1 * p̃i)

            # y < j = ỹ
            result += κtargets(i, ỹ) * (pi * p̃ỹm1 + pỹ * p̃i)
        end

        # i = y = j < ỹ
        result += κtargets(y, y) * pym1 * p̃y

        # i = y < j < ỹ and y < i < j = ỹ
        for ij in (y + 1):(ỹ - 1)
            @inbounds pij = p[ij]
            @inbounds p̃ij = p̃[ij]

            # i = y < j < ỹ
            result += κtargets(y, ij) * (pym1 * p̃ij + pij * p̃y)

            # y < i < j = ỹ
            result += κtargets(ij, ỹ) * (pij * p̃ỹm1 + pỹ * p̃ij)
        end

        # i = ỹ = j
        result += κtargets(ỹ, ỹ) * pỹ * (p̃ỹ - 1)

        # i = y < ỹ < j and i = ỹ < j
        for j in (ỹ + 1):n
            @inbounds pj = p[j]
            @inbounds p̃j = p̃[j]

            # i = y < ỹ < j
            result += κtargets(y, j) * (pym1 * p̃j + pj * p̃y)

            # i = ỹ < j
            result += κtargets(ỹ, j) * (p̃ỹm1 * pj + p̃j * pỹ)
        end
    end

    return result
end

function unsafe_skce_eval_targets(
    ::WhiteKernel,
    p::AbstractVector{<:Real},
    y::Integer,
    p̃::AbstractVector{<:Real},
    ỹ::Integer,
)
    @inbounds res = (y == ỹ) - p[ỹ] - p̃[y] + dot(p, p̃)
    return res
end

function unsafe_skce_eval_targets(κtargets::Kernel, p::Real, y::Bool, p̃::Real, ỹ::Bool)
    noty = !y
    notỹ = !ỹ
    z = κtargets(y, ỹ) - κtargets(noty, ỹ) - κtargets(y, notỹ) + κtargets(noty, notỹ)
    return (y ? 1 - p : p) * (ỹ ? 1 - p̃ : p̃) * z
end
function unsafe_skce_eval_targets(::WhiteKernel, p::Real, y::Bool, p̃::Real, ỹ::Bool)
    return 2 * (y - p) * (ỹ - p̃)
end


text:
function operator_pass(x, state)
    if x isa CSTParser.BinaryOpCall || x isa CSTParser.BinarySyntaxOpCall
        if CSTParser.precedence(x.op) in (8,13,14,16) || x.op.fullspan == 0
            ensure_no_space_after(x.arg1, state, state.offset)
            ensure_no_space_after(x.op, state, state.offset + x.arg1.fullspan)
        else
            ensure_single_space_after(x.arg1, state, state.offset)
            ensure_single_space_after(x.op, state, state.offset + x.arg1.fullspan)
        end
    elseif x isa CSTParser.WhereOpCall
        ensure_single_space_after(x.op, state, state.offset + x.arg1.fullspan)
        n = length(x.args)
        offset = state.offset + x.arg1.fullspan + x.op.fullspan
        for (i,a) in enumerate(x.args)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x isa CSTParser.EXPR{CSTParser.ColonOpCall}
        offset = state.offset
        n = length(x.args)
        for (i,a) in enumerate(x.args)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x isa CSTParser.EXPR{CSTParser.ChainOpCall} || x isa CSTParser.EXPR{CSTParser.Comparison}
        offset = state.offset
        n = length(x.args)
        for (i,a) in enumerate(x.args)
            if i != n
                ensure_single_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function tuple_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.TupleH}
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            i == n && continue
            if a isa CSTParser.PUNCTUATION && a.kind == Tokens.COMMA && !(x.args[i+1] isa CSTParser.PUNCTUATION)
                ensure_single_space_after(a, state, offset)
            elseif !(x.args[i + 1] isa CSTParser.EXPR{CSTParser.Parameters})
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function curly_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.Curly}
        offset = state.offset
        n = length(x)
        for (i, a) in enumerate(x)
            if i != n
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    end
end

function call_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.Call}
        offset = state.offset + x.args[1].fullspan
        n = length(x)
        for (i, a) in enumerate(x)
            i == 1 && continue
            if a isa CSTParser.PUNCTUATION && a.kind == Tokens.COMMA
                ensure_single_space_after(a, state, offset)
            # elseif a isa CSTParser.EXPR{CSTParser.Parameters}
            elseif i != n && !(x.args[i + 1] isa CSTParser.EXPR{CSTParser.Parameters})
                ensure_no_space_after(a, state, offset)
            end
            offset += a.fullspan
        end
    elseif x isa CSTParser.EXPR{CSTParser.Kw}
        ensure_single_space_after(x.args[1], state, state.offset)
        ensure_single_space_after(x.args[2], state, state.offset + x.args[1].fullspan)
    end
end

function forloop_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.For}
        offset = state.offset + x.args[1].fullspan
        for a in x.args[2]
            # convert iter = I into iter in I
            if a isa CSTParser.BinarySyntaxOpCall && CSTParser.is_eq(a.op)
                offset += a.arg1.fullspan
                push!(state.edits, Edit(offset+1:offset+2, "in "))
                offset += a.op.fullspan
                offset += a.arg2.fullspan
            else
                offset += a.fullspan
            end
        end
    end
end

# TODO: move this to CSTParser?
function str_value(x::CSTParser.PUNCTUATION)
    x.kind == Tokens.LPAREN && return "("
    x.kind == Tokens.LBRACE && return "{"
    x.kind == Tokens.LSQUARE && return "["
    x.kind == Tokens.RPAREN && return ")"
    x.kind == Tokens.RBRACE && return "}"
    x.kind == Tokens.RSQUARE && return "]"
    x.kind == Tokens.COMMA && return ","
    x.kind == Tokens.SEMICOLON && return ";"
    x.kind == Tokens.AT_SIGN && return "@"
    x.kind == Tokens.DOT && return "."
    return ""
end

function str_value(x::CSTParser.EXPR)
    s = ""
    for a in x
        s *= str_value(a)
    end
    return s
end

function str_value(x::CSTParser.UnarySyntaxOpCall)
    s = str_value(x.arg1)
    s *= str_value(x.arg2)
    return s
end
str_value(x) = CSTParser.str_value(x)

function doc_pass(x, state)
    if x isa CSTParser.EXPR{CSTParser.MacroCall} && x.args[1] isa CSTParser.EXPR{CSTParser.GlobalRefDoc}
        # Align global docstring to:
        #
        # """
        # doc
        # """
        #
        # If the doc is single quoted i.e. "doc", they will be replaced with triple quotes.
        offset = state.offset + x.args[1].fullspan
        doc = x.args[2]

        val = str_value(doc)

        s = strip(val, ['\n'])
        ds = string("\"\"\"\n", s, "\n", "\"\"\"\n")

        # Check if docstring needs to be edited
        if length(ds) != doc.fullspan || s != val
            # Remove previous docstring
            push!(state.edits, Edit(offset+1:offset+doc.fullspan, ""))
            # Append newly formatted docstring
            push!(state.edits, Edit(offset, ds))
        end
    end
end


text:
# This file was generated, do not modify it. # hide
using MLJ, RDatasets

auto = dataset("ISLR", "Auto")
y, X = unpack(auto, ==(:MPG), col->true)
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=444);


text:
@testset "BigTIFFLong8.tif" begin
filename = joinpath(testbase, "BigTIFFLong8.tif")
    # $ ScanImageTiffReader image shape BigTIFFLong8.tif
    # Shape: 3 x 64 x 64 @ u8
    tsize = ScanImageTiffReader.open(filename) do io
        size(io)
    end
    @test tsize == (3, 64, 64)
    ttype = ScanImageTiffReader.open(filename) do io
        pxtype(io)
    end
    @test ttype == UInt8

    # $ ScanImageTiffReader image bytes BigTIFFLong8.tif
    # 12 kB
    dat = ScanImageTiffReader.open(filename) do io
        data(io)
    end
    @test abs(sizeof(dat)/2^10 - 12.0) < 1e-5
    # @test ScanImageTiffReader.open(filename, length) == size(dat, 3) # TODO: check this assumption

    # $ ScanImageTiffReader descriptions --frame 0 BigTIFFLong8.tif
    #
    @test ScanImageTiffReader.open(description, filename, 1) == ""

    # $ ScanImageTiffReader metadata BigTIFFLong8.tif
    #
    @test ScanImageTiffReader.open(metadata, filename) == ""
end


text:
@testset "test Utils" begin
    include("base.jl")
    include("tiling.jl")
    include("sum_tree.jl")
    include("readers_writer_lock.jl")
end


text:
ref_netcount(A, j, j′) = length(union(A.rowval[A.colptr[j]:A.colptr[j′] - 1]))

ref_partwisenetcount(A, Φ_domain, j, j′, k) = length(intersect(union(A[:, j:(j′ - 1)].rowval), Φ_domain.prm[Φ_domain.spl[k]:Φ_domain.spl[k + 1] - 1]))
ref_partwisevertexcount(At, Φ_domain, j, j′, k) = length(intersect(union(At[:, Φ_domain.prm[Φ_domain.spl[k]:Φ_domain.spl[k + 1] - 1]].rowval), j:(j′ - 1)))

ref_selfnetcount(A, j, j′) = length(setdiff(union(A.rowval[A.colptr[j]:A.colptr[j′] - 1]), union(A.rowval[A.colptr[1]:A.colptr[j] - 1], A.rowval[A.colptr[j′]:A.colptr[end] - 1])))

#ref_primarypincount(A, Φ_domain, j, j′, k) = sum(A[Φ_domain.prm[Φ_domain.spl[k]:Φ_domain.spl[k + 1] - 1], j:j′ - 1] .!= 0)

ref_selfpincount(A, j, j′) = sum(A[j:j′ - 1, j:j′ - 1] .!= 0)
ref_dianetcount(A, j, j′) = ref_netcount(A + UniformScaling(true), j, j′)

@testset "SparseColorArrays" begin

    for hint = ((), (ChainPartitioners.RandomHint(),), (ChainPartitioners.SparseHint(),), (ChainPartitioners.StepHint(),))
        for m = test_dims
            for n = test_dims
                A = dropzeros!(sprand(UInt, m, n, 0.5))
                net = netcount(hint..., A)
                selfnet = selfnetcount(hint..., A)
                @testset "hint=$hint, n = $n, m = $m" begin
                    for (j, j′) in test_points(10, 1:n + 1, 1:n + 1)
                        (j, j′) = minmax(j, j′)
                        @test net[j, j′] == ref_netcount(A, j, j′)
                        @test selfnet[j, j′] == ref_selfnetcount(A, j, j′)

                        for _ = 1:10
                            _j, _j′= rand(filter(((_j, _j′),) -> destep(_j) < destep(_j′), collect(Iterators.product(Δ(j, 1:n + 1), Δ(j′, 1:n + 1)))))
                            (j, j′) = map(destep, (_j, _j′))
                            @test ChainPartitioners.Step(net)(_j, _j′) == ref_netcount(A, j, j′)
                            @test ChainPartitioners.Step(selfnet)(_j, _j′) == ref_selfnetcount(A, j, j′)
                        end
                    end
                end
            end
        end

        for m = test_dims, K = 1:8
            for n = test_dims
                Φ_map = MapPartition(K, rand(1:K, m))
                Φ_domain = convert(DomainPartition, Φ_map)
                n = m
                A = dropzeros!(sprand(UInt, m, n, 0.5))
                At = SparseMatrixCSC(A')
                if nnz(A) > 0
                    args, Ap = partwise(A, Φ_map)
                    partwisenet = partwisecount!(hint..., args..., netcount(hint..., Ap))
                    partwisevertex = partwisecount!(hint..., args..., VertexCount())
                    @testset "hint=$hint, m = $m, n = $n, K = $K" begin
                        for (j, j′, k) in test_points(10, 1:n + 1, 1:n + 1, 1:K)
                            (j, j′) = minmax(j, j′)

                            @test partwisenet[j, j′, k] == ref_partwisenetcount(A, Φ_domain, j, j′, k)
                            @test partwisevertex[j, j′, k] == ref_partwisevertexcount(At, Φ_domain, j, j′, k)

                            for _ = 1:10
                                _j, _j′= rand(filter(((_j, _j′),) -> destep(_j) < destep(_j′), collect(Iterators.product(Δ(j, 1:n + 1), Δ(j′, 1:n + 1)))))
                                _k = rand(Δ(k, 1:K))
                                (j, j′, k) = map(destep, (_j, _j′, _k))
                                @test ChainPartitioners.Step(partwisenet)(_j, _j′, _k) == ref_partwisenetcount(A, Φ_domain, j, j′, k)
                                @test ChainPartitioners.Step(partwisevertex)(_j, _j′, _k) == ref_partwisevertexcount(At, Φ_domain, j, j′, k)
                            end
                        end
                    end
                end
            end
        end

        for m = test_dims
            n = m
            A = dropzeros!(sprand(UInt, m, n, 0.5))
            if nnz(A) > 0
                selfpin = selfpincount(hint..., A)
                dianet = dianetcount(hint..., A)
                @testset "hint=$hint, m = $m, n = $n" begin
                    for (j, j′) in test_points(10, 1:n + 1, 1:n + 1)
                        (j, j′) = minmax(j, j′)

                        @test selfpin[j, j′] == ref_selfpincount(A, j, j′)
                        @test dianet[j, j′] == ref_dianetcount(A, j, j′)

                        for _ = 1:10
                            _j, _j′= rand(filter(((_j, _j′),) -> destep(_j) < destep(_j′), collect(Iterators.product(Δ(j, 1:n + 1), Δ(j′, 1:n + 1)))))
                            (j, j′) = map(destep, (_j, _j′))
                            @test ChainPartitioners.Step(selfpin)(_j, _j′) == ref_selfpincount(A, j, j′)
                            @test ChainPartitioners.Step(dianet)(_j, _j′) == ref_dianetcount(A, j, j′)
                        end
                    end
                end
            end
        end
    end
end


text:
"""
    JuMP.build_variable(::Function, info::JuMP.VariableInfo, set::T) where T <: MOI.AbstractScalarSet

DOCSTRING

# Arguments:
- ``: DESCRIPTION
- `info`: DESCRIPTION
- `set`: DESCRIPTION
"""
function JuMP.build_variable(
    ::Function,
    info::JuMP.VariableInfo,
    set::T,
) where {T<:MOI.AbstractScalarSet}
    return JuMP.VariableConstrainedOnCreation(JuMP.ScalarVariable(info), set)
end

"""
    Optimizer <: MOI.AbstractOptimizer

DOCSTRING

# Arguments:
- `solver::Solver`: DESCRIPTION
- `status::MOI.TerminationStatusCode`: DESCRIPTION
- `options::Options`: DESCRIPTION
"""
mutable struct Optimizer <: MOI.AbstractOptimizer
    solver::LS.MainSolver
    int_vars::Set{Int}
    compare_vars::Set{Int}
end

"""
    Optimizer(model = Model(); options = Options())

DOCSTRING
"""
function Optimizer(model = model(); options = Options())
    return Optimizer(
        solver(model, options = options),
        Set{Int}(),
        Set{Int}()
        )
end

# forward functions from Solver
@forward Optimizer.solver LS.variable!, LS._set_domain!, LS.constraint!, LS.solution
@forward Optimizer.solver LS.max_domains_size, LS.objective!, Base.empty!, LS._inc_cons!
@forward Optimizer.solver LS._best_bound, LS.best_value, LS.is_sat, LS.get_value
@forward Optimizer.solver LS.domain_size, LS.best_values, LS._max_cons, LS.update_domain!
@forward Optimizer.solver LS.get_variable, LS.has_solution, LS.sense, LS.sense!
@forward Optimizer.solver LS.time_info, LS.status

# forward functions from Solver (from Options)
@forward Optimizer.solver LS._verbose, LS.set_option!, LS.get_option

"""
    MOI.get(::Optimizer, ::MOI.SolverName) = begin

DOCSTRING
"""
MOI.get(::Optimizer, ::MOI.SolverName) = "LocalSearchSolvers"

"""
    MOI.set(::Optimizer, ::MOI.Silent, bool = true) = begin

DOCSTRING

# Arguments:
- ``: DESCRIPTION
- ``: DESCRIPTION
- `bool`: DESCRIPTION
"""
MOI.set(::Optimizer, ::MOI.Silent, bool = true) = @debug "TODO: Silent"

"""
    MOI.is_empty(model::Optimizer) = begin

DOCSTRING
"""
MOI.is_empty(model::Optimizer) = LS._is_empty(model.solver)

"""
Copy constructor for the optimizer
"""
MOIU.supports_default_copy_to(::Optimizer, copy_names::Bool) = false
function MOI.copy_to(model::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(model, src; kws...)
end

"""
    MOI.optimize!(model::Optimizer)
"""
MOI.optimize!(optimizer::Optimizer) = solve!(optimizer.solver)

"""
    DiscreteSet(values)
"""
struct DiscreteSet{T <: Number} <: MOI.AbstractScalarSet
    values::Vector{T}
end
DiscreteSet(values) = DiscreteSet(collect(values))
DiscreteSet(values::T...) where {T<:Number} = DiscreteSet(collect(values))

"""
    Base.copy(set::DiscreteSet) = begin

DOCSTRING
"""
Base.copy(set::DiscreteSet) = DiscreteSet(copy(set.values))

"""
    MOI.empty!(opt) = begin

DOCSTRING
"""
MOI.empty!(opt) = empty!(opt)

function MOI.is_valid(optimizer::Optimizer, index::CI{SVF, MOI.Integer})
    return index.value ∈ optimizer.int_vars
end


text:
export test_problem_ncqp, tridiagonal

#
# test the linear program solver on a simple problem
#
using Ipopt, KNITRO
#simple_tests();

function test_problem_ncqp(name::String)
	file_name = name * ".mat";

  dir = dirname(@__FILE__) * "/Problems/"

  A, b, c = get_netlib_problem(dir, file_name);
  Q = tridiagonal(length(c), 0.0,  1.0);

	println("Solving ", file_name, " with the homogeneous algorithm")
	println(size(A,2), " variables and ", size(A,1), " constraints")
	println("Non-zeros: ", length(nonzeros(A)))
	settings = class_settings();

	#println("=================== Linear system solver is matlab ldl ==================")
	println("=================== Linear system solver is MUMPS ==================")
	settings.linear_system_solver = linear_solver_MUMPS();#linear_solver_MATLAB();

	settings.newton_solver = class_homogeneous_newton();

  qp = class_quadratic_program(A, b, c, Q);
  homogeneous_algorithm(qp,settings)
  #trivial_test(A, b, c, Q, 1, file_name, settings, true);

	#settings.newton_solver = class_newton_ip();
	#trivial_test(A, b, c, Q, 1, file_name, settings, true);

	println("=========================================================================")
	println("Calls IPOPT")
	solve_with_JuMP(A, b, c, Q, IpoptSolver(max_iter=1000,print_level=3));
end

max_forks_repo_name: ohinder/Everest


text:
function InitialRTBP(T=Float64)

   μ = parse(BigFloat,"0.012277471")

   q=[parse(BigFloat,"0.994"),parse(BigFloat,"0")]
   p=[parse(BigFloat,"0"),parse(BigFloat,"-2.00158510637908252240537862224")+0.994]

   u0=convert.(T,[q;p])

   return u0,convert(T,μ)

end


text:
function test_leaky()
    return leaky_sigm(1000) > 1.0 &&
           leaky_relu(-10) < 0.0 &&
           leaky_tanh(1000) > 1.0
end


text:
using Uchiyama
using Documenter

using Literate
using Plots

# generate examples
examples = ["event_driven.jl", "hard_spheres.jl",
            "periodic_hard_spheres.jl", "periodic_uchiyama.jl"]
output = joinpath(@__DIR__, "src")
examples_dir = joinpath(@__DIR__, "..", "examples")
for example in examples
    jl_file = joinpath(examples_dir, example)
    Literate.markdown(jl_file, output, documenter=true)
    Literate.notebook(jl_file, output, execute=false)
    #Literate.script(EXAMPLE, OUTPUT)
end

makedocs(;
    modules=[Uchiyama],
    authors="Nathalie Ayi and Pierre Navaro",
    repo="https://github.com/pnavaro/Uchiyama.jl/blob/{commit}{path}#L{line}",
    sitename="Uchiyama.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pnavaro.github.io/Uchiyama.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Uchiyama model in box" => "event_driven.md",
        "Uchiyama model periodic" => "periodic_uchiyama.md",
        "Hard spheres in box" => "hard_spheres.md",
        "Hard spheres periodic" => "periodic_hard_spheres.md",
    ],
)

deploydocs(;
    repo="github.com/pnavaro/Uchiyama.jl",
)


text:
#name: JuKNNApply
#meta.mlname: JuKNN
#meta.mlrole: apply
#description: Custom Julia Apply func for KNN using ScikitLearn
#language: julia
#input: blob model
#input: dataframe df
#input: string namesKeys [original features' names]
#input: string namesValues [new features' names]
#output: dataframe data_out

#>>> Load necessary packages
#using Pkg
#Pkg.add("PyCall")
import PyCall
import BSON: @load
using ScikitLearn

#>>> If original(namesKeys) and new(namesValues) passed, map original names to new
namesValues=split(namesValues, ",")
namesKeys=split(namesKeys, ",")
if length(namesKeys) > 0
    featuresNames = names(df)
    for i in 1:length(namesKeys)
        replace!(featuresNames, namesValues[i]  =>  namesKeys[i])
    end
end
testX=Matrix(select(df, Symbol.(featuresNames)))

#>>> Retrieve saved/trained model
@load file=model trained_model

#>>> Predict using trained model
predY = predict(trained_model, testX)
data_out = DataFrame([predY],[:out])

max_stars_repo_name: DotNetPart/public

max_issues_repo_name: DotNetPart/public

max_forks_repo_name: DotNetPart/public


text:
using OrdinaryDiffEq, ForwardDiff, LinearAlgebra, Test

const nknots = 10
const h = 1.0/(nknots+1)
x = range(0, step=h, length=nknots)
u0 = sin.(π*x)

@inline function f(du,u,p,t)
  du .= zero(eltype(u))
  u₃ = @view u[3:end]
  u₂ = @view u[2:end-1]
  u₁ = @view u[1:end-2]
  @. du[2:end-1] = p[1]*((u₃ - 2*u₂ + u₁)/(h^2.0))
  nothing
end

p_true = [0.42]
jac_proto = Tridiagonal(similar(u0,nknots-1), similar(u0), similar(u0, nknots-1))
prob = ODEProblem(ODEFunction(f,jac_prototype=jac_proto), u0, (0.0,1.0), p_true)
sol_true = solve(prob, Rodas4P(), saveat=0.1)

function loss(p)
  _prob = remake(prob, p=p)
  sol = solve(_prob, Rodas4P(autodiff=false), saveat=0.1)
  sum((sol .- sol_true).^2)
end
@test ForwardDiff.gradient(loss, [1.0])[1] ≈ 0.6662949361011025


text:
cd(@__DIR__)
using Pkg
Pkg.activate(@__DIR__); Pkg.instantiate()

include("hodgkin_huxley_setup.jl")
include("helpers.jl")

system = neuron.sys
@named simulation = ODESystem([D(system.Isyn) ~ 0]; systems = [system])
simplified = structural_simplify(simulation)

sim_exp = ModelingToolkit.build_torn_function(simplified; expression=true)
cleaned = clean_expr(sim_exp, simplified)
#display_expr(cleaned)
#clipboard(format_text(string(prettify(cleaned)))) # can paste readable output below

############################################################################################
# Generated function expression from MTK lightly cleaned up for better readability

function hh_generated!(du, u, p, t)
    let neuron₊Isyn = @inbounds(u[1]),
        neuron₊Vₘ = @inbounds(u[2]),
        neuron₊NaV₊m = @inbounds(u[3]),
        neuron₊NaV₊h = @inbounds(u[4]),
        neuron₊Kdr₊n = @inbounds(u[5]),
        neuron₊cₘ = @inbounds(p[1]),
        neuron₊aₘ = @inbounds(p[2]),
        neuron₊ENa = @inbounds(p[3]),
        neuron₊EK = @inbounds(p[4]),
        neuron₊El = @inbounds(p[5]),
        neuron₊NaV₊gbar = @inbounds(p[6]),
        neuron₊Kdr₊gbar = @inbounds(p[7]),
        neuron₊leak₊g = @inbounds(p[8])

        let
            du[1] = 0
            du[2] =
                (
                    if t > 100.0
                        if t < 200.0
                            0.0004
                        else
                            0.0
                        end
                    else
                        0.0
                    end +
                    -1neuron₊Isyn +
                    -1 * neuron₊aₘ * neuron₊leak₊g * (neuron₊Vₘ + -1neuron₊El) +
                    -1 *
                    neuron₊Kdr₊gbar *
                    neuron₊aₘ *
                    neuron₊Kdr₊n^4.0 *
                    (neuron₊Vₘ + -1neuron₊EK) +
                    -1 *
                    neuron₊NaV₊gbar *
                    neuron₊aₘ *
                    neuron₊NaV₊h *
                    neuron₊NaV₊m^3.0 *
                    (neuron₊Vₘ + -1neuron₊ENa)
                ) *
                inv(neuron₊aₘ) *
                inv(neuron₊cₘ)
            du[3] =
                (1 + -1neuron₊NaV₊m) * if neuron₊Vₘ == -40.0
                    1.0
                else
                    (4.0 + 0.1neuron₊Vₘ) * inv(1.0 + -1 * exp(-4.0 + -0.1neuron₊Vₘ))
                end +
                -4.0 *
                neuron₊NaV₊m *
                exp(-3.6111111111111107 + -0.05555555555555555neuron₊Vₘ)
            du[4] =
                -1 * neuron₊NaV₊h * inv(1.0 + exp(-3.5 + -0.1neuron₊Vₘ)) +
                0.07 * exp(-3.25 + -0.05neuron₊Vₘ) * (1 + -1neuron₊NaV₊h)
            du[5] =
                if neuron₊Vₘ == -55.0
                    0.1
                else
                    (0.55 + 0.01neuron₊Vₘ) * inv(1.0 + -1 * exp(-5.5 + -0.1neuron₊Vₘ))
                end * (1 + -1neuron₊Kdr₊n) +
                -0.125 * neuron₊Kdr₊n * exp(-0.8125 + -0.0125neuron₊Vₘ)
            nothing
        end
    end
end


text:
module zad2

function zad22(funk,xL,xU;ep=1e-8,golden=0.618)
              fa=funk(xL)
              fb=funk(xU)
              x1=xU-((xU-xL)*(golden))
              x2=xL+((xU-xL)*(golden))
              err=Inf
              while err>ep
                  if funk(x2)<funk(x1)
                      xL=x1
                      x1=x2
                      x2=xU-((xU-xL)*(golden))
                  elseif funk(x2)>funk(x1)
                      xU=x2
                      x2=x1
                  x1=xL+((xU-xL)*(golden))
                  else
                      xL=(x1+x2)/2
                      xU=xL
                  end
                 err=2*abs((xU-xL)/(xU+xL))
                  end
              Xmin=(x1+x2)/2
              print("Minimum: ",round(funk(Xmin),digits=2), " for: ", Xmin, " x1: " ,round(x1,digits=2), " x2: ", round(x2,digits=2))

       return(round(funk(Xmin),digits=2),round(Xmin,digits=2))
       end

export zad22

end # module


text:
function write_ascii(A::Array{Float64, 2}, filename::String, ascii_header::Dict)
    f = open(filename, "w")

    write(f, "ncols         $(ascii_header["ncols"])\n")
    write(f, "nrows         $(ascii_header["nrows"])\n")
    write(f, "xllcorner     $(ascii_header["xllcorner"])\n")
    write(f, "yllcorner     $(ascii_header["yllcorner"])\n")
    write(f, "cellsize      $(ascii_header["cellsize"])\n")
    write(f, "nodata_value  $(ascii_header["nodata_value"])\n")

    writedlm(f, round.(A, digits=8), ' ')
    close(f)
end

function parse_ascii_header(path::String)
    cf = init_ascii_header()
    f = open(path, "r")
    for i = 1:6
        key_val = split(readline(f))
        key_val[1] = lowercase(key_val[1])
        cf["$(key_val[1])"] = key_val[2]
    end

    cf
end

function read_ascii(path::String)
    a = readdlm(path, Float64; skipstart = 6)
    a
end


text:
#White Noise kernel

"""
    Noise <: Kernel

Noise kernel (covariance)
```math
k(x,x') = σ²δ(x-x'),
```
where ``δ`` is the Kronecker delta function and ``σ`` is the signal standard deviation.
"""
mutable struct Noise <: Kernel
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Noise(lσ::Float64)

    Create `Noise` with signal standard deviation `exp(lσ)`.
    """
    Noise(lσ::Float64) = new(exp(2 * lσ), [])
end

Statistics.cov(noise::Noise, sameloc::Bool) = sameloc ? noise.σ2 : 0.0
Statistics.cov(noise::Noise, x::VecF64, y::VecF64) = cov(noise, euclidean(x, y) < eps())

get_params(noise::Noise) = Float64[log(noise.σ2) / 2]
get_param_names(noise::Noise) = [:lσ]
num_params(noise::Noise) = 1

function set_params!(noise::Noise, hyp::VecF64)
    length(hyp) == 1 || throw(ArgumentError("Noise kernel only has one parameter"))
    noise.σ2 = exp(2 * hyp[1])
end

dKij_dθp(noise::Noise, X::MatF64, i::Int, j::Int, p::Int, dim::Int) =
    2 * cov(noise, view(X, :, i), view(X, :, j))

@inline function dKij_dθp(noise::Noise, X::MatF64, data::EmptyData, i::Int, j::Int, p::Int, dim::Int)
    return dKij_dθp(noise, X, i, j, p, dim)
end