text:
using DataFrames
using InteractiveUtils: subtypes
using InlineStrings: InlineStrings  # for `eval(repr(...))`
using PowerFlowData
using Tables
using Test

@testset "PowerFlowData.jl" begin
    @testset "Infers" begin
        for T in subtypes(PowerFlowData.Records)
            isabstracttype(T) && continue
            @inferred T()
            @inferred T(1)
        end
    end

    @testset "Tables interface" begin
        struct TestRecs <: PowerFlowData.Records
            x1::Vector{Int}
            yy::Vector{Float64}
        end
        recs = TestRecs([1, 2, 3], [3.14, 1.72, 0.05])

        @test Tables.istable(TestRecs)

        @test Tables.columnaccess(TestRecs)
        @test Tables.columns(recs) === recs
        @test Tables.columnnames(recs) == (:x1, :yy)

        @test recs.x1 == [1, 2, 3]
        @test Tables.getcolumn(recs, :x1) == [1, 2, 3]
        @test Tables.getcolumn(recs, 1) == [1, 2, 3]

        @test length(recs) == 3
        @test size(recs) == (3, 2)

        df = DataFrame(recs)
        @test df isa DataFrame
        @test size(df) == (3, 2)

        for T in (
            Buses, Loads, Generators, Branches, Transformers, AreaInterchanges,
            TwoTerminalDCLines, VSCDCLines, SwitchedShunts, ImpedanceCorrections,
            MultiTerminalDCLines, MultiSectionLineGroups, Zones, InterAreaTransfers,
            Owners, FACTSDevices
        )
            @test T <: PowerFlowData.Records
            @test Tables.istable(T)
        end

        caseid = CaseID(ic=0, sbase=100.0)
        @test caseid[1] == caseid[:ic]
        @test caseid[2] == caseid[:sbase]
        @test NamedTuple(caseid) isa NamedTuple
    end

    @testset "show" begin
        # test of file which has MultiTerminalDCLines as they have their own `show`.
        net = parse_network("testfiles/synthetic_data_v29.raw")

        # CaseID should have a parseable `repr`; `AbstractRows` don't get this for free.
        @test startswith(repr(net.caseid), "CaseID(ic = 0, sbase = 100.0, ")
        @test isequal(eval(Meta.parse(repr(net.caseid))), CaseID(ic=0, sbase=100.0))

        mime = MIME("text/plain")
        context = :compact => true
        @test repr(mime, net; context) == "Network"
        @test repr(mime, net) == strip(
            """
            Network (v30) with 17 data categories:
             $(sprint(show, mime, net.caseid))
             $(sprint(show, mime, net.buses; context))
             $(sprint(show, mime, net.loads; context))
             $(sprint(show, mime, net.generators; context))
             $(sprint(show, mime, net.branches; context))
             $(sprint(show, mime, net.transformers; context))
             $(sprint(show, mime, net.area_interchanges; context))
             $(sprint(show, mime, net.two_terminal_dc; context))
             $(sprint(show, mime, net.vsc_dc; context))
             $(sprint(show, mime, net.switched_shunts; context))
             $(sprint(show, mime, net.impedance_corrections; context))
             $(sprint(show, mime, net.multi_terminal_dc; context))
             $(sprint(show, mime, net.multi_section_lines; context))
             $(sprint(show, mime, net.zones; context))
             $(sprint(show, mime, net.area_transfers; context))
             $(sprint(show, mime, net.owners; context))
             $(sprint(show, mime, net.facts; context))
            """
        )
        @test startswith(repr(mime, net.caseid), "CaseID: (ic = 0, sbase = 100.0, ")

        @test repr(mime, net.buses; context=(:compact => true)) == "Buses with 2 records"
        @test startswith(repr(mime, net.buses), "Buses with 2 records, 11 columns:\n──")

        mt_dc_line = net.multi_terminal_dc.lines[1]
        @test eval(Meta.parse(repr(mt_dc_line))) isa MultiTerminalDCLine
        @test repr(mime, mt_dc_line) == strip(
            """
            $(sprint(show, mime, mt_dc_line.line_id))
            $(sprint(show, mime, mt_dc_line.converters))
            $(sprint(show, mime, mt_dc_line.buses))
            $(sprint(show, mime, mt_dc_line.links))
            """
        )
        line_id = mt_dc_line.line_id
        @test eval(Meta.parse(repr(line_id))) isa DCLineID
    end

    @testset "v30 file" begin
        net1 = parse_network("testfiles/synthetic_data_v30.raw")
        @test net1 isa Network

        caseid = net1.caseid
        @test caseid.ic == 0
        @test caseid.sbase == 100.0

        # Test first and last columns are parsed as expected
        buses = net1.buses
        @test buses.i == [111, 112, 113]
        @test buses.owner == [1, 2, 2]

        # Test string column as expected
        @test buses.name[2] == "D2JK"

        loads = net1.loads
        @test loads.i == [111, 113]  # first col
        @test loads.owner == [1, 2]  # last col
        @test isequal(loads.intrpt, [missing, missing])

        gens = net1.generators
        @test gens.i == [111, -112, 113]
        @test gens.f1 == [1.0, 1.0, 1.0]  # last col guaranteed to exist
        @test isequal(gens.wpf, [missing, missing, missing])
        @test gens.id[1] == "ST"

        branches = net1.branches
        @test branches.i == [111, 111, 112]
        @test branches.j == [112, -113, 113]  # negative numbers should be allowed
        @test branches.f1 == [1.0, 1.0, 1.0]
        @test branches.ckt[1] == "3"

        transformers = net1.transformers
        @test transformers.i == [112, 113]                     #  1st entry of 1st row
        @test transformers.f1 == [1.0, 1.0]                    # last entry of 1st row
        @test transformers.r1_2 == [0.032470, 0.039570]        #  1st entry of 2nd row
        @test transformers.sbase1_2 == [200.0, 200.0]          # last entry of 2nd row (T2)
        @test isequal(transformers.anstar, [missing, 1.01893]) # last entry of 2nd row (T3)
        @test transformers.windv1 == [1.0, 1.0]                #  1st entry of 3rd row
        @test transformers.cx1 == [0.0, 0.0]                   # last entry of 3rd row
        @test transformers.windv2 == [1.0, 1.0]                #  1st entry of 4th row
        @test transformers.nomv2 == [169.0, 169.0]             # last entry of 4th row (T2)
        @test isequal(transformers.cx2, [missing, 0.0])        # last entry of 4th row (T3)
        @test isequal(transformers.windv3, [missing, 13.8])    #  1st entry of 5th row
        @test isequal(transformers.cx3, [missing, 0.0])        # last entry of 5th row
        @test transformers.ckt[1] == "G1"                      # string col

        # v30 testfile has both 2-winding and 3-winding data, so should return all columns
        @test length(Tables.columnnames(transformers)) == fieldcount(Transformers)
        @test size(DataFrame(transformers)) == (2, fieldcount(Transformers))
        @test size(transformers) == (2, fieldcount(Transformers))
        @test length(transformers) == 2

        area_interchanges = net1.area_interchanges
        @test area_interchanges.i == [113]
        @test area_interchanges.isw == [456]
        @test area_interchanges.pdes == [2121.7211]
        @test area_interchanges.ptol == [6.0]
        @test area_interchanges.arname == ["ABC"]

        two_terminal_dc = net1.two_terminal_dc
        @test two_terminal_dc.i == [11]           #  1st entry of 1st row
        @test two_terminal_dc.cccacc == [1.0]     # last entry of 1st row
        @test two_terminal_dc.ipr == [112]        #  1st entry of 2nd row
        @test two_terminal_dc.xcapr == [0.0]      # last entry of 2nd row
        @test two_terminal_dc.ipi == [2222]       #  1st entry of 3nd row
        @test two_terminal_dc.xcapi == [2.0]      # last entry of 3rd row

        vsc_dc = net1.vsc_dc
        @test vsc_dc.name == ["line 1"]   #  1st entry of 1st row
        @test vsc_dc.f4 == [1.0]           # last entry of 1st row
        @test vsc_dc.ibus1 == [1117]       #  1st entry of 2nd row
        @test vsc_dc.rmpct1 == [100.0]     # last entry of 2nd row
        @test vsc_dc.ibus2 == [114]        #  1st entry of 3nd row
        @test vsc_dc.rmpct2 == [100.0]     # last entry of 3nd row

        switched_shunts = net1.switched_shunts
        @test switched_shunts.i == [113]   # first col
        @test switched_shunts.n1 == [1]    # `n1` always present
        @test switched_shunts.b1 == [26.0] # `b1` always present
        @test switched_shunts.n8 == [0]    # `n8` not present; should default to zero
        @test switched_shunts.b8 == [0.0]  # last col; `b8` not present; default to zero

        impedance_corrections = net1.impedance_corrections
        @test impedance_corrections.i == [1]       # first col
        @test impedance_corrections.f11 == [0.0]   # last col; `f11` not present; default to zero

        multi_terminal_dc = net1.multi_terminal_dc
        @test isempty(multi_terminal_dc)

        multi_section_lines = net1.multi_section_lines
        @test multi_section_lines.i == [1, 114]                     # first col
        @test multi_section_lines.id == ["&1", "&2"]                # string col
        @test isequal(multi_section_lines.dum2, [missing, 5])       # `dum2` present only for 1 row of the data
        @test isequal(multi_section_lines.dum9, [missing, missing]) # last col; not present; default to missing

        zones = net1.zones
        @test zones.i == [117, 127, 227]
        @test zones.zoname == ["ABC", "CDEF", "CDEG"]

        area_transfers = net1.area_transfers
        @test area_transfers.arfrom == [1, 1]
        @test area_transfers.ptran == [10.0, -20.0]

        owners = net1.owners
        @test owners.i == [1, 2]
        @test owners.owname == ["ABC", "CDE"]

        facts = net1.facts
        @test facts.n == [1]
        @test facts.vsref == [0]
    end

    @testset "v29 file" begin
        net2 = parse_network("testfiles/synthetic_data_v29.raw")

        caseid = net2.caseid
        @test caseid.ic == 0
        @test caseid.sbase == 100.0

        # Test first and last columns are parsed as expected
        buses = net2.buses
        @test buses.i == [1, 222222]
        @test buses.owner == [1, 7]

        # Test string column as expected
        @test buses.name[2] == "PRPR C D"

        loads = net2.loads
        @test loads.i == [1, 222222]  # first col
        @test loads.owner == [1, 7]   # last col
        @test isequal(loads.intrpt, [missing, missing])

        gens = net2.generators
        @test gens.i == [104]    # first col
        @test gens.f1 == [1.0]   # last col guaranteed to exist
        @test isequal(gens.wpf, [missing])
        @test gens.id[1] == "1" # string col

        branches = net2.branches
        @test branches.i == [1, 2, 222222]
        @test branches.j == [-543210, 9, 333333]  # negative numbers should be allowed
        @test branches.f1 == [1.0, 1.0, 1.0]
        @test branches.ckt[2] == "6"

        transformers = net2.transformers
        @test length(transformers.i) == 3
        @test transformers.i == [42, 4774, 222222]            #  1st entry of 1st row
        @test transformers.f1 == [1.0, 1.0, 1.0]              # last entry of 1st row
        @test transformers.r1_2 == [0.0025, 0.0, 0.0005]      #  1st entry of 2nd row
        @test transformers.sbase1_2 == [100.0, 100.0, 100.0]  # last entry of 2nd row
        @test transformers.windv1 == [1.0, 1.0462, 1.0]       #  1st entry of 3rd row
        @test transformers.cx1 == [0.0, 0.0, 0.0]             # last entry of 3rd row
        # Important to test a row where the 1st character is '0', to get it does not
        # get misinterpreted as the start of a "0 bus" records terminating the section.
        @test transformers.windv2 == [1.0, 1.045, 0.98250]    #  1st entry of 4th row
        @test transformers.nomv2 == [138.0, 240.35, 345.0]    # last entry of 4th row
        @test transformers.ckt == ["K1", "90", "B1"]          # string col

        # v29 testfile has only 2-winding data, so should return only 2-winding columns
        ncols_expected = 43
        @test length(Tables.columnnames(transformers)) == ncols_expected
        @test size(DataFrame(transformers)) == (3, ncols_expected)
        @test size(transformers) == (3, ncols_expected)
        @test length(transformers) == 3

        area_interchanges = net2.area_interchanges
        @test area_interchanges.i == [615, 762]
        @test area_interchanges.isw == [615001, 1234]
        @test area_interchanges.pdes == [32.677, -224.384]
        @test area_interchanges.ptol == [5.0, 5.0]
        @test area_interchanges.arname == ["RE", "OTP"]

        two_terminal_dc = net2.two_terminal_dc
        @test isempty(two_terminal_dc)

        vsc_dc = net2.vsc_dc
        @test isempty(vsc_dc)

        switched_shunts = net2.switched_shunts
        @test switched_shunts.i == [175, 177]      # first col
        @test switched_shunts.n1 == [2, 1]         # `n1` col present for both
        @test switched_shunts.b1 == [19.76, 15.60] # `b1` col present for both
        @test switched_shunts.n2 == [0, 1]         # `n2` col present only for second entry
        @test switched_shunts.b2 == [0.0, 17.69]   # `b2` col present only for second entry
        @test switched_shunts.n8 == [0, 0]         # `n8` col not present for either entry
        @test switched_shunts.b8 == [0.0, 0.0]     # last col; `b8` col not present for either entry

        impedance_corrections = net2.impedance_corrections
        @test impedance_corrections.i == [1, 2]          # first col
        @test impedance_corrections.f11 == [3.34, 1.129] # last col

        multi_terminal_dc = net2.multi_terminal_dc
        @test length(multi_terminal_dc) == 1

        mt_dc = only(multi_terminal_dc.lines)
        line_id = mt_dc.line_id
        @test line_id.i == 1     # first val
        @test line_id.vconvn == 0 # last val

        converters = mt_dc.converters
        @test length(converters) == line_id.nconv == 4
        @test converters.ib == [402, 401, 212, 213]
        @test converters.cnvcod == [3, 3, 1, 4]

        dc_buses = mt_dc.buses
        @test length(dc_buses) == line_id.ndcbs == 5
        @test dc_buses.idc == [1, 2, 3, 4, 5]
        @test dc_buses.owner == [4, 2, 4, 2, 4]

        dc_links = mt_dc.links
        @test length(dc_links) == line_id.ndcln == 4
        @test dc_links.idc == [1, 2, 3, 4]
        @test dc_links.ldc == [0.0, 0.0, 0.0, 0.0]

        multi_section_lines = net2.multi_section_lines
        @test isempty(multi_section_lines)

        zones = net2.zones
        @test zones.i == [1, 9]
        @test zones.zoname == ["ABL", "EFGN"]

        area_transfers = net2.area_transfers
        @test isempty(area_transfers)

        owners = net2.owners
        @test owners.i == [1, 2]
        @test owners.owname == ["EAI", "OUYK"]

        facts = net2.facts
        @test isempty(facts)
    end

    @testset "v33 file" begin
        net33 = parse_network("testfiles/synthetic_data_v33.RAW")
        @test net33 isa Network
        @test net33.version == 33

        caseid = net33.caseid
        @test caseid.ic == 0
        @test caseid.sbase == 100.0
        @test caseid.rev == 33

        buses = net33.buses
        @test buses.i == [1, 500]  # first col
        @test buses.name == ["WINNSBORO 0", "MC CORMICK 0"]  # string col
        @test buses.evlo == [0.9, 0.9]  # last col

        loads = net33.loads
        @test loads.i == [2, 500]  # first col
        @test loads.scale == [1, 1]  # last col
        @test isequal(loads.intrpt, [missing, missing])  # possible last col

        fixed_shunts = net33.fixed_shunts
        @test fixed_shunts.i == [320, 784]  # first col
        @test fixed_shunts.bl == [48.274, 19.939]  # last col

        gens = net33.generators
        @test gens.i == [9, 498]  # first col
        @test gens.wpf == [1.0, 1.0]  # last col

        branches = net33.branches
        @test branches.i == [2, 500]  # first col
        @test branches.met == [1, 1]  # only in v33
        @test branches.f4 == [1.0, 1.0]  # last col

        # v33 has extra columns on rows 1, 3, 4, 5.
        transformers = net33.transformers
        @test transformers.i == [8, 190, 498]                              #  1st entry of 1st row
        @test all(transformers.vecgrp .== "")                  # last entry of 1st row
        @test transformers.r1_2 == [3.55062E-4, 7.65222E-4, 4.00610E-3]    #  1st entry of 2nd row
        @test transformers.sbase1_2 == [100.0, 100.0, 100.0]               # last entry of 2nd row (T2)
        @test isequal(transformers.anstar, [missing, -65.843144, missing]) # last entry of 2nd row (T3)
        @test transformers.windv1 == [1.0, 1.0, 1.0]                       #  1st entry of 3rd row
        @test transformers.cnxa1 == [0.0, 0.0, 0.0]                        # last entry of 3rd row
        @test transformers.windv2 == [1.0, 1.0, 1.0]                       #  1st entry of 4th row
        @test transformers.nomv2 == [345.0, 115.0, 138.0]                  # last entry of 4th row (T2)
        @test isequal(transformers.cnxa2, [missing, 0.0, missing])         # last entry of 4th row (T3)
        @test isequal(transformers.windv3, [missing, 1.0, missing])        #  1st entry of 5th row
        @test isequal(transformers.cnxa3, [missing, 0.0, missing])         # last entry of 5th row
        @test transformers.ckt[1] == "1"                                  # string col

        # v33 testfile has both 2-winding and 3-winding data, so should return all columns
        @test length(Tables.columnnames(transformers)) == fieldcount(Transformers)
        @test size(DataFrame(transformers)) == (3, fieldcount(Transformers))
        @test size(transformers) == (3, fieldcount(Transformers))
        @test length(transformers) == 3

        area_interchanges = net33.area_interchanges
        @test area_interchanges.i == [1]
        @test area_interchanges.arname == ["SouthCarolin"]

        two_terminal_dc = net33.two_terminal_dc
        @test two_terminal_dc.name == ["DC Line 1", "DC Line 1"]  # first col
        @test two_terminal_dc.cccacc == [0.0, 0.0]         # last entry of 1st row
        @test two_terminal_dc.ipr == [2060653, 3008030]    #  1st entry of 2nd row
        @test two_terminal_dc.xcapr == [0.0, 0.0]          # last entry of 2nd row
        @test two_terminal_dc.ipi == [66353, 61477]        #  1st entry of 3nd row
        @test two_terminal_dc.xcapi == [0.0, 0.0]          # last entry of 3rd row

        vsc_dc = net33.vsc_dc
        @test isempty(vsc_dc)

        impedance_corrections = net33.impedance_corrections
        @test impedance_corrections.i == [1, 2, 3]          # first col
        @test impedance_corrections.f6 == [1.03, 0.0, 1.41] # col present on some rows not others
        @test impedance_corrections.f11 == [0.0, 0.0, 0.0]  # last col (missing, default to zero)

        multi_terminal_dc = net33.multi_terminal_dc
        @test length(multi_terminal_dc) == 1

        mt_dc = only(multi_terminal_dc.lines)
        line_id = mt_dc.line_id
        @test line_id.name == "DC Line 1" # first val
        @test line_id.vconvn == 0            # last val

        multi_section_lines = net33.multi_section_lines
        @test multi_section_lines.i == [1]                 # first col
        @test multi_section_lines.id == ["&1"]             # string col
        @test multi_section_lines.met == [1]               # only in v33 data
        @test isequal(multi_section_lines.dum1, [3])       # `dum1` is last col present
        @test isequal(multi_section_lines.dum9, [missing]) # last col; not present; default to missing

        zones = net33.zones
        @test zones.i == [1, 2]
        @test zones.zoname == ["Upstate", "Midlands"]

        area_transfers = net33.area_transfers
        @test isempty(area_transfers)

        owners = net33.owners
        @test owners.i == [1]
        @test owners.owname == ["1"]

        facts = net33.facts
        @test isempty(facts)

        switched_shunts = net33.switched_shunts
        @test switched_shunts.i == [27, 491]      # first col
        @test switched_shunts.adjm == [0, 0]      #   3rd col (only in v33)
        @test switched_shunts.stat == [1, 1]      #   4th col (only in v33)
        @test switched_shunts.n1 == [1, 1]        #  `n1` col always present
        @test switched_shunts.b1 == [80.0, 50.0]  #  `b1` col always present
        @test switched_shunts.n8 == [0, 0]         # `n8` col not present; default to zero
        @test switched_shunts.b8 == [0.0, 0.0]     # `b8` col not present; default to zero
    end

    @testset "issues" begin
        net = parse_network("testfiles/spacezero.raw")
        @test length(net.buses) == 2
        @test length(net.loads) == 1
    end

    @testset "empty network" begin

        empty_net = Network(
            30,
            CaseID(),
            Buses30(),
            Loads(),
            nothing,
            Generators(),
            Branches30(),
            Transformers(),
            AreaInterchanges(),
            TwoTerminalDCLines30(),
            VSCDCLines(),
            SwitchedShunts30(),
            ImpedanceCorrections(),
            MultiTerminalDCLines(),
            MultiSectionLineGroups30(),
            Zones(),
            InterAreaTransfers(),
            Owners(),
            FACTSDevices30(),
        )

        full_net = parse_network("testfiles/synthetic_data_v30.raw")

        @test isempty(empty_net)
        @test !isempty(full_net)
    end
end


text:
"""

```
eulermat(beta, gamma, theta, deg::String)
```

This routine was added because it defines
the rotation in the "human" way, an is thus used
to set the position of the fixed molecules. `deg` can only be `"degree"`, in which
case the angles with be considered in degrees. If no `deg` argument
is provided, radians are used.

That means: `beta` is a counterclockwise rotation around `x` axis.
            `gamma` is a counterclockwise rotation around `y` axis.
            `theta` is a counterclockwise rotation around `z` axis.

"""
function eulermat(beta, gamma, theta, deg::String)

    if deg != "degree"
        error("ERROR: to use radians just omit the last parameter")
    end

    beta = beta * π / 180
    gamma = gamma * π / 180
    theta = theta * π / 180

    return eulermat(beta, gamma, theta)
end

function eulermat(beta::T, gamma::T, theta::T) where {T}
    c1 = cos(beta)
    s1 = sin(beta)
    c2 = cos(gamma)
    s2 = sin(gamma)
    c3 = cos(theta)
    s3 = sin(theta)
    @SMatrix [
        c2*c3 -c2*s3 s2
        (c1*s3+c3*s1*s2) (c1*c3-s1*s2*s3) -c2*s1
        (s1*s3-c1*c3*s2) (c1*s2*s3+c3*s1) c1*c2
    ]
end

ext: jl


text:
## Test pvalue adjustment methods ##
module Test_pval_adjustment

using MultipleTesting
using Base.Test

pi0 = 0.4
pval1 = [0.0, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.2, 0.4, 0.75, 1.0]
ref1 = Dict([
             (Bonferroni, [0.0, 0.001, 0.01, 0.1, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0]),
             (Holm, [0.0, 9e-4, 8e-3, 0.07, 0.3, 0.5, 0.8, 1.0, 1.0, 1.0]),
             (Hochberg, [0.0, 9e-4, 8e-3, 0.07, 0.3, 0.5, 0.8, 1.0, 1.0, 1.0]),
             (Hommel, [0.0, 9e-4, 8e-3, 0.07, 0.3, 0.5, 0.8, 1.0, 1.0, 1.0]),
             (BenjaminiHochberg, [0.0, 0.0005, 0.003333333, 0.025, 0.1, 0.166666667, 0.285714286, 0.5, 0.833333333, 1.0]),
             (BenjaminiYekutieli, [0.0, 0.001464484, 0.009763228, 0.073224206, 0.292896825, 0.488161376, 0.836848073, 1.0, 1.0, 1.0]),
             (BenjaminiLiu, [0.0, 0.0008096761, 0.0063776447, 0.0475542565, 0.1589448656, 0.2047550000, 0.2361600000, 0.2361600000, 0.2361600000, 0.2361600000]), #mutoss
             (Sidak, [ 0.0, 0.0009995501, 0.0099551198, 0.0956179250, 0.4012630608, 0.6513215599, 0.8926258176, 0.9939533824, 0.9999990463, 1.0000000000]), #mutoss
             (ForwardStop, [0.0, 0.0000500025, 0.0003668351, 0.0027877103, 0.0124888271, 0.0279674419, 0.0558497432, 0.1127217283, 0.2542297986, 1.0])
             ])

pval2 = [0.0001, 0.0001, 0.001, 0.01, 0.05, 0.1, 0.1, 0.4, 0.75, 1.0]
ref2 = Dict([
             (Bonferroni, [0.001, 0.001, 0.01, 0.1, 0.5, 1.0, 1.0, 1.0, 1.0, 1.0]),
             (Holm, [0.001, 0.001, 0.008, 0.07, 0.3, 0.5, 0.5, 1.0, 1.0, 1.0]),
             (Hochberg, [9e-4, 9e-4, 8e-3, 0.07, 0.3, 0.4, 0.4, 1.0, 1.0, 1.0]),
             (Hommel, [9e-4, 9e-4, 8e-3, 0.07, 0.25, 0.4, 0.4, 1.0, 1.0, 1.0]),
             (BenjaminiHochberg, [0.0005, 0.0005, 0.003333333, 0.025, 0.1, 0.142857143, 0.142857143, 0.5, 0.833333333, 1.0]),
             (BenjaminiYekutieli, [0.001464484, 0.001464484, 0.009763228, 0.073224206, 0.292896825, 0.418424036, 0.418424036, 1.0, 1.0, 1.0]),
             (BenjaminiLiu, [0.0009995501, 0.0009995501, 0.0063776447, 0.0475542565, 0.1589448656, 0.2047550000, 0.2047550000, 0.2352000000, 0.2352000000, 0.2352000000]), #mutoss
             (Sidak,  [0.0009995501, 0.0009995501, 0.0099551198, 0.0956179250, 0.4012630608, 0.6513215599, 0.6513215599, 0.9939533824, 0.9999990463, 1.0000000000]), #mutoss
             (ForwardStop, [0.0001000050, 0.0001000050, 0.0004001701, 0.0028127115, 0.0125088281, 0.0279841094, 0.0390378817, 0.0980113495, 0.2411539063, 1.0])
             ])

for method in keys(ref1)
    println(" ** ", method)
    @test issubtype(method, PValueAdjustmentMethod)
    @test issubtype(typeof(method()), PValueAdjustmentMethod)
    @test_throws MethodError method(0.1)
    ## no valid p-values as input
    @test_throws DomainError adjust([-1.0, 0.5], method())
    @test_throws DomainError adjust([0.5, 1.5], method())
    ## single p-value is returned unchanged
    if method != ForwardStop  # this test is not valid for ForwardStop
        pval = rand(1)
        @test adjust(pval, method()) == pval
    end
    ## compare with reference values
    @test_approx_eq_eps adjust(pval1, method()) ref1[method] 1e-9
    ## compare with reference values having ties
    @test_approx_eq_eps adjust(pval2, method()) ref2[method] 1e-9
end

end


text:
# ===========================================================================
# All constants needed by tests are here
const unicode_chars = String.(readdlm(path*"/SampleFiles/julia-unicode.csv", '\n')[:,1])
const n_unicode = length(unicode_chars)
const breaking_dict = Dict{String, Any}(
    "0" => rand(Char),
    "1" => randstring(2^rand(2:10)),
   "16" => rand(UInt8),
   "17" => rand(UInt16),
   "18" => rand(UInt32),
   "19" => rand(UInt64),
   "20" => rand(UInt128),
   "32" => rand(Int8),
   "33" => rand(Int16),
   "34" => rand(Int32),
   "35" => rand(Int64),
   "36" => rand(Int128),
   "48" => rand(Float16),
   "49" => rand(Float32),
   "50" => rand(Float64),
   "80" => rand(Complex{UInt8}),
   "81" => rand(Complex{UInt16}),
   "82" => rand(Complex{UInt32}),
   "83" => rand(Complex{UInt64}),
   "84" => rand(Complex{UInt128}),
   "96" => rand(Complex{Int8}),
   "97" => rand(Complex{Int16}),
   "98" => rand(Complex{Int32}),
   "99" => rand(Complex{Int64}),
  "100" => rand(Complex{Int128}),
  "112" => rand(Complex{Float16}),
  "113" => rand(Complex{Float32}),
  "114" => rand(Complex{Float64}),
  "128" => rand(Char, 2^rand(2:6)),
  "129" => [randstring(2^rand(3:8)) for i = 1:rand(4:24)],
  "144" => rand(UInt8, rand(4:24)),
  "145" => rand(UInt16, rand(4:24)),
  "146" => rand(UInt32, rand(4:24)),
  "147" => rand(UInt64, rand(4:24)),
  "148" => rand(UInt128, rand(4:24)),
  "160" => rand(Int8, rand(4:24)),
  "161" => rand(Int16, rand(4:24)),
  "162" => rand(Int32, rand(4:24)),
  "163" => rand(Int64, rand(4:24)),
  "164" => rand(Int128, rand(4:24)),
  "176" => rand(Float16, rand(4:24)),
  "177" => rand(Float32, rand(4:24)),
  "178" => rand(Float64, rand(4:24)),
  "208" => rand(Complex{UInt8}, rand(4:24)),
  "209" => rand(Complex{UInt16}, rand(4:24)),
  "210" => rand(Complex{UInt32}, rand(4:24)),
  "211" => rand(Complex{UInt64}, rand(4:24)),
  "212" => rand(Complex{UInt128}, rand(4:24)),
  "224" => rand(Complex{Int8}, rand(4:24)),
  "225" => rand(Complex{Int16}, rand(4:24)),
  "226" => rand(Complex{Int32}, rand(4:24)),
  "227" => rand(Complex{Int64}, rand(4:24)),
  "228" => rand(Complex{Int128}, rand(4:24)),
  "240" => rand(Complex{Float16}, rand(4:24)),
  "242" => rand(Complex{Float64}, rand(4:24)),
  "241" => rand(Complex{Float32}, rand(4:24))
  )
const NOOF = "
  HI ITS CLOVER LOL
          ,-'-,  `---..
         /             \\
         =,             .
  ______<3.  ` ,+,     ,\\`
 ( \\  + `-”.` .; `     `.\\
 (_/   \\    | ((         ) \\
  |_ ;  \"    \\   (        ,’ |\\
  \\    ,- '💦 (,\\_____,’   / “\\
   \\__---+ }._)              |\\
   / _\\__💧”)/                  +
  ( /    💧” \\                  ++_
   \\)    ,“  |)                ++  ++
   💧     “💧  (                 *    +***
"


text:
# TODO:
# - midpoint rule
# x scattered/interpolation
# x GCV
# x L-curve
# - bounds on λ
# - initial guess for λ? will require mods to RegularizationTools
# - scaled λ? will need to work out equivalency with λ² formulation, and resolve with
#   derivative rather than difference matrix
# - standard deviation
# - relative weights? (not sure if worth implementing)
# - arbitrary weighting

const LA = LinearAlgebra

# might move this to RegularizationTools.jl? JJS 12/6/21
"""
Order d derivative matrix for the provided t vector
"""
function _derivative_matrix(t::AbstractVector, d::Int)
    N = length(t)
    if d == 0
        return Array{Float64}(LA.I, (N, N))
    end
    dt = t[d+1:end] - t[1:end-d]
    V = LA.diagm(1 ./ dt)
    # could use this symbol:  Dᵈ⁻¹ (D\^d\^-\^1) (but a lot of typing...)
    Ddm1 = diff(_derivative_matrix(t, d-1), dims=1)
    D = d * V*Ddm1
    return D
end

### Regularization data smoothing and interpolation
struct RegularizationSmooth{uType,tType,FT,T} <: AbstractInterpolation{FT,T}
    u::uType
    û::uType
    t::tType
    t̂::tType
    d::Int              # derivative degree used to calculate the roughness
    λ::AbstractFloat    # regularization parameter
    alg::Symbol  # how to determine λ: `:fixed`, `:gcv_svd`, `:gcv_tr`, `L_curve`
    Aitp::AbstractInterpolation{FT,T}
    RegularizationSmooth{FT}(u,û,t,t̂,d,λ,alg,Aitp) where FT =
        new{typeof(u),typeof(t),FT,eltype(u)}(u,û,t,t̂,d,λ,alg,Aitp)
end

function RegularizationSmooth(u, t, d=2; λ=1.0, t̂=nothing, alg=:gcv_svd)
    # TBD:  make use of `munge_data`
    u, t = munge_data(u, t)
    λ = float(λ) # in case an integer is given for λ
    N = length(t)
    if isnothing(t̂)
        t̂ = t
        M = Array{eltype(t)}(LA.I, (N, N)) # identity mapping matrix
    else
        M = _mapping_matrix(t̂, t)
    end

    # RT.Γ function, used internally there, does not not account for spacing of x
    #Ψ = RT.setupRegularizationProblem(M, d)
    D = _derivative_matrix(t̂,d)
    Ψ = RT.setupRegularizationProblem(M, D)

    if alg == :fixed
        b̄ = RT.to_standard_form(Ψ, u) # via b̄
        ū = RT.solve(Ψ, b̄, λ)
        û = RT.to_general_form(Ψ, u, ū)
    else
        # note that the provided λ is currently ignored
        result = RT.solve(Ψ, u, alg=alg)
        û = result.x
        λ = result.λ
    end
    # I'm not sure about mathematical theory for performing a posteriori interpolation via
    # Regularization smoothing -- applying Cubic-spline interpolation seems reasonable to
    # me, so setting that up here to use after object-creation; could allow users the option
    # of a posteriori interpolation or not, or even a choice for the interpolation,
    # depending on what is needed. B-spline interpolation of order d+1 also seems to be a
    # reasonable choice, but it is not allowing extrapolation and doesn't show much benefit
    # over cubic-spline. JJS 12/3/21
    Aitp = CubicSpline(û,t̂)
    #Aitp = BSplineInterpolation(û,t̂,d+1,:ArcLen,:Average)
    RegularizationSmooth{true}(u,û,t,t̂,d,λ,alg,Aitp)
end

"""Linear interpolation mapping matrix, which maps `û` to `u`."""
function _mapping_matrix(t̂::AbstractVector,t::AbstractVector)
    N = length(t)
    N̂ = length(t̂)
    # map the scattered points to the appropriate index of the smoothed points
    idx = searchsortedlast.(Ref(t̂),t)
    # allow for "extrapolation"; i.e. for t̂ extremum that are interior to t
    idx[idx.==0] .+= 1
    idx[idx.==N̂] .+= -1
    # create the linear interpolation matrix
    m2 = @. (t - t̂[idx])/(t̂[idx+1] - t̂[idx])
    M = zeros(eltype(t), (N,N̂))
    for i in 1:N
        M[i,idx[i]] = 1 - m2[i]
        M[i,idx[i]+1] = m2[i]
    end
    return M
end

_interpolate(A::RegularizationSmooth{<:AbstractVector{<:Number}}, t::Number) =
    _interpolate(A.Aitp, t)


text:
module NumericalRelativity

# Tensor
#
# (array of functions, coordinates, (indices up or down)) -> tensor
# T[-1,-1] means lower-lower indices and T[1,1] means up-up
# You should be able to find the Rank of tensors
# Contracting is possible? At least with a function contract

using Symbolics

function partial(tensor, coords)
	[(expand_derivatives∘Differential(coord))(tencom) for coord ∈ coords, tencom ∈ tensor]
end

# Metric should be a rank two tensor
#
# partial gives back the symbols ∂T
# christoffel symbols should be (metric) -> christoffel
#
# Riemann Tensor g -> Riemann Tensor
# Ricci Tensor g -> Ricci Tensor
# Ricci Scalar g -> Ricci Scalar

end

max_stars_repo_name: kawazugit/NumericalRelativity

max_issues_repo_name: kawazugit/NumericalRelativity

max_forks_repo_name: kawazugit/NumericalRelativity


text:
@testset "MapFun" begin
    @test (MapFun <: Augmentor.ImageOperation) == false
    @test (MapFun <: Augmentor.Operation) == true
    @testset "constructor" begin
        @test_throws MethodError MapFun()
        @test typeof(@inferred(MapFun((x)->x))) <: MapFun <: Augmentor.Operation
        @test typeof(@inferred(MapFun(identity))) <: MapFun <: Augmentor.Operation
        @test str_show(MapFun(identity)) == "Augmentor.MapFun(identity)"
        @test str_showconst(MapFun(identity)) == "MapFun(identity)"
        @test str_showcompact(MapFun(identity)) == "Map function \"identity\" over image"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(MapFun(identity), nothing)
        @test Augmentor.supports_eager(MapFun) === true
        img_out = map(x->x-Gray{Float32}(0.1), rect)
        imgs = [
            rect,
            view(rect,:,:),
            Augmentor.prepareaffine(rect),
            OffsetArray(rect, 0, 0),
            view(rect, IdentityRange(1:2), IdentityRange(1:3))
        ]
        @testset "single image" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(MapFun(identity), img_in))
                @test res == img_in
                @test eltype(res) <: Gray{N0f8}
                @test typeof(axes(img_in)) <: NTuple{2,Base.OneTo} ? typeof(res) <: Array : typeof(res) <: Array
                res = @inferred(Augmentor.applyeager(MapFun(x->x-Gray{Float32}(0.1)), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
            img = OffsetArray(rgb_rect, -2, -1)
            res = @inferred(Augmentor.applyeager(MapFun(x -> x - RGB(.1,.1,.1)), img))
            @test res ≈ collect(img .- RGB(.1,.1,.1))
            @test typeof(res) <: Array{RGB{Float64}}
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(MapFun) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(MapFun) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(MapFun) === true
        res = @inferred(Augmentor.applylazy(MapFun(identity), rect))
        @test res === mappedarray(identity, rect)
        res = @inferred(Augmentor.applylazy(MapFun(identity), rgb_rect))
        @test res === mappedarray(identity, rgb_rect)
        res = @inferred(Augmentor.applylazy(MapFun(x->x-RGB(.1,.1,.1)), rgb_rect))
        @test res == mappedarray(x->x-RGB(.1,.1,.1), rgb_rect)
        res = @inferred(Augmentor.applylazy(MapFun(x->x-RGB(.1,.1,.1)), OffsetArray(rgb_rect,-2,-1)))
        @test axes(res) == OffsetRange.((-1:0, 0:2))
        @test @inferred(getindex(res,0,0)) isa RGB{Float64}
        @test res == mappedarray(x->x-RGB(.1,.1,.1), OffsetArray(rgb_rect,-2,-1))
        @test typeof(res) <: MappedArrays.ReadonlyMappedArray{ColorTypes.RGB{Float64}}
    end
    @testset "view" begin
        @test Augmentor.supports_view(MapFun) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(MapFun) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(MapFun) === false
    end
end

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

@testset "AggregateThenMapFun" begin
    @test (AggregateThenMapFun <: Augmentor.ImageOperation) == false
    @test (AggregateThenMapFun <: Augmentor.Operation) == true
    @testset "constructor" begin
        @test_throws MethodError AggregateThenMapFun()
        @test_throws MethodError AggregateThenMapFun(x->x)
        @test typeof(@inferred(AggregateThenMapFun(x->x, x->x))) <: AggregateThenMapFun <: Augmentor.Operation
        @test typeof(@inferred(AggregateThenMapFun(mean, identity))) <: AggregateThenMapFun <: Augmentor.Operation
        mean_str = VERSION >= v"1.2.0" ? "mean" : "Statistics.mean"
        @test str_show(AggregateThenMapFun(mean, identity)) == "Augmentor.AggregateThenMapFun($mean_str, identity)"
        @test str_showconst(AggregateThenMapFun(mean, identity)) == "AggregateThenMapFun($mean_str, identity)"
        @test str_showcompact(AggregateThenMapFun(mean, identity)) == "Map result of \"mean\" using \"identity\" over image"
    end
    @testset "eager" begin
        @test_throws MethodError Augmentor.applyeager(AggregateThenMapFun(mean, identity), nothing)
        @test Augmentor.supports_eager(AggregateThenMapFun) === true
        m = mean(rect)
        img_out = map(x->x-m, rect)
        imgs = [
            rect,
            view(rect,:,:),
            Augmentor.prepareaffine(rect),
            OffsetArray(rect, 0, 0),
            view(rect, IdentityRange(1:2), IdentityRange(1:3))
        ]
        @testset "single image" begin
            for img_in in imgs
                res = @inferred(Augmentor.applyeager(AggregateThenMapFun(mean, (x,a)->x), img_in))
                @test res == img_in
                @test eltype(res) <: Gray{N0f8}
                @test typeof(axes(img_in)) <: NTuple{2,Base.OneTo} ? typeof(res) <: Array : typeof(res) <: Array
                res = @inferred(Augmentor.applyeager(AggregateThenMapFun(mean, (x,a)->x-a), img_in))
                @test res == img_out
                @test typeof(res) == typeof(img_out)
            end
            img = OffsetArray(rgb_rect, -2, -1)
            res = @inferred(Augmentor.applyeager(AggregateThenMapFun(mean, (x,a)->x-a), img))
            @test res ≈ collect(img .- mean(rgb_rect))
            @test typeof(res) <: Array{RGB{Float64}}
        end
    end
    @testset "affine" begin
        @test Augmentor.supports_affine(AggregateThenMapFun) === false
    end
    @testset "affineview" begin
        @test Augmentor.supports_affineview(AggregateThenMapFun) === false
    end
    @testset "lazy" begin
        @test Augmentor.supports_lazy(AggregateThenMapFun) === true
        res = @inferred(Augmentor.applylazy(AggregateThenMapFun(mean, (x,a)->x), rect))
        @test res == rect
        @test parent(res) === rect
        @test res isa ReadonlyMappedArray
        res = @inferred(Augmentor.applylazy(AggregateThenMapFun(mean, (x,a)->x-a), rgb_rect))
        @test res == mappedarray(x->x-mean(rgb_rect), rgb_rect)
        @test parent(res) === rgb_rect
        @test typeof(res) <: MappedArrays.ReadonlyMappedArray{ColorTypes.RGB{Float64}}
        img = OffsetArray(rgb_rect, -2, -1)
        res = @inferred(Augmentor.applylazy(AggregateThenMapFun(mean, (x,a)->x-a), img))
        @test res == mappedarray(x->x-mean(rgb_rect), img)
        @test parent(res) === img
        @test typeof(res) <: MappedArrays.ReadonlyMappedArray{ColorTypes.RGB{Float64}}
        @test @inferred(getindex(res,0,0)) isa RGB{Float64}
    end
    @testset "view" begin
        @test Augmentor.supports_view(AggregateThenMapFun) === false
    end
    @testset "stepview" begin
        @test Augmentor.supports_stepview(AggregateThenMapFun) === false
    end
    @testset "permute" begin
        @test Augmentor.supports_permute(AggregateThenMapFun) === false
    end
end


text:
using CUSPARSE
using CUDArt
using Base.Test

m = 25
n = 35
k = 10
blockdim = 5

##############
# test_csrmm #
##############

for (func,funcs,funch,func2,typ) in ((:test_csrmm,:test_csrmm_symm,:test_csrmm_herm,:test_csrmm2,:CudaSparseMatrixCSR),
                                     (:test_cscmm,:test_cscmm_symm,:test_cscmm_herm,:test_cscmm2,:CudaSparseMatrixCSC))
    @eval begin
        function $funcs(elty)
            A = sparse(rand(elty,m,m))
            A = A + A.'
            B = rand(elty,m,n)
            C = rand(elty,m,n)
            alpha = rand(elty)
            beta = rand(elty)
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = Symmetric($typ(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = to_host(d_C)
            C = alpha * A * B + beta * C
            @test_approx_eq(C,h_C)
            d_C = d_A.' * d_B
            h_C = to_host(d_C)
            C = A.' * B
            @test_approx_eq(C,h_C)
            d_B = CudaArray(rand(elty,k,n))
            @test_throws(DimensionMismatch, CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O'))
        end

        function $funch(elty)
            A = sparse(rand(elty,m,m))
            A = A + A'
            B = rand(elty,m,n)
            C = rand(elty,m,n)
            alpha = rand(elty)
            beta = rand(elty)
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = Hermitian($typ(A))
            d_C = CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_C,'O')
            h_C = to_host(d_C)
            C = alpha * A * B + beta * C
            @test_approx_eq(C,h_C)
            d_B = CudaArray(rand(elty,k,n))
            @test_throws(DimensionMismatch, CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O'))
        end

        function $func(elty)
            A = sparse(rand(elty,m,k))
            B = rand(elty,k,n)
            C = rand(elty,m,n)
            alpha = rand(elty)
            beta = rand(elty)
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = $typ(A)
            @test_throws(DimensionMismatch, CUSPARSE.mm!('T',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm!('N',alpha,d_A,d_B,beta,d_B,'O'))
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = alpha * A * B + beta * C
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm('N',d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + beta * C
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm('N',d_A,d_B,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + C
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm('N',alpha,d_A,d_B,'O')
            h_D = to_host(d_D)
            D = alpha * A * B
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm('N',d_A,d_B,'O')
            h_D = to_host(d_D)
            D = A * B
            @test_approx_eq(D,h_D)
            d_D = d_A*d_B
            h_D = to_host(d_D)
            D = A * B
            @test_approx_eq(D,h_D)
            @test_throws(DimensionMismatch, CUSPARSE.mm('T',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm('N',alpha,d_A,d_B,beta,d_B,'O'))
        end
        function $func2(elty)
            A = sparse(rand(elty,m,k))
            B = rand(elty,k,n)
            C = rand(elty,m,n)
            alpha = rand(elty)
            beta = rand(elty)
            d_B = CudaArray(B)
            d_C = CudaArray(C)
            d_A = $typ(A)
            @test_throws(DimensionMismatch, CUSPARSE.mm2!('N','T',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm2!('T','N',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm2!('T','T',alpha,d_A,d_B,beta,d_C,'O'))
            @test_throws(DimensionMismatch, CUSPARSE.mm2!('N','N',alpha,d_A,d_B,beta,d_B,'O'))
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = alpha * A * B + beta * C
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + beta * C
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
            h_D = to_host(d_D)
            D = A * B + C
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
            h_D = to_host(d_D)
            D = alpha * A * B
            @test_approx_eq(D,h_D)
            d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
            h_D = to_host(d_D)
            D = A * B
            @test_approx_eq(D,h_D)
        end
    end
end

function test_bsrmm2(elty)
    A = sparse(rand(elty,m,k))
    B = rand(elty,k,n)
    C = rand(elty,m,n)
    alpha = rand(elty)
    beta = rand(elty)
    d_B = CudaArray(B)
    d_C = CudaArray(C)
    d_A = CudaSparseMatrixCSR(A)
    d_A = CUSPARSE.switch2bsr(d_A,convert(Cint,blockdim))
    @test_throws(DimensionMismatch, CUSPARSE.mm2('N','T',alpha,d_A,d_B,beta,d_C,'O'))
    @test_throws(DimensionMismatch, CUSPARSE.mm2('T','N',alpha,d_A,d_B,beta,d_C,'O'))
    @test_throws(DimensionMismatch, CUSPARSE.mm2('T','T',alpha,d_A,d_B,beta,d_C,'O'))
    @test_throws(DimensionMismatch, CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_B,'O'))
    d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,beta,d_C,'O')
    h_D = to_host(d_D)
    D = alpha * A * B + beta * C
    @test_approx_eq(D,h_D)
    d_D = CUSPARSE.mm2('N','N',d_A,d_B,beta,d_C,'O')
    h_D = to_host(d_D)
    D = A * B + beta * C
    @test_approx_eq(D,h_D)
    d_D = CUSPARSE.mm2('N','N',d_A,d_B,d_C,'O')
    h_D = to_host(d_D)
    D = A * B + C
    @test_approx_eq(D,h_D)
    d_D = CUSPARSE.mm2('N','N',alpha,d_A,d_B,'O')
    h_D = to_host(d_D)
    D = alpha * A * B
    @test_approx_eq(D,h_D)
    d_D = CUSPARSE.mm2('N','N',d_A,d_B,'O')
    h_D = to_host(d_D)
    D = A * B
    @test_approx_eq(D,h_D)
end

types = [Float32,Float64,Complex64,Complex128]
for elty in types
    tic()
    test_csrmm_symm(elty)
    test_cscmm_symm(elty)
    println("mm_symm! took ", toq(), " for ", elty)
    tic()
    test_csrmm_herm(elty)
    test_cscmm_herm(elty)
    println("mm_herm! took ", toq(), " for ", elty)
    tic()
    test_csrmm(elty)
    test_cscmm(elty)
    println("mm took ", toq(), " for ", elty)
    tic()
    test_csrmm2(elty)
    test_cscmm2(elty)
    test_bsrmm2(elty)
    println("mm2 took ", toq(), " for ", elty)
end


text:
"""
    qrdelete!(Q, R)

Delete the left-most column of F = Q*R by updating Q and R.
The right-most column of Q and last column/last row of R is
invalid on exit.
"""
# TODO should only run up to m_eff!
function qrdelete!(Q::AbstractMatrix{T}, R) where T
    n, m = size(Q)
    for i = 1:m-1
        #temp = sqrt(R[i, i+1]^2 + R[i+1, i+1]^2)
        temp = hypot(R[i, i+1], R[i+1, i+1])
        c = R[i, i+1]/temp
        s = R[i+1, i+1]/temp
        R[i, i+1] = temp
        R[i+1, i+1] = T(0)
        if i < m-1
            for j = i+2:m
                temp = c*R[i, j] + s*R[i+1, j]
                R[i+1, j] = -s*R[i, j] + c*R[i+1, j]
                R[i, j] = temp
            end
        end

        for l = 1:n
            temp = c*Q[l, i] + s*Q[l, i+1]
            Q[l, i+1] = -s*Q[l, i] + c*Q[l, i+1]
            Q[l, i] = temp
        end
    end
    # my addition: move it up into the corner
    for j = 2:m
        for i = 1:m-1
            R[i, j-1] = R[i, j]
        end
    end
end


text:
using KDTrees

function run_bench_query_ball()
    dims = 3
    n_points = [10^i for i in 3:5]
    rs = [0.1 0.2 0.3 0.4]

    times = fill(0.0, length(rs), length(n_points))

    n_iters = 2000

    # Compile it
    tree = KDTree(randn(2,2))
    inball(tree, zeros(2), 0.1)

    for (i, r) in enumerate(rs)
        for (j , n_point) in enumerate(n_points)
            data = rand(dims, n_point)
            tree = KDTree(data, 10, true)
            t = time_ns()
            for z = 1:n_iters
                p = rand(dims)
                idxs = inball(tree, p, r)
            end
            t = (time_ns() - float(t)) / 10^9
            times[i,j] =  t / n_iters
        end
    end

    println(times)
end

run_bench_query_ball()

# Removed old benchmarks because apparently the sort in the
# end was taking most of the time making the results useless.


text:
using DataGenerators

function try_generating_for_type(t::Type, supplementalts::Vector{DataType} = Vector{DataType}(); tries::Int = 100, showdatum::Bool=false, showerror::Bool=false)

	println("type $(t):")

	g = type_generator(:TypeGen, t, supplementalts)

	gi = g()

	print("  generating $(tries) data")
	if showdatum
		println()
	end

	for i in 1:100

		x = try
			choose(gi)
		catch exc
			if isa(exc, DataGeneratorTranslators.TypeGenerationException)
				if showerror
					warn("$(exc)")
				else
					print("!")
				end
				continue
			else
				rethrow(exc)
			end
		end

		if showdatum
			println("Type $(typeof(x)); Value: $(x)")
		else
			print(".")
		end

		@test DataGeneratorTranslators.isa_tuple_adjust(x, t)

	end
	println()

	println()

end

@testset "smoke tests" begin

	sd = false
	se = false

	try_generating_for_type(Int16, showdatum=sd, showerror=se)

	try_generating_for_type(Union{Int8, UInt16}, showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{Int32, String}, showdatum=sd, showerror=se)

	try_generating_for_type(Signed, showdatum=sd, showerror=se)

end

@testset "regression tests" begin

	sd = false
	se = false

	try_generating_for_type(Signed, showdatum=sd, showerror=se)

	try_generating_for_type(Any, DataType[Number], showdatum=sd, showerror=se)

	try_generating_for_type(Union{Int8, UInt16}, showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{Int32, String}, showdatum=sd, showerror=se)

	try_generating_for_type(Type{Unsigned}, DataType[Number], showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{TypeVar(:S, Signed, false), TypeVar(:S, Signed, false)}, showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{TypeVar(:S, Signed, true), TypeVar(:T, Signed, true)}, showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{TypeVar(:S, Signed, true), TypeVar(:S, Signed, true)}, showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{TypeVar(:S, Signed, true), Type{TypeVar(:S, Signed, true)}}, showdatum=sd, showerror=se)

	try_generating_for_type(Tuple{TypeVar(:S, Signed, false), Type{TypeVar(:T, Signed, false)}}, showdatum=sd, showerror=se)

end

@testset "dogfood tests" begin

	typeg = type_generator(:TypeGen, Type, DataType[Number;])

	typegi = typeg()

	for i in 1:10

		t = nothing
		while t == nothing
			try
				t = choose(typegi)
				if t == Union{}
					warn("skipping type returned from type generator as it is Union{}")
					t = nothing
					continue
				end
			catch exc
				if isa(exc, DataGeneratorTranslators.TypeGenerationException)
					warn("skipping type returned from type generator owing to: $(exc)")
					continue
				else
					rethrow(exc)
				end
			end
		end

		try_generating_for_type(t, convert(Vector{DataType},DataGenerators.GENERATOR_SUPPORTED_CHOOSE_TYPES), showdatum=false, showerror=false)

	end

end

@testset "edge tests" begin

	sd = false
	se = false

	try_generating_for_type(Tuple, DataType[Int64, String], showdatum=sd, showerror=se)

	try_generating_for_type(DataType, DataType[Number], showdatum=sd, showerror=se)

	try_generating_for_type(Union, DataType[Number], showdatum=sd, showerror=se)

	try_generating_for_type(Type{TypeVar(:T, Integer)}, showdatum=sd, showerror=se)

	try_generating_for_type(Type{TypeVar(:T, Signed)}, DataType[Number], showdatum=sd, showerror=se)

	try_generating_for_type(Type{Any}, DataType[Number], showdatum=sd, showerror=se)

	try_generating_for_type(Type, DataType[Number], showdatum=sd, showerror=se)

end


text:
using Test
using Distributions
include("transforms.jl")
include("legpts.jl")
include("statistics.jl")

@testset "transforms.jl" begin
    @test invBoxCox(boxCox(2, 3), 3) ≈ 2.0 atol = 1e-15
    @test invBoxCox(boxCox(2, 1), 1) ≈ 2.0 atol = 1e-15
    @test boxCox(invBoxCox(2, 3), 3) ≈ 2.0 atol = 1e-15
    @test boxCox(invBoxCox(2, 1), 1) ≈ 2.0 atol = 1e-15
    @test boxCox(invBoxCox([1 2 3], 1), 1) ≈ [1 2 3] atol = 1e-15
    @test boxCox(invBoxCox([1; 2], 1), 1) ≈ [1;2] atol = 1e-15
    @test boxCox(invBoxCox(2., -0.45), -0.45) ≈ 2.0 atol = 1e-15
    @test invBoxCox(boxCox(2., -2), -2) ≈ 2.0 atol = 1e-15
end

@testset "legpts.jl" begin
    f = x -> exp(x)+x^2
    @test_broken int1D(f, -1, 1) ≈ 3.01706 atol = 1e-5
end

@testset "statistics" begin
    log = LogNormal(2, 1)
    y = x-> Distributions.cdf(log, x)-0.5
    med = median(log)
    @test_broken bisection(y, 0, 15) ≈  med atol = 1e-3
end

@testset ""

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
gnorms = mean(luad_norm, 2)
Tdesc = ["5th quantile", "Median", "95th quantile", "IQR", "Variance"]
order = sortperm(gnorms[genes])
genes = genes[order]
pvals = pvals[order,:]
true_ts = true_ts[order,:]
tvals = tvals[order,:,:]

ints = mapslices(x->quantile(x, [0.05,0.95]), tvals, 3)

for i=1:length(Ts)
   println("")
   println("""\\begin{table*}[h!]
\\label{tab:pvals}
\\caption{Posterior predictive model diagnostic -- $(Tdesc[i]).}
      \\begin{tabular}{r||c|ccc}
        \\hline
        Gene ID & Mean expression (counts) & \$T(S_n)\$ & 95th int. \$T(x^{rep})\$ & P-value \\\\
        \\hline""")
   for j=1:rounds*D
       str = join([
        luad_norm_geneids[genes[j]],
        (@sprintf "%.1f" gnorms[genes[j]]),
        (@sprintf "%.2f" true_ts[j,i]),
        (@sprintf "[%.2f, %.2f]" ints[j,i,1] ints[j,i,2]),
        (@sprintf "%.2f" pvals[j,i])], " & ") * " \\\\"
      println(str)
    end
    println("""
\\end{tabular}
\\end{table*}
    """)
end


text:
module ParametricFunctions

using ForwardDiff
using ChainRulesCore

export ParametricFun, BasisFun
export Chebyshev, Fourier

export Parametric, NonParametric

export to_grad, from_grad, params

"""
    ParametricFun

An abstract type representing a parametric function.
"""
abstract type ParametricFun end

include("basis_fun.jl")
include("chebyshev_basis.jl")
include("fourier_basis.jl")

"""
    Parametric

A Trait type denoting functions with trainable parameters.

Functions implementing this trait must have methods for `has_params`, `to_grad`, `from_grad` and `params`.
"""
struct Parametric end

"""
    NonParametric

A trait type denoting functions with no parameters.
"""
struct NonParametric end

"""
    has_params(f)

Whether a functor type has trainable parameters.

Override this to return `Parametric()` for your type to declare it as parametric functions.
"""
has_params(f) = NonParametric()
has_params(::ParametricFun) = Parametric()

"""
    to_grad(f, x)

Return the gradient of `f` relative to its parameters at `x`.

The returned object must be an array (or any object that can be used in  a `QuadratureProblem`).
"""
function to_grad end

"""
        from_grad(f, ∇f)

Create a new instance of the type of `f` from a gradient over its parameters.
"""
function from_grad end

"""
    params(f)

Return the parameters of `f`.
"""
function params end

end


text:
module Api

using Genie, Logging, LoggingExtras

function main()
  Base.eval(Main, :(const UserApp = Api))

  Genie.genie(; context = @__MODULE__)

  Base.eval(Main, :(const Genie = Api.Genie))
  Base.eval(Main, :(using Genie))
end

end


text:
export marginalrule

@marginalrule GammaShapeRate(:out_α_θ) (m_out::Gamma, m_α::PointMass, m_β::PointMass) = begin
    return (out = prod(ProdAnalytical(), GammaShapeRate(mean(m_α), mean(m_β)), m_out), α = m_α, β = m_β)
end


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

name = "Keyutils"
version = v"1.6.1"

# Collection of sources required to build keyutils
sources = [
    "https://git.kernel.org/pub/scm/linux/kernel/git/dhowells/keyutils.git/snapshot/keyutils-$(version).tar.gz" =>
    "3c71dcfc6900d07b02f4e061d8fb218a4ae6519c1d283d6a57b8e27718e2f557",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/keyutils-*/
make -j${nproc}
make install \
    DESTDIR="${prefix}/" \
    SBINDIR="bin" \
    BINDIR="bin" \
    LIBDIR="lib" \
    USRLIBDIR="lib" \
    INCLUDEDIR="include" \
    MANDIR="share/man" \
    SHAREDIR="share/keyutils"

# Fix broken symlink
ln -sf libkeyutils.so.1 ${prefix}/lib/libkeyutils.so
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We are manually disabling
# many platforms that do not seem to work.
platforms = [p for p in supported_platforms() if Sys.islinux(p)]

# The products that we will ensure are always built
products = [
    LibraryProduct("libkeyutils", :libkeyutils),
]

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

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
"""
Returns whether it is applicable and the matched wilds
A truth can match with an expr in many ways
"""
# Returns whether it is applicable and the matched wilds
# a truth can match with an expr in many ways
function _isapplicable(expr::Term, truth::Term, wilds::Array{Term})::Tuple{Bool,Dict{Term,Set{Term}}}
    _result = true
    matched = Dict(wild => Set{Term}() for wild in wilds)  # where the wild needs to take on this value and this value...

    # println(expr)
    # println(truth)
    # println(codomain(expr))
    # println(codomain(truth))
    # println(codomain(expr) ⊂ codomain(truth))
    # println(result(codomain(expr) ⊂ codomain(truth)))
    # println()
    #is_subset = (result(codomain(expr) ⊂ codomain(truth)) == True)
    if (class(truth) <: AVariable) #&& is_subset
        _result = true
        push!(matched[truth], expr)
    elseif expr.class != truth.class
        _result = false
    elseif expr.data != truth.data
        _result = false
    elseif length(expr.args) != length(truth.args)
        _result = false
    else
        for i in 1:length(expr.args)
            inner_result, inner_matched = _isapplicable(expr.args[i], truth.args[i], wilds)
            for wild in wilds
                matched[wild] = union(matched[wild], inner_matched[wild])
            end
            if !inner_result
                _result = false
            end
        end
    end

    for wild in wilds
        if length(matched[wild]) > 1
            _result = false
        end
    end

    return _result, matched
end

function replace(expr::Term, old::Term, new::Term)::Term
    if expr == old
        return new
    else
        new_args = Tuple(replace(arg, old, new) for arg in args(expr))
        return Term(data(expr), class(expr), new_args)
    end
end

function replace(expr::Term, replacements::Dict{Term,Term})::Term
    result = expr
    for (key, value) in replacements
        result = replace(result, key, value)
    end
    return result
end

function get_wilds(expr::Term)::Array{Term}
    result = Term[]

    nodes = preorder_traversal(expr)
    for node in nodes
        if class(node) <: AVariable
            push!(result, node)
        end
    end

    return unique(result)
end

function isapplicable(expr::Term, truth::Term)
    wilds = get_wilds(truth)
    result, matched = _isapplicable(expr, truth, wilds)

    new_matched = Dict{Term,Term}()
    if result
        for wild in wilds
            new_matched[wild] = pop!(matched[wild])
        end
    end

    return result, new_matched
end
export isapplicable

function rewrites(expr::Term, truth_table=truths)::Array{Term}
    result = Term[]

    for name_truth in truth_table
        truth = name_truth.statement
        if class(truth) <: AImplies
            subject = args(truth)[1]
            predicate = args(truth)[2]
            app, matches = isapplicable(expr, subject)
            if app
                subs_right = replace(predicate, matches)
                push!(result, subs_right)
            end
        elseif class(truth) <: AEquivalent
            subject = args(truth)[1]
            predicate = args(truth)[2]
            app, matches = isapplicable(expr, subject)
            if app
                subs_right = replace(predicate, matches)
                push!(result, subs_right)
            end
            subject = args(truth)[2]
            predicate = args(truth)[1]
            app, matches = isapplicable(expr, subject)
            if app
                subs_right = replace(predicate, matches)
                push!(result, subs_right)
            end
        end
    end

    # Repeat recursivley on the arguments
    Threads.@threads for (i, child) in collect(enumerate(expr.args))

        child_rewrites = rewrites(child, truth_table)
        for j in 1:length(child_rewrites)
            # args = copy(expr.args)
            args::Array{Term{<:ABasic}} = collect(expr.args)
            args[i] = child_rewrites[j]
            child_rewrites[j] = Term(expr.data, expr.class, Tuple(args))
        end
        append!(result, child_rewrites)
    end

    return unique(result)
end
export rewrites


text:
# ---
# title: 1574. Shortest Subarray to be Removed to Make Array Sorted
# id: problem1574
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Binary Search
# link: <https://leetcode.com/problems/shortest-subarray-to-be-removed-to-make-array-sorted/description/>
# hidden: true
# ---
#
# Given an integer array `arr`, remove a subarray (can be empty) from `arr` such
# that the remaining elements in `arr` are **non-decreasing**.
#
# A subarray is a contiguous subsequence of the array.
#
# Return  _the length of the shortest subarray to remove_.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,2,3,10,4,2,3,5]
#     Output: 3
#     Explanation: The shortest subarray we can remove is [10,4,2] of length 3. The remaining elements after that will be [1,2,3,3,5] which are sorted.
#     Another correct solution is to remove the subarray [3,10,4].
#
# **Example 2:**
#
#
#
#     Input: arr = [5,4,3,2,1]
#     Output: 4
#     Explanation: Since the array is strictly decreasing, we can only keep a single element. Therefore we need to remove a subarray of length 4, either [5,4,3,2] or [4,3,2,1].
#
#
# **Example 3:**
#
#
#
#     Input: arr = [1,2,3]
#     Output: 0
#     Explanation: The array is already non-decreasing. We do not need to remove any elements.
#
#
# **Example 4:**
#
#
#
#     Input: arr = [1]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `0 <= arr[i] <= 10^9`
#
#
## @lc code=start
using LeetCode

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


text:
module Geo3x3
  function encode(lat, lng, level)
    if level < 1
      return null
    end
    res = ""
    if lng >= 0
      res = "E"
    else
      res = "W"
      lng += 180.0
    end
    lat += 90.0 # 180:the North Pole,  0:the South Pole
    unit = 180.0
    for i = 1:level - 1
      unit /= 3
      x = Int(floor(lng / unit))
      y = Int(floor(lat / unit))
      res = string(res, x + y * 3 + 1)
      lng -= x * unit
      lat -= y * unit
    end
    return res
  end

  function decode(code)
    clen = length(code)
    if clen == 0
      return null
    end
    beg = 0
    flg = false
    c = code[1]
    if c == '-' || c == 'W'
      flg = true
      beg = 1
    elseif c == '+' || c == 'E'
      beg = 1
    end
    unit = 180.0
    lat = 0.0
    lng = 0.0
    level = 1
    for i = beg:clen - 1
      n = Int(code[i + 1]) - 48
      if n < 1 || n > 9
        break
      end
      unit /= 3.0
      n -= 1
      lng += n % 3 * unit
      lat += floor(n / 3) * unit
      level += 1
    end
    lat += unit / 2.0
    lng += unit / 2.0
    lat -= 90.0
    if flg
      lng -= 180.0
    end
    return (lat, lng, level, unit)
  end
end


text:
# U: the dimension of the universe
# D: the dimension of the manifold
# N: the number of vertices
# T: the type of the coordinates
# C: the complimentary dimension (should always be U-D)
struct SphereChart{P,T}
    center::P
    radius::T
end

volume(p::SphereChart) = 4π * p.radius^2

coordtype(ch::Type{SphereChart{P,T}}) where {P,T} = T
coordtype(p::SphereChart) = coordtype(typeof(p))

dimension(p::Type{SphereChart{P,T}}) where {P,T} = 2
dimension(p::SphereChart) = dimension(typeof(p))

universedimension(::Type{SphereChart{P,T}}) where {P,T} = 3
universedimension(p::SphereChart) = universedimension(typeof(p))

struct ThetaPhiPlane{T} end
domain(ch::SphereChart{P,T}) where {P,T} = ThetaPhiPlane{T}()

struct SphereNeighborhood{P,Q,T}
    sphere::SphereChart{P,T}
    paramters::Q
    cartesian::P
    tangents::Tuple{P,P}
end

function neighborhood(ch::SphereChart, uv)
    ϕ, Θ = uv

    cosϕ, sinϕ = cos(ϕ), sin(ϕ)
    cosΘ, sinΘ = cos(Θ), sin(Θ)

    p = @SVector [cosϕ*sinΘ, sinϕ*sinΘ, cosΘ]

    tΘ = @SVector [cosϕ*cosΘ, sinϕ*cosΘ, -sinΘ]
    tϕ = @SVector [-sinϕ*sinΘ/abs(sinΘ), cosϕ*sinΘ/abs(sinΘ), 0]

    SphereNeighborhood(ch, uv, p, (tΘ, tϕ))
end

cartesian(nbd::SphereNeighborhood) = nbd.cartesian

function quadpoints(ch::SphereChart, L)
    UW = zip(legendre(L+1, -1.0, 1.0)...)
    dp = 2π / (2L+1)
    #P = linspace(0, 2π-dp, 2L+1)
    P = range(0, stop=2π-dp, length=2L+1)
    [(neighborhood(ch, (p,acos(u))), dp*w) for (u,w) in UW for p in P]
end


text:
#############################################################################
#############################################################################
#
# This file implements polynomial division
#
#############################################################################
#############################################################################

"""  Modular algorithm.
f divide by g

f = q*g + r

p is a prime
"""

function divide(num::Polynomial, den::Polynomial)
    function division_function(p::Int)
        f, g = mod(num,p), mod(den,p)
        degree(f) < degree(num) && return nothing
        iszero(g) && throw(DivideError())
        q = Polynomial()
        prev_degree = degree(f)
        while degree(f) ≥ degree(g)
            h = Polynomial( (leading(f) ÷ leading(g))(p) )  #syzergy
            f = mod((f - h*g), p)
            q = mod((q + h), p)
            prev_degree == degree(f) && break
            prev_degree = degree(f)
        end

        @assert iszero( mod((num  - (q*g + f)),p))
        return q, f
    end
    return division_function
end

function divide(num::PolynomialModP, den::PolynomialModP)
    p = num.mod
    @assert num.mod == den.mod
    f, g = num, den
    if degree(g) > degree(f)
        output =  [zero(PolynomialModP,p) , f]
        return output
    end

    iszero(g) && throw(DivideError())
    iszero(f) && return zero(PolynomialModP,p), zero(PolynomialModP,p)

    q = PolynomialModP(f.mod)
    prev_degree = degree(f)
    while degree(f) ≥ degree(g)
        h = (leading(f) ÷ leading(g))(p) #syzergy
        hg = - h*g
        f = f+hg
        q = q+h
        prev_degree == degree(f) && break
        prev_degree = degree(f)
    end

    @assert iszero(num - (q*g + f))

    return q,f
end

"""
The quotient from polynomial division. Returns a function of an integer.
"""
÷(num::Polynomial, den::Polynomial)  = (p::Int) -> first(divide(num,den)(p))
÷(num::PolynomialModP, den::PolynomialModP) = first(divide(num,den))

"""
The remainder from polynomial division. Returns a function of an integer.
"""
rem(num::Polynomial, den::Polynomial)  = (p::Int) -> last(divide(num,den)(p))
rem(num::PolynomialModP, den::PolynomialModP)  = last(divide(num,den))


text:
fn = "0002.sdf"
v_header, data, grids, units = open(deserialize, "0002.jl")

@testset "SDF header" begin
    h = header(fn)
    saved_props = setdiff(propertynames(h),
        [:sdf_magic, :endianness, :summary_location, :summary_size,
        :first_block_location, :nblocks, :block_header_length,
        :string_length])
    for p in saved_props
        @test getproperty(h, p) == getindex(v_header, string(p))
    end
end


text:
struct Ensemble{D} <: MHSampler
    n_walkers::Int
    proposal::D
end

# Define the first sampling step.
# Return a 2-tuple consisting of the initial sample and the initial state.
# In this case they are identical.
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::DensityModel,
    spl::Ensemble;
    init_params = nothing,
    kwargs...,
)
    if init_params === nothing
        transitions = propose(rng, spl, model)
    else
        transitions = [Transition(model, x) for x in init_params]
    end

    return transitions, transitions
end

# Define the other sampling steps.
# Return a 2-tuple consisting of the next sample and the the next state.
# In this case they are identical, and for each walker they are either a new proposal
# (if accepted) or the previous proposal (if not accepted).
function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::DensityModel,
    spl::Ensemble,
    params_prev::Vector{<:Transition};
    kwargs...,
)
    # Generate a new proposal. Accept/reject happens at proposal level.
    transitions = propose(rng, spl, model, params_prev)
    return transitions, transitions
end

#
# Initial proposal
#
function propose(rng::Random.AbstractRNG, spl::Ensemble, model::DensityModel)
    # Make the first proposal with a static draw from the prior.
    static_prop = StaticProposal(spl.proposal.proposal)
    mh_spl = MetropolisHastings(static_prop)
    return [propose(rng, mh_spl, model) for _ in 1:spl.n_walkers]
end

#
# Every other proposal
#
function propose(
    rng::Random.AbstractRNG,
    spl::Ensemble,
    model::DensityModel,
    walkers::Vector{<:Transition},
)
    new_walkers = similar(walkers)

    others = 1:(spl.n_walkers - 1)
    for i in 1:spl.n_walkers
        walker = walkers[i]
        idx = mod1(i + rand(rng, others), spl.n_walkers)
        other_walker = walkers[idx]
        new_walkers[i] = move(rng, spl, model, walker, other_walker)
    end

    return new_walkers
end

#####################################
# Basic stretch move implementation #
#####################################
struct StretchProposal{P, F<:AbstractFloat} <: Proposal{P}
    proposal :: P
    stretch_length::F
end

StretchProposal(p) = StretchProposal(p, 2.0)

function move(
    rng::Random.AbstractRNG,
    spl::Ensemble{<:StretchProposal},
    model::DensityModel,
    walker::Transition,
    other_walker::Transition,
)
    # Calculate intermediate values
    proposal = spl.proposal
    n = length(walker.params)
    a = proposal.stretch_length
    z = ((a - 1) * rand(rng) + 1)^2 / a
    alphamult = (n - 1) * log(z)

    # Make new parameters
    y = @. other_walker.params + z * (walker.params - other_walker.params)

    # Construct a new walker
    new_walker = Transition(model, y)

    # Calculate accept/reject value.
    alpha = alphamult + new_walker.lp - walker.lp

    if -Random.randexp(rng) <= alpha
        return new_walker
    else
        return walker
    end
end

#########################
# Elliptical slice step #
# #########################

# struct EllipticalSlice{E} <: ProposalStyle
#     ellipse::E
# end

# function move(
#     # spl::Ensemble,
#     spl::Ensemble{Proposal{T,P}},
#     model::DensityModel,
#     walker::Transition,
#     other_walker::Transition,
# ) where {T<:EllipticalSlice,P}
#     # Calculate intermediate values
#     proposal = spl.proposal
#     n = length(walker.params)
#     nu = rand(proposal.type.ellipse)

#     u = rand()
#     y = walker.lp - Random.randexp()

#     theta = 2 * π * rand()

#     theta_min = theta - 2.0*π
#     theta_max = theta

#     f = walker.params
#     while true
#         stheta, ctheta = sincos(theta)

#         f_prime = f .* ctheta + nu .* stheta

#         new_walker = Transition(model, f_prime)

#         if new_walker.lp > y
#             return new_walker
#         else
#             if theta < 0
#                 theta_min = theta
#             else
#                 theta_max = theta
#             end

#             theta = theta_min + (theta_max - theta_min) * rand()
#         end
#     end
# end

#####################
# Slice and stretch #
#####################
# struct EllipticalSliceStretch{E, S<:Stretch} <: ProposalStyle
#     ellipse::E
#     stretch::S
# end

# EllipticalSliceStretch(e) = EllipticalSliceStretch(e, Stretch(2.0))

# function move(
#     # spl::Ensemble,
#     spl::Ensemble{Proposal{T,P}},
#     model::DensityModel,
#     walker::Transition,
#     other_walker::Transition,
# ) where {T<:EllipticalSliceStretch,P}
#     # Calculate intermediate values
#     proposal = spl.proposal
#     n = length(walker.params)
#     nu = rand(proposal.type.ellipse)

#     # Calculate stretch step first
#     subspl = Ensemble(spl.n_walkers, Proposal(proposal.type.stretch, proposal.proposal))
#     walker = move(subspl, model, walker, other_walker)

#     u = rand()
#     y = walker.lp - Random.randexp()

#     theta = 2 * π * rand()

#     theta_min = theta - 2.0*π
#     theta_max = theta

#     f = walker.params

#     i = 0
#     while true
#         i += 1

#         stheta, ctheta = sincos(theta)

#         f_prime = f .* ctheta + nu .* stheta

#         new_walker = Transition(model, f_prime)

#         # @info "Slice step" i f f_prime y new_walker.lp theta theta_max theta_min

#         if new_walker.lp > y
#             return new_walker
#         else
#             if theta < 0
#                 theta_min = theta
#             else
#                 theta_max = theta
#             end

#             theta = theta_min + (theta_max - theta_min) * rand()
#         end
#     end
# end


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

name = "OpenSlide"
version = v"3.4.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/openslide/openslide/releases/download/v3.4.1/openslide-3.4.1.tar.gz", "fed08fab8a9b1ded95a34e196652291127ebe392c11f9bc13d26e760295a102d")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd openslide-3.4.1/
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = filter(!Sys.isapple, supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libopenslide", :libopenslide)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="Zlib_jll", uuid="83775a58-1f1d-513f-b197-d71354ab007a"))
    Dependency(PackageSpec(name="libpng_jll", uuid="b53b4c65-9356-5827-b1ea-8c7a1a84506f"))
    Dependency(PackageSpec(name="JpegTurbo_jll", uuid="aacddb02-875f-59d6-b918-886e6ef4fbf8"))
    Dependency(PackageSpec(name="Libtiff_jll", uuid="89763e89-9b03-5906-acba-b20f662cd828"))
    Dependency(PackageSpec(name="OpenJpeg_jll", uuid="643b3616-a352-519d-856d-80112ee9badc"))
    Dependency(PackageSpec(name="gdk_pixbuf_jll", uuid="da03df04-f53b-5353-a52f-6a8b0620ced0"))
    Dependency(PackageSpec(name="XML2_jll", uuid="02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"))
    Dependency(PackageSpec(name="SQLite_jll", uuid="76ed43ae-9a5d-5a62-8c75-30186b810ce8"))
    Dependency(PackageSpec(name="Cairo_jll", uuid="83423d85-b0ee-5818-9007-b63ccbeb887a"))
    Dependency(PackageSpec(name="Glib_jll", uuid="7746bdde-850d-59dc-9ae8-88ece973131d"), v"2.59.0"; compat="2.59")
    BuildDependency(PackageSpec(name="Xorg_xorgproto_jll", uuid="c4d99508-4286-5418-9131-c86396af500b"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
module TestGraphGenerators
using Test

using Catlab.Graphs.BasicGraphs, Catlab.Graphs.GraphGenerators
using Catlab.CategoricalAlgebra: homomorphisms

# Path graphs
#------------

n = 5
g = path_graph(Graph, n)
@test (nv(g), ne(g)) == (n, n-1)
g = path_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 2(n-1))

# Cycle graphs
#-------------

g = cycle_graph(Graph, n)
@test (nv(g), ne(g)) == (n, n)
g = cycle_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 2n)

# Complete graphs
#----------------

for T in (Graph, SymmetricGraph)
  g = complete_graph(T, n)
  @test (nv(g), ne(g)) == (n, n*(n-1))
end
for T in (ReflexiveGraph, SymmetricReflexiveGraph)
  g = complete_graph(T, n)
  @test (nv(g), ne(g)) == (n, n*n)
end

# Star graphs
#------------

g = star_graph(Graph, n)
@test (nv(g), ne(g)) == (n, n-1)
@test length(unique(src(g,e) for e in edges(g))) == 1
g = star_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 2(n-1))

# Wheel graphs
#-------------

g = wheel_graph(Graph, n)
@test (nv(g), ne(g)) == (n, 2(n-1))
triangle = Graph(3)
add_edges!(triangle, [1,2,1], [2,3,3])
@test length(homomorphisms(triangle, g)) == n-1

g = wheel_graph(SymmetricGraph, n)
@test (nv(g), ne(g)) == (n, 4(n-1))
triangle = cycle_graph(SymmetricGraph, 3)
@test length(homomorphisms(triangle, g)) == 6(n-1) # == 3! * (n-1)

# Parallel arrows
#----------------

g = parallel_arrows(Graph, n)
@test (nv(g), ne(g)) == (2, n)
@test all(==(1), src(g))
@test all(==(2), tgt(g))

# Erdos-Renyi graphs
#-------------------

g = erdos_renyi(Graph, 100, 0.0)
@test (nv(g), ne(g)) == (100, 0)
g = erdos_renyi(Graph, 100, 0)
@test (nv(g), ne(g)) == (100, 0)
g = erdos_renyi(Graph, 100, 5)
@test (nv(g), ne(g)) == (100, 5)
g = erdos_renyi(Graph, 100, 1.0)
@test g == complete_graph(Graph, 100)

# Expected Degree Graph
#----------------------

g = expected_degree_graph(Graph, 5 .* ones(100))
@test nv(g) == 100

# Watts Strogatz Graph
#---------------------

g = watts_strogatz(Graph, 100, 4, 0.2)
@test (nv(g), ne(g)) == (100, 200)

end

max_forks_repo_name: epatters/CompCat


text:
using OrdinaryDiffEq, Test, Random
Random.seed!(123)

using OrdinaryDiffEq, DiffEqOperators, LinearAlgebra
A = 0.01*rand(3, 3)
rn = (du, u, p, t) -> begin
    mul!(du, A, u)
end
u0 = rand(3)
prob = ODEProblem(ODEFunction(rn, jac_prototype=JacVecOperator{Float64}(rn, u0; autodiff=false)), u0, (0, 10.))
@test_nowarn sol = solve(prob, TRBDF2(autodiff=false));
@test_nowarn sol = solve(prob, TRBDF2(autodiff=false, linsolve=LinSolveGMRES()));
@test_nowarn sol = solve(prob, TRBDF2(autodiff=false, linsolve=LinSolveGMRES(), smooth_est=false));
@test_nowarn sol = solve(prob, TRBDF2(autodiff=false, linsolve=LinSolveGMRES(Pl=lu(A)), smooth_est=false));
@test_nowarn sol = solve(prob, TRBDF2(autodiff=false, linsolve=LinSolveGMRES(Pr=lu(A)), smooth_est=false));
@test_nowarn sol = solve(prob, TRBDF2(autodiff=false, linsolve=LinSolveGMRES(Pl=lu(A), Pr=lu(A)), smooth_est=false));


text:
import Dates

struct HolidayRoundingNBD <: HolidayRoundingConvention end
function apply(::HolidayRoundingNBD, date::Dates.Date, calendar::Calendar)::Dates.Date
    while is_holiday(calendar, date)
        date += Dates.Day(1)
    end
    date
end
Base.show(io::IO, ::HolidayRoundingNBD) = print(io, "NBD")

struct HolidayRoundingPBD <: HolidayRoundingConvention end
function apply(::HolidayRoundingPBD, date::Dates.Date, calendar::Calendar)::Dates.Date
    while is_holiday(calendar, date)
        date -= Dates.Day(1)
    end
    date
end
Base.show(io::IO, ::HolidayRoundingPBD) = print(io, "PBD")

struct HolidayRoundingNBDSM <: HolidayRoundingConvention end
function apply(::HolidayRoundingNBDSM, date::Dates.Date, calendar::Calendar)::Dates.Date
    new_date = date
    while is_holiday(calendar, new_date)
        new_date += Dates.Day(1)
    end

    if Dates.month(new_date) != Dates.month(date)
        new_date = date
        while is_holiday(calendar, new_date)
            new_date -= Dates.Day(1)
        end
    end

    new_date
end
Base.show(io::IO, ::HolidayRoundingNBDSM) = print(io, "NBDSM")

struct HolidayRoundingPBDSM <: HolidayRoundingConvention end
function apply(::HolidayRoundingPBDSM, date::Dates.Date, calendar::Calendar)::Dates.Date
    new_date = date
    while is_holiday(calendar, new_date)
        new_date -= Dates.Day(1)
    end

    if Dates.month(new_date) != Dates.month(date)
        new_date = date
        while is_holiday(calendar, new_date)
            new_date += Dates.Day(1)
        end
    end

    new_date
end
Base.show(io::IO, ::HolidayRoundingPBDSM) = print(io, "PBDSM")

struct HolidayRoundingNR <: HolidayRoundingConvention end
apply(::HolidayRoundingNR, date::Dates.Date, ::Calendar) = date
Base.show(io::IO, ::HolidayRoundingNR) = print(io, "NR")

const HOLIDAY_ROUNDING_MAPPINGS = Dict("NR" => HolidayRoundingNR(), "NBD" => HolidayRoundingNBD(), "PBD" => HolidayRoundingPBD(), "NBDSM" => HolidayRoundingNBDSM(), "PBDSM" => HolidayRoundingPBDSM())


text:
using Test
using JSON3

using Pinecone

APIKEY = ENV["PINECONE_API_KEY"]
CLOUDENV="us-west1-gcp"

context = Pinecone.init(APIKEY, CLOUDENV)

INDEX = "filter-example"
NAMESPACE = "mynamespace"
index = Pinecone.Index(INDEX);

testdict = Dict{String, Any}("genre"=>"documentary", "year"=>2019);
v1 = [0.3,0.11,0.3,0.3,0.3,0.3,0.3,0.3,0.4,0.3]
v2 = [0.9, 0.8, 0.7, 0.6, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3]

@testset verbose = true "Test upsert()" begin
    testdict = Dict{String, Any}("genre"=>"documentary", "year"=>2019);
    meta = [Dict{String,Any}("foo"=>"bar"), Dict{String,Any}("bar"=>"baz")]
    testvector = Pinecone.PineconeVector("testid", v1, testdict)

    @testset "Regular upsert()" begin
       result = Pinecone.upsert(context, index, [testvector], "testnamespace")
       @test result !== nothing
       @test typeof(result) == String
       #test no namespace

       result = Pinecone.upsert(context, index, [testvector])

       @test result !== nothing
       @test typeof(result) == String
    end
    @testset "Upsert with invalid args (ArgumentError)" begin
       #test Arg checking exceptions
       @test_throws ArgumentError Pinecone.upsert(context, index, ["zipA"], [v1, v2])
       @test_throws ArgumentError Pinecone.upsert(context, index, ["zipA", "zipB"], [v2])
       @test_throws ArgumentError Pinecone.upsert(context, index, ["zipA", "zipB"], [v1, v2], [Dict{String,Any}("foo"=>"bar")])
       # test too many vectors on insert
       @test_throws ArgumentError Pinecone.upsert(context, index, ["zipA", "zipB"], [rand(Float32,10) for i in 1:1001])
       largevec = [testvector for i in 1:1001]
       @test_throws ArgumentError Pinecone.upsert(context, index, largevec)
       @test_throws ArgumentError Pinecone.upsert(context, index, largevec, "ns")
    end

    @testset "Upsert with zipped vectors and ids" begin
       result = Pinecone.upsert(context, index, ["zipA", "zipB"], [v1, v2])
       @test result !== nothing
       @test typeof(result) == String

       result = Pinecone.upsert(context, index, ["zipA", "zipB"], [v1, v2], meta)
       @test result !== nothing
       @test typeof(result) == String

       result = Pinecone.upsert(context, index, ["zipA", "zipB"], [v1, v2], meta, "mynamespace")
       @test result !== nothing
       @test typeof(result) == String
    end

    @testset "Large upsert()" begin
       context = Pinecone.init(APIKEY, CLOUDENV)
       largeindex = Pinecone.Index("testlargeindex")
       testrows= Vector{Pinecone.PineconeVector}()
       for i in 1:1000
          testvector = Pinecone.PineconeVector("testid_$i", rand(Float32, 10), testdict)
          push!(testrows, testvector)
       end
       result = Pinecone.upsert(context, largeindex, testrows, "testnamespace")
       @test result !== nothing
       @test typeof(result) == String
    end
end


text:
figure_directory = "figures/"

# Discontinuous Galerkin Intro visualizations
include("dg_projection.jl")
save(figure_directory * "DG_projections.png", fig)
include("dg_derivative.jl")
save(figure_directory * "DG_derivatives.png", fig)
include("dg_split_form_tendencies.jl")
save(figure_directory * "DG_split_form_tendencies.png", fig)

# Held-Suarez visualizations
include("held_suarez_statistics.jl")
save(figure_directory * "Held_Suarez_Statistics.png", fig)

include("held_suarez_resolution.jl")
save(figure_directory * "Held_Suarez_Resolutions.png", fig)

include("small_held_suarez_statistics.jl")
save(figure_directory * "Small_Held_Suarez_Statistics.png", fig)

include("held_suarez_snapshots.jl")
save(figure_directory * "surface_fields.png", fig)

include("held_suarez_3d_snapshots.jl")
save(figure_directory * "vertical_velocity.png", fig)

include("annulus_plot.jl")
save(figure_directory * "annulus_and_gridpoints.png", fig)

include("cubed_sphere_warp.jl")
save(figure_directory * "cubed_sphere.png", fig)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
module MaterialModels

using Polynomials
using LinearAlgebra
using Dierckx

using ..Geometry
using ..CrossSection
using ..Mesh

export steel, ConstitutiveLaw

struct ConstitutiveLaw

    stress_strain::Dierckx.Spline1D
    E_tan_strain::Dierckx.Spline1D

end

function steel(σy, σy1, σu, σf, ϵy, ϵy1, ϵu, ϵf, n)

    # define elastic and yielding portion of the steel stress-strain curve

    seg_1 = norm([σy, ϵy])
    seg_2 = norm([(σy-σy1), (ϵy-ϵy1)])

    θ1 = rad2deg(atan(σy/ϵy))
    θ2 = rad2deg(atan((σy-σy1) / (ϵy-ϵy1)))

    ΔL = [seg_1, seg_2]
    θ = [θ1, θ2]  #degrees
    n_ep = [n[1], n[2]]          #hard code this for now
    radius = [0.001]  #hard code this for now
    n_radius = [400]     #hard code this for now
    closed_or_open = 1

   elastic_plastic = CrossSection.Feature(ΔL, θ, n_ep, radius, n_radius, closed_or_open)
   ϵ, σ = CrossSection.get_xy_coordinates(elastic_plastic)

    # #define plateau region
    # plateau_slope = (σy1 - σy) / (ϵy1 - ϵy)
    # ϵ_plateau = (ϵy + (ϵy1 - ϵy)/n[2]):(ϵy1 - ϵy)/n[2]:(ϵy1 - (ϵy1 - ϵy)/n[2])
    # σ_plateau = plateau_slope .* (ϵ_plateau .- ϵy) .+ σy

    # ϵ = [ϵ; ϵ_plateau]
    # σ = [σ; σ_plateau]

    #define plateau and hardening part of steel stress-strain curve
    p = fit([ϵy1, ϵu, ϵf],[σy1, σu, σf], 2)

    ϵ_range = ϵy1:(ϵf-ϵy1)./n[3]:ϵf

    σ_fit = zeros(length(ϵ_range))

    for i=1:length(ϵ_range)
        σ_fit[i] = p(ϵ_range[i])
    end

    ϵ = [ϵ; ϵ_range[2:end]]
    σ = [σ; σ_fit[2:end]]

    #calculate tangent elastic modulus
    num_points = length(ϵ)
    E_tan = zeros(Float64,  num_points)
    order = 1
    for i = 1:num_points

        if i < num_points - 2

            x = ϵ[i:i+2]
            x0 = x[1]

            stencil = Mesh.calculate_weights(order, x0, x)

            E_tan[i] = sum(stencil .* σ[i:i+2])

        elseif i > num_points - 2

            x = ϵ[i-2:i]
            x0 = x[3]

            stencil = Mesh.calculate_weights(order, x0, x)

            E_tan[i] = sum(stencil .* σ[i-2:i])

        end

        #tough to round the corner at yielding, revisit this with better model of proportional limit
        if E_tan[i] > E_tan[1]

            E_tan[i] = E_tan[i-1]

        end

    end

    stress_strain = Spline1D(ϵ, σ)
    E_tan_strain = Spline1D(ϵ, E_tan)

    steel_law = ConstitutiveLaw(stress_strain, E_tan_strain)

    return steel_law

end

end


text:
"helper for testing"
function quick_plot(graph, xs, ys)
    nv(graph)==length(xs) == length(ys) || error("need 1 position per vertex")
    scatter(xs, ys; markeralpha=0, text=string.(vertices(graph)))

    weights_mat = weights(graph)
    # now draw connections
    for edge in edges(graph)
        lxs = [xs[src(edge)], xs[dst(edge)]]
        lys = [ys[src(edge)], ys[dst(edge)]]
        w = 5*weights_mat[edge.src, edge.dst]
        plot!(lxs, lys; linewidth=w, alpha=0.7, legend=false)
    end
end

function quick_plot(graph, xs, ys, paths)
    nv(graph)==length(xs) == length(ys) || error("need 1 position per vertex")
    scatter(xs, ys; markeralpha=0, text=string.(vertices(graph)))

    weights_mat = weights(graph)
    # now draw connections
    for edge in edges(graph)
        lxs, lys = paths[edge]
        w = 5*weights_mat[edge.src, edge.dst]
        plot!(lxs, lys; linewidth=w, alpha=0.7, legend=false)
    end
end

function quick_plot_solve_paths(layout, graph)
    xs, ys, paths = solve_positions(layout, graph)
    quick_plot(graph, xs, ys, paths)
end
function quick_plot_solve_direct(layout, graph)
    xs, ys, _ = solve_positions(layout, graph)
    quick_plot(graph, xs, ys)
end

@testset "quick_plot" begin
    ref_filename =  joinpath(@__DIR__, "references", "test_utils", "$quick_plot.png")
    @plottest quick_plot(SimpleDiGraph(Edge.([1=>2, 2=>3])), [1,2,5], [1,2,3]) ref_filename true 0.05
    paths = Dict(Edge(1, 2) => ([1, 2], [1, 2]), Edge(2, 3) => ([2, 5], [2, 3]))
    @plottest quick_plot(SimpleDiGraph(Edge.([1=>2, 2=>3])), [1,2,5], [1,2,3], paths) ref_filename true 0.05
end

function test_example(layout, graph_name, tol=0.05)
    @testset "$graph_name" begin
        @testset "$graph_name direct" begin
            graph = getfield(Examples, graph_name)
            ref_filename = joinpath(@__DIR__, "references", string(typeof(layout)), "direct", "$graph_name.png")
            mkpath(dirname(ref_filename))
            @plottest quick_plot_solve_direct(layout, graph) ref_filename true tol
        end
        @testset "$graph_name paths" begin
            graph = getfield(Examples, graph_name)
            ref_filename = joinpath(@__DIR__, "references", string(typeof(layout)), "paths", "$graph_name.png")
            mkpath(dirname(ref_filename))
            @plottest quick_plot_solve_paths(layout, graph) ref_filename true tol
        end
    end
end

@testset "$layout Demos" for layout in (Zarate(),)
    test_example(layout, :cross)
    test_example(layout, :loop)
    test_example(layout, :medium_pert)
    test_example(layout, :sankey_3twos)
    test_example(layout, :two_lines, 0.02)
    test_example(layout, :xcross)
    test_example(layout, :tree, 0.07)
    #test_example(layout, :large_depgraph)  # too big
    #test_example(layout, :extra_large_depgraph)  # too big
end


text:
using Documenter, MonteCarlo

makedocs(
    modules = [MonteCarlo],
    format = :html,
    sitename = "MonteCarlo.jl",
    pages = [
        "Home" => "index.md",
        "Manual" => [
            "Getting started" => "manual/gettingstarted.md",
            "Examples" => "manual/examples.md"
        ],
        "Physical models" => [
            "Ising model" => "models/ising.md",
            "Attractive Hubbard model" => "models/hubbardattractive.md"
        ],
        "Monte Carlo flavors" => [
            "MC" => "flavors/mc.md",
            "DQMC" => "flavors/dqmc.md"
        ],
        "Lattices" => "lattices.md",
        "Customize" => "customize.md",
        "Methods" => [
            "General" => "methods/general.md"
        ],
        "Interfaces" => [
            "MC" => "interfaces/MC.md",
            "DQMC" => "interfaces/DQMC.md"
        ]
    ]
    # assets = ["assets/custom.css", "assets/custom.js"]
)

deploydocs(
    repo   = "github.com/crstnbr/MonteCarlo.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    julia  = "release",
    osname = "linux"
)


text:
"""
    PaulisI

This module contains `PauliI <: AbstractPauli` and supporting code.
"""
module PaulisI

import ..AbstractOps: op_index, _AbstractOp
import ..AbstractPaulis: AbstractPauli

export PauliI

"""
    struct PauliI <: AbstractPauli

Another implementation of `AbstractPauli`.

This implementation encodes the Pauli operators via the index `0-3`
in a single byte.
"""
struct PauliI <: AbstractPauli{Complex{Int}}
    ind::Int
end

####
#### Constructors
####

const I = PauliI(0)
const X = PauliI(1)
const Y = PauliI(2)
const Z = PauliI(3)

# const PAULIS = (I, X, Y, Z)

# Indexing into static Tuple seems faster
"""
    Pauli(ind::Union{Integer, Symbol, AbstractString, AbstractChar})

Return a `Pauli` indexed by `[0, 3]` or a representation of `[I, X, Y, Z]`.
"""
PauliI(s::Union{Symbol, AbstractString, AbstractChar}) = _AbstractOp(PauliI, s)

Base.copy(p::PauliI) = p

op_index(p::PauliI) = p.ind

####
#### Compare / predicates
####

Base.:(==)(p1::PauliI, p2::PauliI) = p1.ind == p2.ind

## This is maybe not really necessary, but this is stricter than the fallback method.
## Also, this would have prevented a perf regression bug.
Base.isone(p::PauliI) = p === one(PauliI)

####
#### Algebra / mathematical operations
####

## We try two versions of the multiplication table.
## Vector and matrix versions have similar performance.

const pauli_i_mult =
    ((0, 1, 2, 3),
     (1, 0, 3, 2),
     (2, 3, 0, 1),
     (3, 2, 1, 0))

const pauli_i_mult_1d =
    (0, 1, 2, 3,
     1, 0, 3, 2,
     2, 3, 0, 1,
     3, 2, 1, 0)

using ..Paulis

"""
    *(p1::PauliI, p2::PauliO)

Multiplication that returns a `PauliI`, but ignores the phase.
"""
Base.:*(p1::PauliI, p2::PauliI) = PauliI(@inbounds pauli_i_mult_1d[p1.ind * 4 + p2.ind + 1])

#Base.:*(p1::PauliI, p2::PauliI) = PauliI(pauli_i_mult[p1.ind+1][p2.ind+1])
## This one allocates for some reason
#Base.:*(p1::PauliI, p2::PauliI) = PauliI(op_index(Pauli(p1.ind) * Pauli(p2.ind)))

end # module PaulisI


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.

# TODO add doc
depends_on(ex::Expr, s::Symbol) = any(a -> depends_on(a, s), ex.args)
depends_on(ex::Symbol, s::Symbol) = (ex == s)
depends_on(ex, s::Symbol) = false
function depends_on(ex1, ex2)
    @assert isa(ex2, Expr)
    @assert ex2.head == :tuple
    return any(s -> depends_on(ex1, s), ex2.args)
end

function is_dependent(idxvars, idxset, i)
    for (it, idx) in enumerate(idxvars)
        it == i && continue
        depends_on(idxset, idx) && return true
    end
    return false
end

function has_dependent_sets(idxvars, idxsets)
    # check if any index set depends on a previous index var
    for i in 2:length(idxsets)
        for v in idxvars[1:(i-1)]
            if depends_on(idxsets[i], v)
                return true
            end
        end
    end
    return false
end

valid_array_index_set(::Base.OneTo) = true
valid_array_index_set(r::UnitRange) = first(r) == 1
valid_array_index_set(s) = false

"""
    generate_container(T, indexvars, indexsets, requestedtype)

Return a tuple, the first element of which is *code* that generates a container
for objects of type `T` given the index variables, index sets, and
`requestedtype`. `requestedtype` may be one of `:Array`, `:DenseAxisArray`,
`:SparseAxisArray`, or `:Auto`. Return error-producing code if requested type is
incompatible. For the case of `:Auto`, the following rules are used to determine
the appropriate container:

1. If all index sets are either explicit `1:B` objects for any `B` or symbols
   which refer to objects of type `Base.OneTo`, then an `Array` is generated of
   the appropriate size. Types of symbols/expressions are not known at compile
   time, so we defer to type-safe functions to check the `Base.OneTo` condition.

2. If condition (1) does not hold, and the index sets are independent (the index
   variable for one set does not appear in the definition of another), then an
   `DenseAxisArray` is generated of the appropriate size.

3. Otherwise, generate an empty `SparseAxisArray{T,N,NTuple{N,Any}}`.

The second element of the return tuple is a `Bool`, `true` if the container type
automatically checks for duplicate terms in the index sets and `false` otherwise.

### Examples

    generate_container(VariableRef, [:i,:j], [:(1:N), :(1:T)], :Auto)
    # Returns code equivalent to:
    # :(Array{VariableRef}(length(1:N), length(1:T))

    generate_container(VariableRef, [:i,:j], [:(1:N), :(2:T)], :Auto)
    # Returns code equivalent to:
    # :(JuMP.Containers.DenseAxisArray(undef, 1:N, 2:T))

    generate_container(VariableRef, [:i,:j], [:(1:N), :(S)], :Auto)
    # Returns code that generates an Array if S is of type Base.OneTo,
    # otherwise an DenseAxisArray.

    generate_container(VariableRef, [:i,:j], [:(1:N), :(1:j)], :Auto)
    # Returns code equivalent to:
    # :(Containers.SparseAxisArray(Dict{NTuple{N,Any},VariableRef}()))
"""
function generate_container(T, indexvars, indexsets, requestedtype)
    has_dependent = has_dependent_sets(indexvars, indexsets)
    onetosets = falses(length(indexsets))
    for (i, indexset) in enumerate(indexsets)
        s = Meta.isexpr(indexset, :escape) ? indexset.args[1] : indexset
        if Meta.isexpr(s, :call) &&
           length(s.args) == 3 &&
           s.args[1] == :(:) &&
           s.args[2] == 1
            onetosets[i] = true
        end
    end

    if requestedtype == :SparseAxisArray ||
       (requestedtype == :Auto && has_dependent)
        N = length(indexvars)
        @assert N == length(indexsets)
        return :(JuMP.Containers.SparseAxisArray(Dict{NTuple{$N,Any},$T}())),
        false
    end

    sizes = Expr(:tuple, [:(length($rng)) for rng in indexsets]...)
    arrayexpr = :(Array{$T}(undef, $sizes...))

    if requestedtype == :Array
        has_dependent && return :(error(
            "Unable to create requested Array because index sets are dependent.",
        )),
        true
        if all(onetosets)
            return arrayexpr, true
        else
            # all sets must be one-based intervals
            condition = Expr(:&&)
            for (i, indexset) in enumerate(indexsets)
                if !onetosets[i]
                    push!(
                        condition.args,
                        Expr(
                            :call,
                            :(JuMP.Containers.valid_array_index_set),
                            indexset,
                        ),
                    )
                end
            end
            return :(
                $condition ||
                    error("Index set for array is not one-based interval.");
                $arrayexpr
            ),
            true
        end
    end

    axisexpr = :(JuMP.Containers.DenseAxisArray{$T}(undef))
    append!(axisexpr.args, indexsets)

    if requestedtype == :DenseAxisArray
        has_dependent && return :(error(
            "Unable to create requested DenseAxisArray because index sets are dependent.",
        )),
        true
        return axisexpr, true
    end

    @assert requestedtype == :Auto
    if all(onetosets) # definitely Array
        return arrayexpr, true
    end

    # Fallback, we have to decide once we know the types of the symbolic sets
    condition = Expr(:&&)
    for (i, indexset) in enumerate(indexsets)
        if !onetosets[i]
            push!(condition.args, Expr(:call, :isa, indexset, :(Base.OneTo)))
        end
    end

    # TODO: check type stability
    return :(
        if $condition
            $arrayexpr
        else
            $axisexpr
        end
    ), true
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode