text:
# Unit Test for constraint type

using COSMO, Test, SparseArrays, Random

# This test is precision agnostic and will be run with TestFloat precision
if !@isdefined(UnitTestFloats)
    UnitTestFloats = [Float64] #if not run in full test setup, just do it for one float type
end

@testset "Constraints" begin

for TestFloat in UnitTestFloats

  rng = Random.MersenneTwister(1872381)

  A_int = 4
  b_int = 2

  A_Float = TestFloat(4.)
  b_Float = TestFloat(2.)

  A_UInt = UInt64(4)
  b_UInt = UInt64(2)

  A_sparse = sprand(rng, TestFloat, 10, 2, 0.4)
  b_sparse = sparse(rand(rng, TestFloat, 10))

  A_rvec = TestFloat[1 2 3 4]
  b_rvec = TestFloat(1)

  A_colvec = TestFloat[1; 2; 3; 4]
  b_colvec = TestFloat[4; 3; 2; 1]

  A_mat = rand(rng, TestFloat, 10, 10)
  b_mat = sparse(rand(rng, TestFloat, 10, 1))

  @testset "Constraints (T = $(TestFloat))" begin

    @testset "Constructors" begin
      @test typeof(COSMO.Constraint(A_int, b_int, COSMO.ZeroSet)) <: COSMO.Constraint{Float64} #integers are always converted to Float64
      @test typeof(COSMO.Constraint(A_Float, b_Float, COSMO.ZeroSet)) <: COSMO.Constraint{TestFloat}
      @test typeof(COSMO.Constraint(A_UInt, b_UInt, COSMO.ZeroSet)) <: COSMO.Constraint{Float64}
      @test typeof(COSMO.Constraint(A_sparse, b_sparse, COSMO.ZeroSet)) <: COSMO.Constraint{TestFloat}
      @test typeof(COSMO.Constraint(A_rvec, b_rvec, COSMO.ZeroSet)) <: COSMO.Constraint{TestFloat}
      @test typeof(COSMO.Constraint(A_colvec, b_colvec, COSMO.ZeroSet)) <: COSMO.Constraint{TestFloat}
      @test typeof(COSMO.Constraint(A_mat, b_mat, COSMO.ZeroSet)) <: COSMO.Constraint{TestFloat}
    end

    @testset "Indices" begin
      A = rand(rng, TestFloat, 3, 3)
      b = rand(rng, TestFloat, 3)
      ind = 3:5
      dim = 10
      cs = COSMO.Constraint(A, b, COSMO.ZeroSet, dim, ind)

      @test cs.A[:, ind] == A
     end

    @testset "Merge Constraints" begin
      A1 = rand(rng, TestFloat, 10, 10)
      b1 = rand(rng, TestFloat, 10)
      A2 = rand(rng, TestFloat, 10, 10)
      b2 = rand(rng, TestFloat, 10)
      Am = [A1; A2]
      bm = [b1; b2]
      c1 = COSMO.Constraint(A1, b1, COSMO.ZeroSet)
      c2 = COSMO.Constraint(A2, b2, COSMO.ZeroSet)
      cArr = [c1; c2]
      cm = COSMO.Constraint(Am, bm, COSMO.ZeroSet)
      COSMO.merge_constraints!(cArr)

      @test cArr[1].A == cm.A
      @test cArr[1].b == cm.b
      @test cArr[1].convex_set.dim == cm.convex_set.dim
      @test typeof(cArr[1].convex_set) == typeof(cm.convex_set)

      c1 = COSMO.Constraint(A1, b1, COSMO.Nonnegatives)
      c2 = COSMO.Constraint(A2, b2, COSMO.Nonnegatives)
      cArr = [c1; c2]
      cm = COSMO.Constraint(Am, bm, COSMO.Nonnegatives)
      COSMO.merge_constraints!(cArr)

      @test cArr[1].A == cm.A
      @test cArr[1].b == cm.b
      @test cArr[1].convex_set.dim == cm.convex_set.dim
      @test typeof(cArr[1].convex_set) == typeof(cm.convex_set)
    end
  end
end
end
nothing


text:
"""
    MattssonAlmquistCarpenter2014Optimal()

Coefficients of the optimal SBP operators with nonuniform grid given in
- Mattsson, Almquist, Carpenter (2014)
  Optimal diagonal-norm SBP operators.
  Journal of Computational Physics 264, pp. 91-111.
"""
struct MattssonAlmquistCarpenter2014Optimal <: SourceOfCoefficients end

function Base.show(io::IO, source::MattssonAlmquistCarpenter2014Optimal)
    if get(io, :compact, false)
        summary(io, source)
    else
        print(io,
            "Mattsson, Almquist, Carpenter (2014) \n",
            "  Optimal diagonal-norm SBP operators. \n",
            "  Journal of Computational Physics 264, pp. 91-111. \n")
    end
end

struct MattssonAlmquistCarpenter2014OptimalGrid{T,Grid} <: AbstractArray{T,1}
    xmin::T
    xmax::T
    d1::T
    d2::T
    d3::T
    uniform_grid::Grid
end

function MattssonAlmquistCarpenter2014OptimalGrid(xmin::T, xmax::T, d1::T, d2::T, d3::T, N::Int) where T
    @argcheck xmin < xmax
    @argcheck N > 6
    @argcheck d1 > 0
    @argcheck d2 > 0
    @argcheck d3 > 0

    uniform_grid = range(xmin+d1+d2+d3, stop=xmax-d1-d2-d3, length=N-6)
    MattssonAlmquistCarpenter2014OptimalGrid{T,typeof(uniform_grid)}(xmin, xmax, d1, d2, d3, uniform_grid)
end

function MattssonAlmquistCarpenter2014OptimalGrid(xmin, xmax, d1, d2, d3, N)
    MattssonAlmquistCarpenter2014OptimalGrid(promote(xmin, xmax, d1, d2, d3)..., N)
end

function Base.length(grid::MattssonAlmquistCarpenter2014OptimalGrid)
    length(grid.uniform_grid) + 6
end

function Base.getindex(grid::MattssonAlmquistCarpenter2014OptimalGrid, i::Int)
    N = length(grid)
    @boundscheck begin
        @argcheck i > 0
        @argcheck i <= N
    end

    if i == 1
        grid.xmin
    elseif i == 2
        grid.xmin + grid.d1
    elseif i == 3
        grid.xmin + grid.d1 + grid.d2
    elseif i == N
        grid.xmax
    elseif i == N-1
        grid.xmax - grid.d1
    elseif i == N-2
        grid.xmax - grid.d1 - grid.d2
    else
        @inbounds x = getindex(grid.uniform_grid, i-3)
        x
    end
end

Base.size(grid::MattssonAlmquistCarpenter2014OptimalGrid) = (length(grid),)
Base.step(grid::MattssonAlmquistCarpenter2014OptimalGrid) = step(grid.uniform_grid)

function construct_grid(::MattssonAlmquistCarpenter2014Optimal, accuracy_order, xmin, xmax, N)
    @argcheck N > 6
    T = promote_type(typeof(xmin), typeof(xmax))

    if accuracy_order == 2
        d1 = T(78866488858096586513//Int128(10)^20)
        d2 = T(95915098594220826013//Int128(10)^20)
        d3 = T(1)
    elseif accuracy_order == 4
        d1 = T(72181367003646814327//Int128(10)^20)
        d2 = T(13409118421582217252//Int128(10)^19)
        d3 = T(12898797485951900258//Int128(10)^19)
    elseif accuracy_order == 6
        d1 = T(51670081689316731234//Int128(10)^20)
        d2 = T(98190527037374634269//Int128(10)^20)
        d3 = T(10868393364992957832//Int128(10)^19)
    elseif accuracy_order == 8
        d1 = T(41669687672575697416//Int128(10)^20)
        d2 = T(78703773886730090312//Int128(10)^20)
        d3 = T(92685925671601406028//Int128(10)^20)
    else
        throw(ArgumentError("Accuracy order $accuracy_order not implemented/derived."))
    end

    h = (xmax-xmin) / (2*(d1+d2+d3) + N - 7)

    MattssonAlmquistCarpenter2014OptimalGrid(xmin, xmax, d1*h, d2*h, d3*h, N)
end

function first_derivative_coefficients(source::MattssonAlmquistCarpenter2014Optimal,
                                       order::Int, T=Float64, mode=FastMode())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,3}(SVector(T(-50000000000000000000//33743097329453577701),
                                                    T(55932483188770411252//33743097329453577701),
                                                    T(-5932483188770411252//33743097329453577701) )),
            # d2
            DerivativeCoefficientRow{T,1,3}(SVector(T(-13983120797192602813//24439920504708372824),
                                                    T(0),
                                                    T(13983120797192602813//24439920504708372824) )),
            # d3
            DerivativeCoefficientRow{T,1,4}(SVector(T(2966241594385205626//46639404052015171765),
                                                    T(-27966241594385205626//46639404052015171765),
                                                    T(0),
                                                    T(5000000000000000000//9327880810403034353) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(1//2))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(33743097329453577701//Int128(10)^20),
                                T(97759682018833491296//Int128(10)^20),
                                T(93278808104030343530//Int128(10)^20) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 4
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,5}(SVector(T(-50000000000000000000//21427296612044126417),
                                                    T(66884898686930380508//21427296612044126417),
                                                    T(-25171531878753856238//21427296612044126417),
                                                    T(10997619816825822803//21427296612044126417),
                                                    T(-13554933125011735296//107136483060220632085) )),
            # d2
            DerivativeCoefficientRow{T,1,5}(SVector(T(-16721224671732595127//28093989712218483700),
                                                    T(0),
                                                    T(92214436948640491071//112375958848873934800),
                                                    T(-16206184326771260201//56187979424436967400),
                                                    T(17707075979581024571//280939897122184837000) )),
            # d3
            DerivativeCoefficientRow{T,1,5}(SVector(T(12585765939376928119//71722939624706300000),
                                                    T(-92214436948640491071//143445879249412600000),
                                                    T(0),
                                                    T(1636075617843355867//2868917584988252000),
                                                    T(-14760875822281158529//143445879249412600000) )),
            # d4
            DerivativeCoefficientRow{T,1,6}(SVector(T(-10997619816825822803//109173230217361308360),
                                                    T(16206184326771260201//54586615108680654180),
                                                    T(-8180378089216779335//10917323021736130836),
                                                    T(0),
                                                    T(17180591347196107273//27293307554340327090),
                                                    T(-625000000000000000//8187992266302098127) )),
            # d5
            DerivativeCoefficientRow{T,1,7}(SVector(T(2259155520835289216//82365134292746683125),
                                                    T(-17707075979581024571//247095402878240049375),
                                                    T(14760875822281158529//98838161151296019750),
                                                    T(-34361182694392214546//49419080575648009875),
                                                    T(0),
                                                    T(800000000000000000//1186057933815552237),
                                                    T(-100000000000000000//1186057933815552237) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(2//3), T(-1//12))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(21427296612044126417//100000000000000000000),
                                T(280939897122184837//250000000000000000),
                                T(717229396247063//500000000000000),
                                T(2729330755434032709//2500000000000000000),
                                T(395352644605184079//400000000000000000) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,7}(SVector(T(-3125000000000000000//944357158252257333),
                                                    T(22223596967296279011//5036571510678705776),
                                                    T(-5854697895099786871//3777428633009029332),
                                                    T(28083754862953532289//50365715106787057760),
                                                    T(-1592329083817967435//15109714532036117328),
                                                    T(-15653772860347171721//1510971453203611732800),
                                                    T(23913677348563239189//5036571510678705776000) )),
            # d2
            DerivativeCoefficientRow{T,1,7}(SVector(T(-22223596967296279011//26989195119035671001),
                                                    T(0),
                                                    T(89405599296515541581//80967585357107013003),
                                                    T(-28597427787314667440//80967585357107013003),
                                                    T(57056178538117177397//809675853571070130030),
                                                    T(41320613074890940489//8096758535710701300300),
                                                    T(-5124091837453295329//1619351707142140260060) )),
            # d3
            DerivativeCoefficientRow{T,1,7}(SVector(T(5854697895099786871//27278567870198137125),
                                                    T(-89405599296515541581//109114271480792548500),
                                                    T(0),
                                                    T(27653905086569037761//36371423826930849500),
                                                    T(-6077915680998146409//36371423826930849500),
                                                    T(83784382166533084621//10911427148079254850000),
                                                    T(42099567773838744673//10911427148079254850000) )),
            # d4
            DerivativeCoefficientRow{T,1,7}(SVector(T(-84251264588860596867//1043526904157157775600),
                                                    T(714935694682866686//2608817260392894439),
                                                    T(-82961715259707113283//104352690415715777560),
                                                    T(0),
                                                    T(75419218459746682761//104352690415715777560),
                                                    T(-14034899831339963049//104352690415715777560),
                                                    T(3512738257179466361//260881726039289443900) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(1592329083817967435//98680905919946100728),
                                                    T(-57056178538117177397//986809059199461007280),
                                                    T(18233747042994439227//98680905919946100728),
                                                    T(-75419218459746682761//98680905919946100728),
                                                    T(0),
                                                    T(18687868497479752801//24670226479986525182),
                                                    T(-15119380469839682133//98680905919946100728),
                                                    T(625000000000000000//37005339719979787773) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(15653772860347171721//10037581831426163456000),
                                                    T(-41320613074890940489//10037581831426163456000),
                                                    T(-29162680879405877//3493763254934272000),
                                                    T(14034899831339963049//100375818314261634560),
                                                    T(-18687868497479752801//25093954578565408640),
                                                    T(0),
                                                    T(4696526232232696867//6273488644641352160),
                                                    T(-5859375000000000//39209304029008451),
                                                    T(1953125000000000//117627912087025353) )),
            # d7
            DerivativeCoefficientRow{T,1,10}(SVector(T(-71741032045689717567//99943556356761752125000),
                                                    T(5124091837453295329//1998871127135235042500),
                                                    T(-42099567773838744673//9994355635676175212500),
                                                    T(-7025476514358932722//499717781783808760625),
                                                    T(15119380469839682133//99943556356761752125),
                                                    T(-75144419715723149872//99943556356761752125),
                                                    T(0),
                                                    T(600000000000000000//799548450854094017),
                                                    T(-120000000000000000//799548450854094017),
                                                    T(40000000000000000//2398645352562282051) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(3//4), T(-3//20), T(1//60))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(944357158252257333//6250000000000000000),
                                T(80967585357107013003//100000000000000000000),
                                T(218228542961585097//200000000000000000),
                                T(2608817260392894439//2500000000000000000),
                                T(12335113239993262591//12500000000000000000),
                                T(39209304029008451//39062500000000000),
                                T(799548450854094017//800000000000000000) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    #= FIXME: There seem to be some errors in the coefficients...
    elseif order == 8
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,8}(SVector(T(-25000000000000000000//6081611055353751439),
                                                    T(66670790901888837033//12163222110707502878),
                                                    T(-10997015095317519523//6081611055353751439),
                                                    T(61752567584332553851//121632221107075028780),
                                                    T(-32312350944133128873//1216322211070750287800),
                                                    T(-16967490160001675093//608161105535375143900),
                                                    T(3031405594112644741//24326444221415005756000),
                                                    T(32350133072942893419//12163222110707502878000) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(-1257939450979034661//1230864766727295094),
                                                    T(0),
                                                    T(86688767821045233147//65235832636546639982),
                                                    T(-24298087640343350527//65235832636546639982),
                                                    T(39549469619698650847//652358326365466399820),
                                                    T(2076352837148473751//652358326365466399820),
                                                    T(-4065343087310119561//4077239539784164998875),
                                                    T(-8167205603296830489//13047166527309327996400) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(10997015095317519523//43865207099050505477),
                                                    T(-86688767821045233147//87730414198101010954),
                                                    T(0),
                                                    T(41032546292236417573//43865207099050505477),
                                                    T(-21014872891771196683//87730414198101010954),
                                                    T(38572177503610408523//877304141981010109540),
                                                    T(-16637807199883547459//87730414198101010954000),
                                                    T(-11580435163412724923//219326035495252527385000) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(-61752567584332553851//973889517710795427990),
                                                    T(24298087640343350527//97388951771079542799),
                                                    T(-82065092584472835146//97388951771079542799),
                                                    T(0),
                                                    T(81102837324727866266//97388951771079542799),
                                                    T(-20423896795865859484//97388951771079542799),
                                                    T(776602705809652689//21641989282462120622),
                                                    T(1285505569126131521//541049732061553015550) )),
            # d5
            DerivativeCoefficientRow{T,1,9}(SVector(T(32312350944133128873//10072514376844677230000),
                                                    T(-39549469619698650847//1007251437684467723000),
                                                    T(21014872891771196683//100725143768446772300),
                                                    T(-40551418662363933133//50362571884223386150),
                                                    T(0),
                                                    T(3202162603743798001//4029005750737870892),
                                                    T(-1969916799269047253//10072514376844677230),
                                                    T(37220336417235830241//1007251437684467723000),
                                                    T(-25000000000000000//7050760063791274061) )),
            # d6
            DerivativeCoefficientRow{T,1,10}(SVector(T(16967490160001675093//4988436332888823941700),
                                                    T(-296621833878353393//142526752368252112620),
                                                    T(-38572177503610408523//997687266577764788340),
                                                    T(10211948397932929742//49884363328888239417),
                                                    T(-80054065093594950025//99768726657776478834),
                                                    T(0),
                                                    T(39951539656523293130//49884363328888239417),
                                                    T(-2222198748535843633//11085414073086275426),
                                                    T(40000000000000000000//1047571629906653027757),
                                                    T(-1250000000000000000//349190543302217675919) )),
            # d7
            DerivativeCoefficientRow{T,1,11}(SVector(T(-3031405594112644741//200106059975821710280000),
                                                    T(4065343087310119561//6253314374244428446250),
                                                    T(16637807199883547459//100053029987910855140000),
                                                    T(-6989424352286874201//200106059975821710280),
                                                    T(1969916799269047253//10005302998791085514),
                                                    T(-3995153965652329313//5002651499395542757),
                                                    T(0),
                                                    T(40008167342759928887//50026514993955427570),
                                                    T(-1000000000000000000//5002651499395542757),
                                                    T(4000000000000000000//105055681487306397897),
                                                    T(-125000000000000000//35018560495768799299) )),
            # d8
            DerivativeCoefficientRow{T,1,12}(SVector(T(-32350133072942893419//99994066100338390832000),
                                                    T(8167205603296830489//19998813220067678166400),
                                                    T(11580435163412724923//249985165250845977080000),
                                                    T(-11569550122135183689//4999703305016919541600),
                                                    T(-37220336417235830241//999940661003383908320),
                                                    T(19999788736822592697//99994066100338390832),
                                                    T(-40008167342759928887//49997033050169195416),
                                                    T(0),
                                                    T(5000000000000000000//6249629131271149427),
                                                    T(-1250000000000000000//6249629131271149427),
                                                    T(5000000000000000000//131242211756694137967),
                                                    T(-156250000000000000//43747403918898045989) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(6081611055353751439//50000000000000000000),
                                T(32617916318273319991//50000000000000000000),
                                T(43865207099050505477//50000000000000000000),
                                T(97388951771079542799//100000000000000000000),
                                T(1007251437684467723//1000000000000000000),
                                T(49884363328888239417//50000000000000000000),
                                T(5002651499395542757//5000000000000000000),
                                T(6249629131271149427//6250000000000000000) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    =#
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

ext: jl


text:
using Test, MLToolkit
using Distributions: Beta

dims = (10, 5)

@testset "Activations" begin
    @testset "break_stick_ibp" begin
        ν = ones(2, 3) ./ 2
        π = break_stick_ibp(ν)
        @test π ≈ [0.5 0.5 0.5; 0.25 0.25 0.25]
    end
    @testset "break_logstick_ibp" begin
        for _ = 1:10
            ν = rand(Beta(2, 1), dims...)
            logν = log.(ν)
            @test all(log.(break_stick_ibp(ν)) .≈ break_logstick_ibp(logν))
        end
    end
end


text:
using Keras
using Flux
using PyCall
@pyimpory numpy

X = numpy.load("X.npy")
Y = numpy.load("Y.npy")

model = Keras.load("structure.json", "weights.h5")

max_stars_repo_path: examples/sentiment analysis using glove embeddings/main.jl


text:
"""

    ransac(setting; t, w=0.5, m=0, k=0, d=0, confidence=0.99)

Run the RANSAC (1981) algorithm for the given regression setting

# Arguments
 - `setting::RegressionSetting`: RegressionSetting object with a formula and a dataset.
 - `t::Float64`: The threshold distance of a sample point to the regression hyperplane to determine if it fits the model well.
 - `w::Float64`: The probability of a sample point being inlier, default=0.5.
 - `m::Int`: The number of points to sample to estimate the model parameter for each iteration. If set to 0, defaults to picking p points which is the minimum required.
 - `k::Int`: The number of iterations to run. If set to 0, is calculated according to the formula given in the paper based on outlier probability and the sample set size.
 - `d::Int`: The number of close data points required to accept the model. Defaults to number of data points multiplied by inlier ratio.
 - `confidence::Float64`: Required to determine the number of optimum iterations if k is not specified.

 # Output
- `["outliers"]`: Array of indices of outliers.

 # Examples
```julia-repl
julia> df = DataFrame(y=[0,1,2,3,3,4,10], x=[0,1,2,2,3,4,2])
julia> reg = createRegressionSetting(@formula(y ~ x), df)
julia> ransac(reg, t=0.8, w=0.85)
Dict{String,Array{Int64,1}} with 1 entry:
  "outliers" => [7]
```

# References
Martin A. Fischler & Robert C. Bolles (June 1981). "Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography"
Comm. ACM. 24 (6): 381–395.
"""
function ransac(setting::RegressionSetting; t::Float64, w::Float64=0.5, m::Int=0, k::Int=0, d::Int=0, confidence::Float64=0.99)
    X, y = @extractRegressionSetting setting
    return ransac(X, y, t=t, w=w, m=m, k=k, d=d, confidence=confidence)
end

function ransac(X::Array{Float64,2}, y::Array{Float64,1}; t::Float64, w::Float64=0.5, m::Int=0, k::Int=0, d::Int=0, confidence::Float64=0.99)

    n, p = size(X)
    if d == 0
        d = Int(floor(n * w))
    end

    if k == 0
        k = Int(ceil(log(1 - confidence) / log(1 - w^d)))
    end

    if m == 0
        m = p
    end

    iteration_inlier_indices = zeros(Int, n)

    maximum_count = d - 1
    maximum_inlier_indices = zeros(Int, n)
    minimum_error = Inf

    for iteration in 1:k
        inliers_count = 0
        sampled_indices = sample(1:n, m, replace=false)
        ols_sampled_points = ols(X[sampled_indices, :], y[sampled_indices])
        betas = coef(ols_sampled_points)

        e = abs.(y - X * betas) ./ norm([1; betas[2:end]], 2)

        iteration_inlier_indices = filter(i -> e[i] < t, 1:n)
        inliers_count = length(iteration_inlier_indices)

        if inliers_count >= d
            ols_inliers = ols(X[iteration_inlier_indices, :], y[iteration_inlier_indices])
            error_vector = residuals(ols_inliers)
            iteration_error = norm(error_vector, 2) / length(error_vector)

            if iteration_error < minimum_error
                maximum_inlier_indices = iteration_inlier_indices
                maximum_count = inliers_count
                minimum_error = iteration_error
            end
        end
    end

    result = Dict(
        "outliers" => setdiff(1:n, maximum_inlier_indices)
    )
    return result
end

max_stars_repo_name: jbytecode/LinRegOutliers

max_issues_repo_name: jbytecode/LinRegOutliers

max_forks_repo_name: jbytecode/LinRegOutliers


text:
"""
    TextConfig

    TextConfig(; kw...)
    TextConfig(face, namepixels, namepos, timepixels, timepos, fcolor, bcolor)

Text configuration for printing timestep and grid name on the image.

# Arguments / Keywords

- `font`: `String` font name.
- `namepixels` and `timepixels`: the pixel size of the font.
- `timepos` and `namepos`: tuples that set the label positions, in `Int` pixels.
- `fcolor` and `bcolor`: the foreground and background colors, as `ARGB32`.
"""
struct TextConfig{F,NPi,NPo,TPi,TPo,FC,BC}
    face::F
    namepixels::NPi
    namepos::NPo
    timepixels::TPi
    timepos::TPo
    fcolor::FC
    bcolor::BC
end
function TextConfig(;
    font=autofont(), namepixels=12, timepixels=12,
    namepos=(3timepixels + namepixels, timepixels),
    timepos=(2timepixels, timepixels),
    fcolor=ARGB32(1.0), bcolor=ZEROCOL,
)
    if font isa AbstractString
        face = FreeTypeAbstraction.findfont(font)
        face isa Nothing && _fontnotfounderror(font)
    else
        _fontnotstring(font)
    end
    TextConfig(face, namepixels, namepos, timepixels, timepos, fcolor, bcolor)
end

function autofont()
    fonts = if Sys.islinux()
        ("cantarell", "sans-serif", "Bookman")
    else
        ("arial", "sans-serif")
    end
    for font in fonts
        face = FreeTypeAbstraction.findfont(font)
        face isa Nothing || return font
    end
    _nodefaultfonterror(fonts)
end

@noinline _fontnotstring(font) = throw(ArgumentError("font $font is not a String"))

@noinline _fontnotfounderror(font) =
    throw(ArgumentError(
        """
        Font "$font" wasn't be found in this system. Specify an existing font name
        with the `font` keyword, or use `text=nothing` to display no text."
        """
    ))
@noinline _nodefaultfonterror(font) =
    error(
        """
        Your system does not contain the default font $font. Specify an existing font
        name `String` with the keyword-argument `font`, for the `Output` or `ImageConfig`.
        """
    )

# Render time `name` and `t` as text onto the image, following config settings.
function _rendertime! end

function _rendertext!(img, config::TextConfig, name, t)
    _rendername!(img, config::TextConfig, name)
    _rendertime!(img, config::TextConfig, t)
    img
end
_rendertext!(img, config::Nothing, name, t) = nothing

# Render `name` as text on the image following config settings.
function _rendername!(img, config::TextConfig, name)
    renderstring!(img, name, config.face, config.namepixels, config.namepos...;
        fcolor=config.fcolor, bcolor=config.bcolor
    )
    img
end
_rendername!(img, config::TextConfig, name::Nothing) = img
_rendername!(img, config::Nothing, name) = img
_rendername!(img, config::Nothing, name::Nothing) = img

# Render time `t` as text on the image following config settings.
function _rendertime!(img, config::TextConfig, t)
    renderstring!(img, string(t), config.face, config.timepixels, config.timepos...;
        fcolor=config.fcolor, bcolor=config.bcolor
    )
    img
end
_rendertime!(img, config::Nothing, t) = img
_rendertime!(img, config::TextConfig, t::Nothing) = img
_rendertime!(img, config::Nothing, t::Nothing) = img


text:
using Test
using ModelingToolkit, Graphs, DiffEqJump

import ModelingToolkit: value

#################################
#  testing for Jumps / all dgs
#################################
@parameters k1 k2
@variables t S(t) I(t) R(t)
j₁ = MassActionJump(k1, [0 => 1], [S => 1])
j₂ = MassActionJump(k1, [S => 1], [S => -1])
j₃ = MassActionJump(k2, [S =>1, I => 1], [S => -1, I => 1])
j₄ = MassActionJump(k2, [S => 2, R => 1], [R => -1])
j₅ = ConstantRateJump(k1*I, [R ~ R + 1])
j₆ = VariableRateJump(k1*k2/(1+t)*S, [S ~ S - 1, R ~ R + 1])
eqs = [j₁,j₂,j₃,j₄,j₅,j₆]
@named js = JumpSystem(eqs, t, [S,I,R],[k1,k2])
S  = value(S); I = value(I); R = value(R)
k1 = value(k1); k2 = value(k2)
# eq to vars they depend on
eq_sdeps  = [Variable[], [S], [S,I], [S,R], [I], [S]]
eq_sidepsf = [Int[], [1], [1,2], [1,3], [2], [1]]
eq_sidepsb = [[2,3,4,6], [3,5],[4]]
deps = equation_dependencies(js)
@test all(i -> isequal(Set(eq_sdeps[i]),Set(deps[i])), 1:length(eqs))
depsbg = asgraph(js)
@test depsbg.fadjlist == eq_sidepsf
@test depsbg.badjlist == eq_sidepsb

# eq to params they depend on
eq_pdeps    = [[k1],[k1],[k2],[k2],[k1],[k1,k2]]
eq_pidepsf = [[1],[1],[2],[2],[1],[1,2]]
eq_pidepsb = [[1,2,5,6],[3,4,6]]
deps = equation_dependencies(js, variables=parameters(js))
@test all(i -> isequal(Set(eq_pdeps[i]),Set(deps[i])), 1:length(eqs))
depsbg2 = asgraph(js, variables=parameters(js))
@test depsbg2.fadjlist == eq_pidepsf
@test depsbg2.badjlist == eq_pidepsb

# var to eqs that modify them
s_eqdepsf = [[1,2,3,6],[3],[4,5,6]]
s_eqdepsb = [[1],[1],[1,2],[3],[3],[1,3]]
ne        = 8
bg        = BipartiteGraph(ne, s_eqdepsf, s_eqdepsb)
deps2     = variable_dependencies(js)
@test isequal(bg,deps2)

# eq to eqs that depend on them
eq_eqdeps = [[2,3,4,6],[2,3,4,6],[2,3,4,5,6],[4],[4],[2,3,4,6]]
dg = SimpleDiGraph(6)
for (eqidx,eqdeps) in enumerate(eq_eqdeps)
    for eqdepidx in eqdeps
       add_edge!(dg, eqidx, eqdepidx)
    end
end
dg3 = eqeq_dependencies(depsbg,deps2)
@test dg == dg3

# var to vars that depend on them
var_vardeps = [[1,2,3],[1,2,3],[3]]
ne = 7
dg = SimpleDiGraph(3)
for (vidx,vdeps) in enumerate(var_vardeps)
    for vdepidx in vdeps
       add_edge!(dg, vidx, vdepidx)
    end
end
dg4 = varvar_dependencies(depsbg,deps2)
@test dg == dg4

#####################################
#       testing for ODE/SDEs
#####################################
@parameters k1 k2
@variables t S(t) I(t) R(t)
D = Differential(t)
eqs = [D(S) ~ k1 - k1*S - k2*S*I - k1*k2/(1+t)*S,
       D(I) ~ k2*S*I,
       D(R) ~ -k2*S^2*R/2 + k1*I + k1*k2*S/(1+t)]
noiseeqs = [S,I,R]
@named os = ODESystem(eqs, t, [S,I,R], [k1,k2])
deps = equation_dependencies(os)
S  = value(S); I = value(I); R = value(R)
k1 = value(k1); k2 = value(k2)
eq_sdeps  = [[S,I], [S,I], [S,I,R]]
@test all(i -> isequal(Set(eq_sdeps[i]),Set(deps[i])), 1:length(deps))

@parameters k1 k2
@variables t S(t) I(t) R(t)
@named sdes = SDESystem(eqs, noiseeqs, t, [S,I,R], [k1,k2])
deps = equation_dependencies(sdes)
@test all(i -> isequal(Set(eq_sdeps[i]),Set(deps[i])), 1:length(deps))

deps = variable_dependencies(os)
s_eqdeps = [[1],[2],[3]]
@test deps.fadjlist == s_eqdeps

#####################################
#       testing for nonlin sys
#####################################
@variables x y z
@parameters σ ρ β

eqs = [0 ~ σ*(y-x),
       0 ~ ρ-y,
       0 ~ y - β*z]
@named ns = NonlinearSystem(eqs, [x,y,z],[σ,ρ,β])
deps = equation_dependencies(ns)
eq_sdeps = [[x,y],[y],[y,z]]
@test all(i -> isequal(Set(deps[i]),Set(value.(eq_sdeps[i]))), 1:length(deps))

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
"""
    struct Density{M,B}
        μ::M
        base::B
    end

For measures μ and ν with μ≪ν, the density of μ with respect to ν (also called
the Radon-Nikodym derivative dμ/dν) is a function f defined on the support of ν
with the property that for any measurable a ⊂ supp(ν), μ(a) = ∫ₐ f dν.

Because this function is often difficult to express in closed form, there are
many different ways of computing it. We therefore provide a formal
representation to allow comptuational flexibilty.
"""
struct Density{M,B,L}
    μ::M
    base::B
    log::L
end

export 𝒹

export log𝒹

log𝒹(μ, base) = Density(μ, base, Val{true}())

"""
    𝒹(μ::AbstractMeasure, base::AbstractMeasure; log=false)

Compute the Radom-Nikodym derivative (or its log, if `log=false`) of μ with
respect to `base`.
"""
function 𝒹(μ::AbstractMeasure, base::AbstractMeasure; log = false)
    return Density(μ, base, Val(log))
end

(f::Density{M,B,Val{true}})(x) where {M,B} = logdensity(f.μ, f.base, x)

(f::Density{M,B,Val{false}})(x) where {M,B} = density(f.μ, f.base, x)

"""
    struct DensityMeasure{F,B} <: AbstractMeasure
        density :: F
        base    :: B
    end

A `DensityMeasure` is a measure defined by a density with respect to some other
"base" measure
"""
struct DensityMeasure{F,B,L} <: AbstractMeasure
    f::F
    base::B
    log::L
end

function Pretty.tile(μ::DensityMeasure{F,B,Val{L}}) where {F,B,L}
    result = Pretty.literal("DensityMeasure ∫(")
    result *= Pretty.pair_layout(Pretty.tile(μ.f), Pretty.tile(μ.base); sep = ", ")
    result *= Pretty.literal("; log = ")
    result *= Pretty.tile(L)
    result *= Pretty.literal(")")
end

function Base.rand(rng::AbstractRNG, T::Type, d::DensityMeasure)
    x = rand(rng, T, d.base)
    WeightedMeasure(d.f(x), Dirac(x))
end

basemeasure(μ::DensityMeasure) = μ.base

logdensity(μ::DensityMeasure{F,B,Val{true}}, x) where {F,B} = μ.f(x)

density(μ::DensityMeasure{F,B,Val{false}}, x) where {F,B} = μ.f(x)

logdensity(μ::DensityMeasure{F,B,Val{false}}, x) where {F,B} = log(density(μ, x))

export ∫

"""
    ∫(f, base::AbstractMeasure)

Define a new measure in terms of a density `f` over some measure `base`.
"""
∫(f, base::AbstractMeasure) = DensityMeasure(f, base, Val(false))

∫(μ::AbstractMeasure, base::AbstractMeasure) = ∫exp(log𝒹(μ, base), base)

export ∫exp

"""
    ∫exp(f, base::AbstractMeasure; log=false)

Define a new measure in terms of a density `f` over some measure `base`.
"""
∫exp(f, μ) = DensityMeasure(f, μ, Val{true}())

# TODO: `density` and `logdensity` functions for `DensityMeasure`

@inline function logdensity(μ::T, ν::T, x) where {T<:AbstractMeasure}
    μ == ν && return 0.0
    invoke(logdensity, Tuple{AbstractMeasure,AbstractMeasure,typeof(x)}, μ, ν, x)
end

@inline function logdensity(μ::AbstractMeasure, ν::AbstractMeasure, x)
    α = basemeasure(μ)
    β = basemeasure(ν)

    # If α===μ and β===ν, The recursive call would be exactly the same as the
    # original one. We need to break the recursion.
    if α === μ && β === ν
        @warn """
        No method found for logdensity(μ, ν, x) where
        typeof(μ) == $(typeof(μ))
        typeof(ν) == $(typeof(ν))

        Returning NaN. If this is incorrect, please add a method
        logdensity(μ::$(typeof(μ)), ν::$(typeof(ν)), x)
        """
        return NaN
    end

    # Infinite or NaN results occur when outside the support of α or β,
    # and also when one measure is singular wrt the other. Computing the base
    # measures first is often much cheaper, and allows the numerically-intensive
    # computation to "fall through" in these cases.
    # TODO: Add tests to check that NaN cases work properly
    ℓ = logdensity(α, β, x)
    isfinite(ℓ) || return ℓ

    ℓ += logdensity(μ, x)
    ℓ -= logdensity(ν, x)

    return ℓ
end

function logpdf(d::AbstractMeasure, x)
    _logpdf(d, basemeasure(d), x)
end

@inline function _logpdf(d::AbstractMeasure, β::AbstractMeasure, x, ℓ = zero(Float64))
    # @show d
    # @show x
    Δℓ = logdensity(d, x)
    newℓ = ℓ + Δℓ
    d == β && return newℓ
    # @show Δℓ, d
    _logpdf(β, basemeasure(β), x, newℓ)
end

logdensity(::Lebesgue, ::Lebesgue, x) = 0.0

# logdensity(::Lebesgue{ℝ}, ::Lebesgue{ℝ}, x) = zero(x)

export density

density(μ, ν::AbstractMeasure, x) = exp(logdensity(μ, ν, x))

density(μ, x) = exp(logdensity(μ, x))


text:
## Kirchhoff-Love plate
function (op::Operator{:∫κᵢⱼMᵢⱼdΩ})(ap::T,k::AbstractMatrix{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    D = op.D
    ν = op.ν
    for ξ in 𝓖
        _,_,_,B₁₁,B₁₂,B₂₂ = get∇²𝝭(ap,ξ)
        𝑤 = get𝑤(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                k[I,J] += D*(B₁₁[i]*B₁₁[j] + ν*(B₁₁[i]*B₂₂[j] + B₂₂[i]*B₁₁[j]) + B₂₂[i]*B₂₂[j] + 2*(1-ν)*B₁₂[i]*B₁₂[j])*𝑤
            end
        end
    end
end

function (op::Operator{:∫wqdΩ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            f[I] += N[i]*ξ.q*𝑤
        end
    end
end

function (op::Operator{:∫wVdΓ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N = get𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            f[I] += N[i]*ξ.V*𝑤
        end
    end
end

function (op::Operator{:∫θₙMₙₙdΓ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁,n₂ = get𝒏(ap)
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        _,B₁,B₂ = get∇𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            f[I] -= (B₁[i]*n₁+B₂[i]*n₂)*ξ.M*𝑤
        end
    end
end

function (op::Operator{:∫∇wMdΓ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁,n₂ = get𝒏(ap)
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        _,B₁,B₂ = get∇𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            f[I] -= (B₁[i]*ξ.M₁+B₂[i]*ξ.M₂)*𝑤
        end
    end
end

function (op::Operator{:∫∇𝑛vθdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁,n₂ = get𝒏(ap)
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        _,B₁,B₂ = get∇𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            θᵢ = B₁[i]*n₁+B₂[i]*n₂
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                θⱼ = B₁[j]*n₁+B₂[j]*n₂
                k[I,J] += op.α*θᵢ*θⱼ*𝑤
            end
            f[I] += op.α*θᵢ*ξ.θ*𝑤
        end
    end
end

function (op::Operator{:∫MₙₙθdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁,n₂,s₁,s₂ = get𝒏(ap)
    α = op.α
    D = op.D
    ν = op.ν
    D₁₁ = D*(n₁^2+ν*n₂^2)
    D₁₂ = 2*D*n₁*n₂*(1-ν)
    D₂₂ = D*(ν*n₁^2+n₂^2)
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        _,B₁,B₂,B₁₁,B₁₂,B₂₂ = get∇²𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            θᵢ = B₁[i]*n₁ + B₂[i]*n₂
            Mᵢ = D₁₁*B₁₁[i] + D₁₂*B₁₂[i] + D₂₂*B₂₂[i]
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                θⱼ = B₁[j]*n₁ + B₂[j]*n₂
                Mⱼ = D₁₁*B₁₁[j] + D₁₂*B₁₂[j] + D₂₂*B₂₂[j]
                k[I,J] += (-Mᵢ*θⱼ-θᵢ*Mⱼ+α*θᵢ*θⱼ)*𝑤
            end
            f[I] += (-Mᵢ+α*θᵢ)*ξ.θ*𝑤
        end
    end
end

function (op::Operator{:∫VgdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁,n₂,s₁,s₂ = get𝒏(ap)
    α = op.α
    D = op.D
    ν = op.ν
    D₁₁₁ = -D*(n₁ + n₁*s₁*s₁ + ν*n₂*s₁*s₂)
    D₁₁₂ = -D*(n₂ + n₂*s₁*s₁ + 2*n₁*s₁*s₂ + (n₂*s₂*s₂ - n₂*s₁*s₁ - n₁*s₁*s₂)*ν)
    D₁₂₂ = -D*(n₁ + n₁*s₂*s₂ + 2*n₂*s₁*s₂ + (n₁*s₁*s₁ - n₁*s₂*s₂ - n₂*s₁*s₂)*ν)
    D₂₂₂ = -D*(n₂ + n₂*s₂*s₂ + ν*n₁*s₁*s₂)
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        N,_,_,_,_,_,B₁₁₁,B₁₁₂,B₁₂₂,B₂₂₂ = get∇³𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            Vᵢ = D₁₁₁*B₁₁₁[i] + D₁₁₂*B₁₁₂[i] + D₁₂₂*B₁₂₂[i] + D₂₂₂*B₂₂₂[i]
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                Vⱼ = D₁₁₁*B₁₁₁[j] + D₁₁₂*B₁₁₂[j] + D₁₂₂*B₁₂₂[j] + D₂₂₂*B₂₂₂[j]
                k[I,J] += (-Vᵢ*N[j]-N[i]*Vⱼ+α*N[i]*N[j])*𝑤
            end
            f[I] += (-Vᵢ+α*N[i])*ξ.g*𝑤
        end
    end
end

function (op::Operator{:∫M̃ₙₙθdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁ = 𝓖[1].n₁
    n₂ = 𝓖[1].n₂
    s₁ = 𝓖[1].s₁
    s₂ = 𝓖[1].s₂
    # α = op.α
    D = op.D
    ν = op.ν
    D₁₁ = D*(n₁^2+ν*n₂^2)
    D₁₂ = 2*D*n₁*n₂*(1-ν)
    D₂₂ = D*(ν*n₁^2+n₂^2)
    for ξ in 𝓖
        𝑤 = get𝑤(ap,ξ)
        _,B₁,B₂,B₁₁,B₁₂,B₂₂ = get∇²𝝭(ap,ξ)
        B̄₁₁,B̄₁₂,B̄₂₂ = get∇̄²𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            θᵢ = B₁[i]*n₁ + B₂[i]*n₂
            Mᵢ = D₁₁*B₁₁[i] + D₁₂*B₁₂[i] + D₂₂*B₂₂[i]
            M̄ᵢ = D₁₁*B̄₁₁[i] + D₁₂*B̄₁₂[i] + D₂₂*B̄₂₂[i]
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                θⱼ = B₁[j]*n₁ + B₂[j]*n₂
                Mⱼ = D₁₁*B₁₁[j] + D₁₂*B₁₂[j] + D₂₂*B₂₂[j]
                # k[I,J] += (-Mᵢ*θⱼ-θᵢ*Mⱼ+α*θᵢ*θⱼ)*ξ.𝑤
                # k[I,J] += (-Mᵢ*θⱼ-θᵢ*Mⱼ+M̄ᵢ*θⱼ)*ξ.𝑤
                k[I,J] += M̄ᵢ*θⱼ*ξ.𝑤
            end
            # f[I] += (-Mᵢ+α*θᵢ)*ξ.θ*ξ.𝑤
            f[I] += (-Mᵢ+M̄ᵢ)*ξ.θ*ξ.𝑤
        end
    end
end

function (op::Operator{:∫ṼgdΓ})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; 𝓖 = ap.𝓖
    n₁ = 𝓖[1].n₁
    n₂ = 𝓖[1].n₂
    s₁ = 𝓖[1].s₁
    s₂ = 𝓖[1].s₂
    # α = op.α
    D = op.D
    ν = op.ν
    D₁₁₁ = -D*(n₁ + n₁*s₁*s₁ + ν*n₂*s₁*s₂)
    D₁₁₂ = -D*(n₁*s₁*s₂ + (n₂*s₂*s₂ + n₂)*ν)
    D₁₂₁ = -D*(n₂ + n₂*s₁*s₁ + n₁*s₁*s₂ + (-n₂ - n₂*s₁*s₁ - n₁*s₁*s₂)*ν)
    D₁₂₂ = -D*(n₁ + n₁*s₂*s₂ + n₂*s₁*s₂ + (-n₁ - n₁*s₂*s₂ - n₂*s₁*s₂)*ν)
    D₂₂₁ = -D*(n₂*s₁*s₂ + (n₁*s₁*s₁ + n₁)*ν)
    D₂₂₂ = -D*(n₂ + n₂*s₂*s₂ + ν*n₁*s₁*s₂)
    for ξ in 𝓖
        N,B₁₁₁,B₁₁₂,B₁₂₁,B₁₂₂,B₂₂₁,B₂₂₂ = get∇∇̃²𝝭(ap,ξ)
        B̄₁₁₁,B̄₁₁₂,B̄₁₂₁,B̄₁₂₂,B̄₂₂₁,B̄₂₂₂ = get∇∇̄²𝝭(ap,ξ)
        for (i,xᵢ) in enumerate(𝓒)
            I = xᵢ.id
            Vᵢ = D₁₁₁*B₁₁₁[i] + D₁₁₂*B₁₁₂[i] + D₁₂₁*B₁₂₁[i] + D₁₂₂*B₁₂₂[i] + D₂₂₁*B₂₂₁[i] + D₂₂₂*B₂₂₂[i]
            V̄ᵢ = D₁₁₁*B̄₁₁₁[i] + D₁₁₂*B̄₁₁₂[i] + D₁₂₁*B̄₁₂₁[i] + D₁₂₂*B̄₁₂₂[i] + D₂₂₁*B̄₂₂₁[i] + D₂₂₂*B̄₂₂₂[i]
            for (j,xⱼ) in enumerate(𝓒)
                J = xⱼ.id
                Vⱼ = D₁₁₁*B₁₁₁[j] + D₁₁₂*B₁₁₂[j] + D₁₂₁*B₁₂₁[j] + D₁₂₂*B₁₂₂[j] + D₂₂₁*B₂₂₁[j] + D₂₂₂*B₂₂₂[j]
                # k[I,J] += (-Vᵢ*N[j]-N[i]*Vⱼ+α*N[i]*N[j])*ξ.𝑤
                k[I,J] += (-Vᵢ*N[j]-N[i]*Vⱼ+V̄ᵢ*N[j])*ξ.𝑤
            end
            # f[I] += (-Vᵢ+α*N[i])*ξ.g*ξ.𝑤
            f[I] += (-Vᵢ+V̄ᵢ)*ξ.g*ξ.𝑤
        end
    end
end

function (op::Operator{:wΔMₙₛ})(ap::T,f::AbstractVector{Float64}) where T<:AbstractElement{:Poi1}
    𝓒 = ap.𝓒; ξ = ap.𝓖[1]
    N = get𝝭(ap,ξ)
    for (i,xᵢ) in enumerate(𝓒)
        I = xᵢ.id
        f[I] -= N[i]*ξ.ΔM
    end
end

function (op::Operator{:ΔMₙₛg})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; ξ = ap.𝓖[1]
    D = op.D
    ν = op.ν
    α = op.α
    Δn₁s₁ = ξ.Δn₁s₁
    Δn₁s₂n₂s₁ = ξ.Δn₁s₂n₂s₁
    Δn₂s₂ = ξ.Δn₂s₂
    D₁₁ = - D*(Δn₁s₁+Δn₂s₂*ν)
    D₁₂ = - D*(1-ν)*Δn₁s₂n₂s₁
    D₂₂ = - D*(Δn₁s₁*ν+Δn₂s₂)
    N,_,_,B₁₁,B₁₂,B₂₂ = get∇²𝝭(ap,ξ)
    for (i,xᵢ) in enumerate(𝓒)
        I = xᵢ.id
        ΔMₙₛᵢ = D₁₁*B₁₁[i] + D₁₂*B₁₂[i] + D₂₂*B₂₂[i]
        for (j,xⱼ) in enumerate(𝓒)
            J = xⱼ.id
            ΔMₙₛⱼ = D₁₁*B₁₁[j] + D₁₂*B₁₂[j] + D₂₂*B₂₂[j]
            k[I,J] += ΔMₙₛᵢ*N[j] + N[i]*ΔMₙₛⱼ + α*N[i]*N[j]
        end
        f[I] += (ΔMₙₛᵢ + α*N[i])*ξ.g
    end
end

function (op::Operator{:ΔM̃ₙₛg})(ap::T,k::AbstractMatrix{Float64},f::AbstractVector{Float64}) where T<:AbstractElement
    𝓒 = ap.𝓒; ξ = ap.𝓖[1]
    D = op.D
    ν = op.ν
    # α = op.α
    Δn₁s₁ = ξ.Δn₁s₁
    Δn₁s₂n₂s₁ = ξ.Δn₁s₂n₂s₁
    Δn₂s₂ = ξ.Δn₂s₂
    D₁₁ = - D*(Δn₁s₁+Δn₂s₂*ν)
    D₁₂ = - D*(1-ν)*Δn₁s₂n₂s₁
    D₂₂ = - D*(Δn₁s₁*ν+Δn₂s₂)
    N,_,_,B₁₁,B₁₂,B₂₂ = get∇²𝝭(ap,ξ)
    B̄₁₁,B̄₁₂,B̄₂₂ = get∇̄²𝝭(ap,ξ)
    for (i,xᵢ) in enumerate(𝓒)
        I = xᵢ.id
        ΔMₙₛᵢ = D₁₁*B₁₁[i] + D₁₂*B₁₂[i] + D₂₂*B₂₂[i]
        ΔM̄ₙₛᵢ = D₁₁*B̄₁₁[i] + D₁₂*B̄₁₂[i] + D₂₂*B̄₂₂[i]
        for (j,xⱼ) in enumerate(𝓒)
            J = xⱼ.id
            ΔMₙₛⱼ = D₁₁*B₁₁[j] + D₁₂*B₁₂[j] + D₂₂*B₂₂[j]
            # k[I,J] += ΔMₙₛᵢ*N[j] + N[i]*ΔMₙₛⱼ + α*N[i]*N[j]
            # k[I,J] += ΔMₙₛᵢ*N[j] + N[i]*ΔMₙₛⱼ + ΔM̄ₙₛᵢ*N[j]
            k[I,J] += ΔM̄ₙₛᵢ*N[j]
        end
        # f[I] += (ΔMₙₛᵢ + α*N[i])*ξ.g
        f[I] += (ΔMₙₛᵢ + ΔM̄ₙₛᵢ)*ξ.g
    end
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
## Imports

using BenchmarkTools
using Flux
using Graphs
using InferOpt
using MultiAgentPathFinding
using PythonCall
using ProgressMeter
using SparseArrays
using UnicodePlots
using Base.Threads;
nthreads();

## Test

rail_generators = pyimport("flatland.envs.rail_generators")
line_generators = pyimport("flatland.envs.line_generators")
rail_env = pyimport("flatland.envs.rail_env")

rail_generator = rail_generators.sparse_rail_generator(; max_num_cities=3)
line_generator = line_generators.sparse_line_generator()

pyenv = rail_env.RailEnv(;
    width=30,
    height=30,
    number_of_agents=20,
    rail_generator=rail_generator,
    line_generator=line_generator,
    random_seed=11,
)
pyenv.reset();

mapf = flatland_mapf(pyenv);

## Data generation

nb_instances = 100

instances = Vector{typeof(flatland_mapf(pyenv))}(undef, nb_instances);
@showprogress "Generating instances: " for k in 1:nb_instances
    pyenv.reset()
    instances[k] = flatland_mapf(pyenv)
end

solutions = Vector{Solution}(undef, nb_instances);
@threads for k in 1:nb_instances
    @info "Instance $k solved by thread $(threadid())"
    solutions[k] = large_neighborhood_search(instances[k]; N=10, steps=100, progress=false)
end

solutions_naive = Vector{Solution}(undef, nb_instances);
@showprogress for k in 1:nb_instances
    solutions_naive[k] = cooperative_astar(instances[k])
end

A = nb_agents(instances[1])
X = [edges_embedding(mapf) for mapf in instances];
Y = [solution_to_mat(solution, mapf) for (solution, mapf) in zip(solutions, instances)];

function maximizer(θ; mapf)
    edge_weights = -θ
    solution = independent_dijkstra(mapf, edge_weights)
    ŷ = solution_to_mat(solution, mapf)
    return ŷ
end

## Initialization

turn_negative(z) = -exp.(z) .- 1.0
repeat_agents(z::AbstractArray) = repeat(z; outer=(1, A))

encoder = Chain(Dense(size(X[1], 1), 1), vec, turn_negative, repeat_agents)
par = Flux.params(encoder)

fenchel_young_loss = FenchelYoungLoss(Perturbed(maximizer; ε=0.1, M=5));

opt = ADAGrad();

k = 1
θ = encoder(X[k])
maximizer(encoder(X[k]); mapf=instances[k])
fenchel_young_loss(encoder(X[k]), Y[k]; mapf=instances[k]) / nb_instances

## Training

nb_epochs = 100
losses = Float64[]
@showprogress "Training -" for epoch in 1:nb_epochs
    l = 0.0
    @showprogress "Epoch $epoch/$nb_epochs -" for k in 1:nb_instances
        gs = gradient(par) do
            l += fenchel_young_loss(encoder(X[k]), Y[k]; mapf=instances[k])
        end
        Flux.update!(opt, par, gs)
    end
    push!(losses, l / nb_instances)
end;

println(lineplot(losses))

encoder[1].weight

## Eval

solutions_pred = Vector{Solution}(undef, nb_instances);
@showprogress for k in 1:nb_instances
    mapf = instances[k]
    solution = independent_dijkstra(mapf, -encoder(X[k]))
    feasibility_search!(solution, mapf)
    solutions_pred[k] = solution
end

for k = 1:nb_instances
    mapf = instances[k]
    sol_naive = solutions_naive[k]
    sol_pred = solutions_pred[k]
    sol = solutions[k]
    @info "Solution comparison" flowtime(sol_naive, mapf) flowtime(sol_pred, mapf) flowtime(sol, mapf)
end


text:
using FactCheck
using Reactive
using Compat

step() = Reactive.run(1)
queue_size() = Base.n_avail(Reactive._messages)
number() = round(Int, rand()*1000)

## Basics

facts("Basic checks") do

    a = Signal(number())
    b = map(x -> x*x, a)

    context("map") do

        # Lift type
        #@fact typeof(b) --> Reactive.Lift{Int}

        # type conversion
        push!(a, 1.0)
        step()
        @fact value(b) --> 1
        # InexactError to be precise
        push!(a, 2.1, (n,x,err) -> @fact n --> a)
        step()

        @fact value(b) --> 1

        push!(a, number())
        step()
        @fact value(b) --> value(a)^2

        push!(a, -number())
        step()
        @fact value(b) --> value(a)^2

        ## Multiple inputs to Lift
        c = map(+, a, b, typ=Int)
        @fact value(c) --> value(a) + value(b)

        push!(a, number())
        step()
        @fact value(c) --> value(a) + value(b)

        push!(b, number())
        step()
        @fact value(c) --> value(a) + value(b)
    end

    context("merge") do

        ## Merge
        d = Signal(number())
        e = merge(d, b, a)

        # precedence to d
        @fact value(e) --> value(d)

        push!(a, number())
        step()
        # precedence to b over a -- a is older.
        @fact value(e) --> value(a)

        c = map(_->_, a) # Make a younger than b
        f = merge(d, c, b)
        push!(a, number())
        step()
        @fact value(f) --> value(c)
    end

    context("foldp") do

        ## foldl over time
        push!(a, 0)
        step()
        f = foldp(+, 0, a)
        nums = round(Int, rand(100)*1000)
        map(x -> begin push!(a, x); step() end, nums)

        @fact sum(nums) --> value(f)
    end

    context("filter") do
        # filter
        g = Signal(0)
        pred = x -> x % 2 != 0
        h = filter(pred, 1, g)
        j = filter(x -> x % 2 == 0, 1, g)

        @fact value(h) --> 1
        @fact value(j) --> 0

        push!(g, 2)
        step()
        @fact value(h) --> 1

        push!(g, 3)
        step()
        @fact value(h) --> 3
    end

    context("sampleon") do
        # sampleon
        g = Signal(0)

        push!(g, number())
        step()
        i = Signal(true)
        j = sampleon(i, g)
        # default value
        @fact value(j) --> value(g)
        push!(g, value(g)-1)
        step()
        @fact value(j) --> value(g)+1
        push!(i, true)
        step()
        @fact value(j) --> value(g)
    end

    context("droprepeats") do
        # droprepeats
        count = s -> foldp((x, y) -> x+1, 0, s)

        k = Signal(1)
        l = droprepeats(k)

        @fact value(l) --> value(k)
        push!(k, 1)
        step()
        @fact value(l) --> value(k)
        push!(k, 0)
        step()
        #println(l.value, " ", value(k))
        @fact value(l) --> value(k)

        m = count(k)
        n = count(l)

        seq = [1, 1, 1, 0, 1, 0, 1, 0, 0]
        map(x -> begin push!(k, x); step() end, seq)

        @fact value(m) --> length(seq)
        @fact value(n) --> 6
    end

    context("filterwhen") do
        # filterwhen
        b = Signal(false)
        n = Signal(1)
        dw = filterwhen(b, 0, n)
        @fact value(dw) --> 0
        push!(n, 2)
        step()
        @fact value(dw) --> 0
        push!(b, true)
        step()
        @fact value(dw) --> 0
        push!(n, 1)
        step()
        @fact value(dw) --> 1
        push!(n, 2)
        step()
        @fact value(dw) --> 2
        dw = filterwhen(b, 0, n)
        @fact value(dw) --> 2
    end

    context("push! inside push!") do
        a = Signal(0)
        b = Signal(1)
        Reactive.preserve(map(x -> push!(a, x), b))

        @fact value(a) --> 0

        step()
        @fact value(a) --> 1

        push!(a, 2)
        step()
        @fact value(a) --> 2
        @fact value(b) --> 1

        push!(b, 3)
        step()
        @fact value(b) --> 3
        @fact value(a) --> 2

        step()
        @fact value(a) --> 3
    end

    context("previous") do
        x = Signal(0)
        y = previous(x)
        @fact value(y) --> 0

        push!(x, 1)
        step()

        @fact value(y) --> 0

        push!(x, 2)
        step()

        @fact value(y) --> 1

        push!(x, 3)
        step()

        @fact value(y) --> 2
        @fact queue_size() --> 0
    end

    context("delay") do
        x = Signal(0)
        y = delay(x)
        @fact value(y) --> 0

        push!(x, 1)
        step()

        @fact queue_size() --> 1
        @fact value(y) --> 0

        step()
        @fact value(y) --> 1
    end
end


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

# %%
# privateとpublicの区別は`Base.propertynames`の定義で行える。
@doc propertynames

# %%
module O

"""
抽象平面ベクトルの型
"""
abstract type AbstractPlanarVector{T} end

"""タプル版和集合(重複を削除)"""
tupleunion(x) = Tuple(unique(x))
tupleunion(x, y) = Tuple(unique((x..., y...)))
tupleunion(x, y, z...) = tupleunion(tupleunion(x, y), z...)

"""
抽象平面ベクトルのpublic propertiesは x と y

タブ補完でもx, yしか表示されなくなる。
"""
function Base.propertynames(p::AbstractPlanarVector, private=false)
    public_properties = (:x, :y)
    private && return tupleunion(public_properties, fieldnames(p))
    public_properties
end

"""
抽象平面ベクトルは (x, y) の形式で表示される
"""
Base.show(io::IO, p::AbstractPlanarVector) = show(io, (p.x, p.y))

"""
AbstractPlanarVector{T}型の抽象平面ベクトルの成分の型は T
"""
Base.eltype(x::AbstractPlanarVector{T}) where T = T

"""
平面ベクトルの型

x, yをfieldsに持ち、それらがそのままpublic propertiesになる。
"""
struct PlanarVector{T} <: AbstractPlanarVector{T} x::T; y::T end

"""
標準基底の型

x, yをfieldsとして持たず、型のパラメータ i で何番目の基底ベクトルであるかを指定。
"""
struct CanonBasis{T, i} <: AbstractPlanarVector{T} end
function CanonBasis{T}(i) where T
    @assert i == 1 || i == 2
    CanonBasis{T, i}()
end
CanonBasis(i) = CanonBasis{Int}(i)

"""
標準基底のpublic propertiesの取得メソッド

x, yをfieldsとして持たないので定義してやる必要がある。

* i == 1 ⇒ (x, y) = (1, 0)
* i == 2 ⇒ (x, y) = (0, 1)
"""
function Base.getproperty(::CanonBasis{T, i}, name::Symbol) where {T, i}
    name === :x && return ifelse(i == 1, one(T), zero(T))
    name === :y && return ifelse(i == 2, one(T), zero(T))
    error("type CanonBasis has no public property $name")
end

########## ベクトルの演算達

# 抽象平面ベクトル p, q について、 +p, -p, p + q, p - q を定義
for op in (:+, :-)
    @eval Base.$op(p::AbstractPlanarVector) = PlanarVector($op(p.x), $op(p.y))
    @eval Base.$op(p::AbstractPlanarVector, q::AbstractPlanarVector) =
        PlanarVector($op(p.x, q.x), $op(p.y, q.y))
end

# 抽象平面ベクトル p とスカラー a について、 a * p, p * a, a \ p, p / a を定義
Base.:*(a, p::AbstractPlanarVector) = PlanarVector(a * p.x, a * p.y)
Base.:*(p::AbstractPlanarVector, a) = PlanarVector(p.x * a, p.y * a)
Base.:\(a, p::AbstractPlanarVector) = PlanarVector(a \ p.x, a \ p.y)
Base.:/(p::AbstractPlanarVector, a) = PlanarVector(p.x / a, p.y / a)

# LinearAlgebra の dot の定義の準備
using LinearAlgebra

##### 抽象ベクトルの内積を定義

# まず、一般的な定義式で内積を定義
LinearAlgebra.dot(p::AbstractPlanarVector, q::AbstractPlanarVector) =
    conj(p.x) * q.x + conj(p.y) * q.y

# 次に、標準基底の場合に特殊化した定義を行う。多重ディスパッチを本質的に使っている。
LinearAlgebra.dot(p::CanonBasis{T, i}, q::AbstractPlanarVector{U}) where {T, i, U} =
    (P = promote_type(T, U); P(i == 1 ? q.x : q.y))
LinearAlgebra.dot(p::AbstractPlanarVector{U}, q::CanonBasis{T, i}) where {U, T, i} =
    (P = promote_type(T, U); P(conj(i == 1 ? p.x : p.y)))

# Juliaの多重ディスパッチでは上の2つの場合のintersectionの場合の定義もしておく必要がある。
LinearAlgebra.dot(p::CanonBasis{T, i}, q::CanonBasis{U, j}) where {T, i, U, j} =
    (P = promote_type(T, U); i == j ? one(P) : zero(P))

end

# %%
p = O.PlanarVector(2+im, 3+im)

# %%
q = O.PlanarVector(-5.0, 10.0)

# %%
p + q

# %%
eltype(p + q)

# %%
-3p + 4q

# %%
O.dot(p, q)

# %%
e1 = O.CanonBasis{Float64}(1)

# %%
e2 = O.CanonBasis(2)

# %%
O.dot(e1, e2)

# %%
O.dot(O.CanonBasis{ComplexF64}(1), e1)

# %%
e1 - e2

# %%
3*e1 - e2/2

# %%
p.x, p.y

# %%
e1.x, e1.y

# %%
getx(p::O.AbstractPlanarVector) = p.x
gety(p::O.AbstractPlanarVector) = p.y

# %%
getx(p), gety(p)

# %%
@code_typed getx(p)

# %%
@code_typed gety(q)

# %%
@code_typed getx(e2)

# %%
@code_typed gety(e2)

# %%
@code_typed O.dot(p, q)

# %%
@code_typed O.dot(e1, q)

# %%
@code_typed O.dot(p, e2)

# %%
@code_llvm debuginfo=:none O.dot(p, e2)

# %%
@code_typed O.dot(e2, e2)

# %%
@code_typed O.dot(e1, e2)

# %%
@code_llvm debuginfo=:none O.dot(e1, e2)

# %%
using BenchmarkTools
n = 10^6
v = [O.PlanarVector(randn(2)...) for _ in 1:n]
a = O.PlanarVector(1, 0)
e = O.CanonBasis(1)

@show sum(Base.Fix1(O.dot, a), v) == sum(Base.Fix1(O.dot, e), v)
@btime sum($(Base.Fix1(O.dot, a)), $v)
@btime sum($(Base.Fix1(O.dot, e)), $v)

# %%
@code_native debuginfo=:none O.dot(a, q)

# %%
@code_native debuginfo=:none O.dot(e, q)

# %%

max_stars_repo_path: 0005/PlanarVector example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/PlanarVector example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/PlanarVector example.jl

max_forks_repo_name: genkuroki/public


text:
export iauNut00a
"""
Nutation, IAU 2000A model (MHB2000 luni-solar and planetary nutation
with free core nutation omitted).

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  canonical model.

Given:
   date1,date2   double   TT as a 2-part Julian Date (Note 1)

Returned:
   dpsi,deps     double   nutation, luni-solar + planetary (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The nutation components in longitude and obliquity are in radians
   and with respect to the equinox and ecliptic of date.  The
   obliquity at J2000.0 is assumed to be the Lieske et al. (1977)
   value of 84381.448 arcsec.

   Both the luni-solar and planetary nutations are included.  The
   latter are due to direct planetary nutations and the
   perturbations of the lunar and terrestrial orbits.

   3. The function computes the MHB2000 nutation series with the
   associated corrections for planetary nutations.  It is an
   implementation of the nutation part of the IAU 2000A precession-
   nutation model, formally adopted by the IAU General Assembly in
   2000, namely MHB2000 (Mathews et al. 2002), but with the free
   core nutation (FCN - see Note 4) omitted.

   4. The full MHB2000 model also contains contributions to the
   nutations in longitude and obliquity due to the free-excitation
   of the free-core-nutation during the period 1979-2000.  These FCN
   terms, which are time-dependent and unpredictable, are NOT
   included in the present function and, if required, must be
   independently computed.  With the FCN corrections included, the
   present function delivers a pole which is at current epochs
   accurate to a few hundred microarcseconds.  The omission of FCN
   introduces further errors of about that size.

   5. The present function provides classical nutation.  The MHB2000
   algorithm, from which it is adapted, deals also with (i) the
   offsets between the GCRS and mean poles and (ii) the adjustments
   in longitude and obliquity due to the changed precession rates.
   These additional functions, namely frame bias and precession
   adjustments, are supported by the SOFA functions iauBi00  and
   iauPr00.

   6. The MHB2000 algorithm also provides "total" nutations, comprising
   the arithmetic sum of the frame bias, precession adjustments,
   luni-solar nutation and planetary nutation.  These total
   nutations can be used in combination with an existing IAU 1976
   precession implementation, such as iauPmat76,  to deliver GCRS-
   to-true predictions of sub-mas accuracy at current dates.
   However, there are three shortcomings in the MHB2000 model that
   must be taken into account if more accurate or definitive results
   are required (see Wallace 2002):

      (i) The MHB2000 total nutations are simply arithmetic sums,
         yet in reality the various components are successive Euler
         rotations.  This slight lack of rigor leads to cross terms
         that exceed 1 mas after a century.  The rigorous procedure
         is to form the GCRS-to-true rotation matrix by applying the
         bias, precession and nutation in that order.

      (ii) Although the precession adjustments are stated to be with
         respect to Lieske et al. (1977), the MHB2000 model does
         not specify which set of Euler angles are to be used and
         how the adjustments are to be applied.  The most literal
         and straightforward procedure is to adopt the 4-rotation
         epsilon_0, psi_A, omega_A, xi_A option, and to add DPSIPR
         to psi_A and DEPSPR to both omega_A and eps_A.

      (iii) The MHB2000 model predates the determination by Chapront
         et al. (2002) of a 14.6 mas displacement between the
         J2000.0 mean equinox and the origin of the ICRS frame.  It
         should, however, be noted that neglecting this displacement
         when calculating star coordinates does not lead to a
         14.6 mas change in right ascension, only a small second-
         order distortion in the pattern of the precession-nutation
         effect.

   For these reasons, the SOFA functions do not generate the "total
   nutations" directly, though they can of course easily be
   generated by calling iauBi00, iauPr00 and the present function
   and adding the results.

   7. The MHB2000 model contains 41 instances where the same frequency
   appears multiple times, of which 38 are duplicates and three are
   triplicates.  To keep the present code close to the original MHB
   algorithm, this small inefficiency has not been corrected.

Called:
   iauFal03     mean anomaly of the Moon
   iauFaf03     mean argument of the latitude of the Moon
   iauFaom03    mean longitude of the Moon's ascending node
   iauFame03    mean longitude of Mercury
   iauFave03    mean longitude of Venus
   iauFae03     mean longitude of Earth
   iauFama03    mean longitude of Mars
   iauFaju03    mean longitude of Jupiter
   iauFasa03    mean longitude of Saturn
   iauFaur03    mean longitude of Uranus
   iauFapa03    general accumulated precession in longitude

References:

   Chapront, J., Chapront-Touze, M. & Francou, G. 2002,
   Astron.Astrophys. 387, 700

   Lieske, J.H., Lederle, T., Fricke, W. & Morando, B. 1977,
   Astron.Astrophys. 58, 1-16

   Mathews, P.M., Herring, T.A., Buffet, B.A. 2002, J.Geophys.Res.
   107, B4.  The MHB_2000 code itself was obtained on 9th September
   2002 from ftp//maia.usno.navy.mil/conv2000/chapter5/IAU2000A.

   Simon, J.-L., Bretagnon, P., Chapront, J., Chapront-Touze, M.,
   Francou, G., Laskar, J. 1994, Astron.Astrophys. 282, 663-683

   Souchay, J., Loysel, B., Kinoshita, H., Folgueira, M. 1999,
   Astron.Astrophys.Supp.Ser. 135, 111

   Wallace, P.T., "Software for Implementing the IAU 2000
   Resolutions", in IERS Workshop 5.1 (2002)

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauNut00a(double date1, double date2, double *dpsi, double *deps)

function iauNut00a(date1::Real, date2::Real)
   # Allocate return values
   ref_dpsi = Ref{Float64}(0.0)
   ref_deps = Ref{Float64}(0.0)

   ccall((:iauNut00a, libsofa_c), Cvoid,
        (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
        convert(Float64, date1), convert(Float64, date2),
        ref_dpsi, ref_deps)

   return ref_dpsi[], ref_deps[]
end


text:
export System,
    JacobianConfig, JacobianDiffResult,
    evaluate_and_jacobian!, evaluate_and_jacobian

"""
    System(polys [, variables])

Construct a system of polynomials from the given polynomials `polys`.
"""
struct System{T}
    polys::Vector{Polynomial{T}}
end

function System(polys::Vector{<:MP.AbstractPolynomialLike}, variables=_variables(polys))
    System([Polynomial(p, variables) for p in polys])
end
_variables(polys) = sort!(union(Iterators.flatten(MP.variables.(polys))), rev=true)

Base.length(F::System) = length(F.polys)
Base.getindex(F::System, i) = getindex(F.polys, i)

"""
    nvariables(F::System)

Returns the number of variables of `F`.
"""
nvariables(F::System) = size(exponents(F[1]), 1)

"""

    JacobianConfig(F::Vector{Polynomial{T}}, [x::AbstractVector{S}])

A data structure with which the jacobian of a `Vector` `F` of `Polynomial`s can be
evaluated efficiently. Note that `x` is only used to determine the
output type of `F(x)`.

    JacobianConfig(F::Vector{Polynomial{T}}, [S])

Instead of a vector `x` a type can also be given directly.
"""
mutable struct JacobianConfig{T1, T2}
    polys::Vector{PolyConfig{T1}}
    differences::Matrix{UInt8}
    differences_values::Matrix{T2}
end

function JacobianConfig(f::System{T}, ::AbstractArray{S}) where {T, S}
    JacobianConfig(f, S)
end
JacobianConfig(f::System{T}) where T = JacobianConfig(f, T)
function JacobianConfig(F::System{T}, ::Type{S}) where {T, S}
    diffs, diffs_values, big_lookups, reduced_exponents = differences(F.polys, S)
    polys = broadcast(PolyConfig, F.polys, reduced_exponents, big_lookups, S)

    JacobianConfig(polys, diffs, diffs_values)
end

"""
    config(F::System, x)

Construct a `JacobianConfig` for the evaluation of `F` with values like `x`.
"""
function config(f::System{T}, x::AbstractVector{S}) where {S, T}
    JacobianConfig(f, typeof(one(T) * one(S) + one(T) * one(S)))
end

function differences(F::Vector{<:Polynomial}, ::Type{S}) where {S}
    reduced_exponents = map(F) do f
        convert.(UInt16, max.(f.exponents .- 1, 0))
    end
    differences, big_lookups = computetables(reduced_exponents)
    differences_values = similar(differences, S)

    differences, differences_values, big_lookups, reduced_exponents
end

"""
    evaluate(F, x, cfg::JacobianConfig [, precomputed=false])

Evaluate the system `F` at `x` using the precomputated values in `cfg`.
Note that this is usually signifcant faster than `map(f -> evaluate(f, x), F)`.
The return vector is constructed using `similar(x, T)`.

### Example
```julia
cfg = JacobianConfig(F)
evaluate(F, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
`x`.
"""
function evaluate(G::System, x::AbstractVector, cfg::JacobianConfig{T}, precomputed=false) where T
    evaluate!(similar(x, T, length(G)), G, x, cfg, precomputed)
end

"""
    evaluate!(u, F, x, cfg::JacobianConfig [, precomputed=false])

Evaluate the system `F` at `x` using the precomputated values in `cfg`
and store the result in `u`.
Note that this is usually signifcant faster than `map!(u, f -> evaluate(f, x), F)`.

### Example
```julia
cfg = JacobianConfig(F)
evaluate!(u, F, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
`x`.
"""
function evaluate!(u::AbstractVector, G::System, x::AbstractVector, cfg::JacobianConfig, precomputed=false)
    if !precomputed
        fillvalues!(cfg.differences_values, x, cfg.differences)
        for i=1:length(cfg.polys)
            fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
            u[i] = _evaluate(x, cfg.polys[i])
        end
    else
        for i=1:length(cfg.polys)
            u[i] = _evaluate(x, cfg.polys[i])
        end
    end
    u
end

"""
    jacobian(u, F, x, cfg::JacobianConfig [, precomputed=false])

Evaluate the jacobian of `F` at `x` using the precomputated values in `cfg`. The return
matrix is constructed using `similar(x, T, m, n)`.

### Example
```julia
cfg = JacobianConfig(F)
jacobian(F, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
`x`.
"""
function jacobian(g::System, x::AbstractVector, cfg::JacobianConfig{T}, precomputed=false) where T
    u = similar(x, T, (length(g), length(x)))
    jacobian!(u, g, x, cfg, precomputed)
    u
end

"""
    jacobian!(u, F, x, cfg::JacobianConfig [, precomputed=false])

Evaluate the jacobian of `F` at `x` using the precomputated values in `cfg`
and store the result in `u`.

### Example
```julia
cfg = JacobianConfig(F)
jacobian!(u, F, x, cfg)
```

With `precomputed=true` we rely on the previous intermediate results in `cfg`. Therefore
the result is only correct if you previouls called `evaluate`, or `jacobian` with the same
`x`.
"""
function jacobian!(u::AbstractMatrix, G::System, x::AbstractVector, cfg::JacobianConfig, precomputed=false)
    if !precomputed
        fillvalues!(cfg.differences_values, x, cfg.differences)
        for i=1:length(G)
            fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
            gradient_row!(u, x, cfg.polys[i], i)
        end
    else
        for i=1:length(G)
            gradient_row!(u, x, cfg.polys[i], i)
        end
    end
    u
end

"""
    evaluate_and_jacobian!(u, U, F, x, cfg::JacobianConfig)

Evaluate `F` and its Jacobian at `x` using the precomputated values in `cfg`
and store the result in `u` and `U` (Jacobian).

### Example
```julia
evaluate_and_jacobian!(u, U, F, x, config(F, x))
```
"""
function evaluate_and_jacobian!(u::AbstractVector, U::AbstractMatrix, G::System, x::AbstractVector, cfg::JacobianConfig)
    fillvalues!(cfg.differences_values, x, cfg.differences)
    for i=1:length(G)
        fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
        gradient_row!(U, x, cfg.polys[i], i)
        u[i] = _evaluate(x, cfg.polys[i])
    end
end

function evaluate_and_jacobian(G::System, x::AbstractVector, cfg::JacobianConfig{T}) where T
    u = similar(x, T, length(G))
    U = similar(x, T, (length(G), length(x)))

    evaluate_and_jacobian!(u, U, G, x, cfg)
    u, U
end

"""
    JacobianDiffResult(cfg::GradientConfig)

During the computation of the jacobian ``J_F(x)`` we compute nearly everything we need for the evaluation of
``F(x)``. `JacobianDiffResult` allocates memory to hold both values.
This structure also signals `jacobian!` to store ``F(x)`` and ``J_F(x)``.

### Example

```julia
cfg = JacobianConfig(F, x)
r = JacobianDiffResult(cfg)
jacobian!(r, F, x, cfg)

value(r) == map(f -> f(x), F)
jacobian(r) == jacobian(F, x, cfg)
```

    JacobianDiffResult(value::AbstractVector, jacobian::AbstractMatrix)

Allocate the memory to hold the value and the jacobian by yourself.
"""
mutable struct JacobianDiffResult{T, AV<:AbstractVector{T}, AM<:AbstractMatrix{T}}
    value::AV
    jacobian::AM
end

function JacobianDiffResult(cfg::JacobianConfig{T}) where T
    JacobianDiffResult{T, Vector{T}, Matrix{T}}(
        zeros(T, length(cfg.polys)),
        zeros(T, length(cfg.polys), size(cfg.differences, 1)))
end

@static if VERSION < v"0.7-"
    function JacobianDiffResult(value::AbstractVector{T}, jacobian::AbstractMatrix{T}) where T
        JacobianDiffResult{T, typeof(value), typeof(jacobian)}(value, jacobian)
    end
end

value(r::JacobianDiffResult) = r.value
jacobian(r::JacobianDiffResult) = r.jacobian

"""
    jacobian!(r::JacobianDiffResult, F, x, cfg::JacobianConfig)

Compute ``F(x)`` and the jacobian of `F` at `x` at once using the precomputated values in `cfg`
and store thre result in `r`. This is faster than computing both values separetely.

### Example
```julia
cfg = GradientConfig(g)
r = GradientDiffResult(cfg)
gradient!(r, g, x, cfg)

value(r) == g(x)
gradient(r) == gradient(g, x, cfg)
```
"""
function jacobian!(r::JacobianDiffResult, G::System, x::AbstractVector, cfg::JacobianConfig{T}) where T
    fillvalues!(cfg.differences_values, x, cfg.differences)
    for i=1:length(G)
        fillreduced_values!(cfg.polys[i], G[i], x, cfg.differences_values)
        gradient_row!(r.jacobian, x, cfg.polys[i], i)
        r.value[i] = _evaluate(x, cfg.polys[i])
    end

    r
end


text:
module ToyAD
  using Blink, JSON
  import Base.+, Base.(.+), Base.*, Base.(.*), Base.show

  # Automatic Differentiation
  include("Utils.jl")
  include("AD.jl")
  include("ADPlot.jl")
  include("ADOperations.jl")
  include("Backprop.jl")

  # types
  export AD

  # plotting graph
  export plot

  # operations on AD type
  export +, .+, *, .*, relu, sigmoid, tanh, softmaxLoss

  # Backprop
  export backprop, resetGrad

end # module


text:
@testset "gamma and friends" begin
    @testset "digamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            @test digamma(convert(elty, 9)) ≈ convert(elty, 2.140641477955609996536345)
            @test digamma(convert(elty, 2.5)) ≈ convert(elty, 0.7031566406452431872257)
            @test digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041107679516822)
            @test digamma(convert(elty, 7e-4)) ≈ convert(elty, -1429.147493371120205005198)
            @test digamma(convert(elty, 7e-5)) ≈ convert(elty, -14286.29138623969227538398)
            @test digamma(convert(elty, 7e-6)) ≈ convert(elty, -142857.7200612932791081972)
            @test digamma(convert(elty, 2e-6)) ≈ convert(elty, -500000.5772123750382073831)
            @test digamma(convert(elty, 1e-6)) ≈ convert(elty, -1000000.577214019968668068)
            @test digamma(convert(elty, 7e-7)) ≈ convert(elty, -1428572.005785942019703646)
            @test digamma(convert(elty, -0.5)) ≈ convert(elty, .03648997397857652055902367)
            @test digamma(convert(elty, -1.1)) ≈ convert(elty,  10.15416395914385769902271)

            @test digamma(convert(elty, 0.1)) ≈ convert(elty, -10.42375494041108)
            @test digamma(convert(elty, 1/2)) ≈ convert(elty, -γ - log(4))
            @test digamma(convert(elty, 1)) ≈ convert(elty, -γ)
            @test digamma(convert(elty, 2)) ≈ convert(elty, 1 - γ)
            @test digamma(convert(elty, 3)) ≈ convert(elty, 3/2 - γ)
            @test digamma(convert(elty, 4)) ≈ convert(elty, 11/6 - γ)
            @test digamma(convert(elty, 5)) ≈ convert(elty, 25/12 - γ)
            @test digamma(convert(elty, 10)) ≈ convert(elty, 7129/2520 - γ)
        end
    end

    @testset "trigamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            @test trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
            @test trigamma(convert(elty, 0.1)) ≈ convert(elty, 101.433299150792758817)
            @test trigamma(convert(elty, 1/2)) ≈ convert(elty, π^2/2)
            @test trigamma(convert(elty, 1)) ≈ convert(elty, π^2/6)
            @test trigamma(convert(elty, 2)) ≈ convert(elty, π^2/6 - 1)
            @test trigamma(convert(elty, 3)) ≈ convert(elty, π^2/6 - 5/4)
            @test trigamma(convert(elty, 4)) ≈ convert(elty, π^2/6 - 49/36)
            @test trigamma(convert(elty, 5)) ≈ convert(elty, π^2/6 - 205/144)
            @test trigamma(convert(elty, 10)) ≈ convert(elty, π^2/6 - 9778141/6350400)
        end
    end

    @testset "invdigamma" begin
        @testset "$elty" for elty in (Float32, Float64)
            for val in [0.001, 0.01, 0.1, 1.0, 10.0]
                @test abs(invdigamma(digamma(convert(elty, val))) - convert(elty, val)) < 1e-8
            end
        end
        @test abs(invdigamma(2)) == abs(invdigamma(2.))
    end

    @testset "polygamma" begin
        @test polygamma(20, 7.) ≈ -4.644616027240543262561198814998587152547
        @test polygamma(20, Float16(7.)) ≈ -4.644616027240543262561198814998587152547
    end

    @testset "eta" begin
        @test eta(1) ≈ log(2)
        @test eta(2) ≈ pi^2/12
        @test eta(Float32(2)) ≈ eta(2)
        @test eta(Complex{Float32}(2)) ≈ eta(2)
    end
    @testset "gamma, loggamma, logabsgamma (complex argument)" begin
        if Base.Math.libm == "libopenlibm"
            @test gamma.(Float64[1:25;]) == gamma.(1:25)
        else
            @test gamma.(Float64[1:25;]) ≈ gamma.(1:25)
        end
        for elty in (Float32, Float64)
            @test gamma(convert(elty,1/2)) ≈ convert(elty,sqrt(π))
            @test gamma(convert(elty,-1/2)) ≈ convert(elty,-2sqrt(π))
            @test logabsgamma(convert(elty,-1/2))[1] ≈ convert(elty,log(abs(gamma(-1/2))))
        end
        @test loggamma(1.4+3.7im) ≈ -3.7094025330996841898 + 2.4568090502768651184im
        @test loggamma(1.4+3.7im) ≈ log(gamma(1.4+3.7im))
        @test loggamma(-4.2+0im) ≈ logabsgamma(-4.2)[1] - 5pi*im
        @test SpecialFunctions.factorial(3.0) == gamma(4.0) == factorial(3)
        for x in (3.2, 2+1im, 3//2, 3.2+0.1im)
            @test SpecialFunctions.factorial(x) == gamma(1+x)
        end
        @test logfactorial(0) == logfactorial(1) == 0
        @test logfactorial(2) == loggamma(3)
        # Ensure that the domain of logfactorial matches that of factorial (issue #21318)
        @test_throws DomainError logfactorial(-3)
        @test_throws DomainError loggamma(-4.2)
        @test_throws MethodError logfactorial(1.0)
    end

    # loggamma & logabsgamma test cases (from Wolfram Alpha)
    @test loggamma(-300im) ≅ -473.17185074259241355733179182866544204963885920016823743 - 1410.3490664555822107569308046418321236643870840962522425im
    @test loggamma(3.099) ≅ loggamma(3.099+0im) ≅ 0.786413746900558058720665860178923603134125854451168869796
    @test loggamma(1.15) ≅ loggamma(1.15+0im) ≅ -0.06930620867104688224241731415650307100375642207340564554
    @test logabsgamma(0.89)[1] ≅ loggamma(0.89+0im) ≅ 0.074022173958081423702265889979810658434235008344573396963
    @test loggamma(0.91) ≅ loggamma(0.91+0im) ≅ 0.058922567623832379298241751183907077883592982094770449167
    @test loggamma(0.01) ≅ loggamma(0.01+0im) ≅ 4.599479878042021722513945411008748087261001413385289652419
    @test loggamma(-3.4-0.1im) ≅ -1.1733353322064779481049088558918957440847715003659143454 + 12.331465501247826842875586104415980094316268974671819281im
    @test loggamma(-13.4-0.1im) ≅ -22.457344044212827625152500315875095825738672314550695161 + 43.620560075982291551250251193743725687019009911713182478im
    @test loggamma(-13.4+0.0im) ≅ conj(loggamma(-13.4-0.0im)) ≅ -22.404285036964892794140985332811433245813398559439824988 - 43.982297150257105338477007365913040378760371591251481493im
    @test loggamma(-13.4+8im) ≅ -44.705388949497032519400131077242200763386790107166126534 - 22.208139404160647265446701539526205774669649081807864194im
    @test logabsgamma(1+exp2(-20))[1] ≅ loggamma(1+exp2(-20)+0im) ≅ -5.504750066148866790922434423491111098144565651836914e-7
    @test loggamma(1+exp2(-20)+exp2(-19)*im) ≅ -5.5047799872835333673947171235997541985495018556426e-7 - 1.1009485171695646421931605642091915847546979851020e-6im
    @test loggamma(-300+2im) ≅ -1419.3444991797240659656205813341478289311980525970715668 - 932.63768120761873747896802932133229201676713644684614785im
    @test loggamma(300+2im) ≅ 1409.19538972991765122115558155209493891138852121159064304 + 11.4042446282102624499071633666567192538600478241492492652im
    @test loggamma(1-6im) ≅ -7.6099596929506794519956058191621517065972094186427056304 - 5.5220531255147242228831899544009162055434670861483084103im
    @test loggamma(1-8im) ≅ -10.607711310314582247944321662794330955531402815576140186 - 9.4105083803116077524365029286332222345505790217656796587im
    @test loggamma(1+6.5im) ≅ conj(loggamma(1-6.5im)) ≅ -8.3553365025113595689887497963634069303427790125048113307 + 6.4392816159759833948112929018407660263228036491479825744im
    @test loggamma(1+1im) ≅ conj(loggamma(1-1im)) ≅ -0.6509231993018563388852168315039476650655087571397225919 - 0.3016403204675331978875316577968965406598997739437652369im
    @test loggamma(-pi*1e7 + 6im) ≅ -5.10911758892505772903279926621085326635236850347591e8 - 9.86959420047365966439199219724905597399295814979993e7im
    @test loggamma(-pi*1e7 + 8im) ≅ -5.10911765175690634449032797392631749405282045412624e8 - 9.86959074790854911974415722927761900209557190058925e7im
    @test loggamma(-pi*1e14 + 6im) ≅ -1.0172766411995621854526383224252727000270225301426e16 - 9.8696044010873714715264929863618267642124589569347e14im
    @test loggamma(-pi*1e14 + 8im) ≅ -1.0172766411995628137711690403794640541491261237341e16 - 9.8696044010867038531027376655349878694397362250037e14im
    @test loggamma(2.05 + 0.03im) ≅ conj(loggamma(2.05 - 0.03im)) ≅ 0.02165570938532611215664861849215838847758074239924127515 + 0.01363779084533034509857648574107935425251657080676603919im
    @test loggamma(2+exp2(-20)+exp2(-19)*im) ≅ 4.03197681916768997727833554471414212058404726357753e-7 + 8.06398296652953575754782349984315518297283664869951e-7im

    @testset "loggamma for non-finite arguments" begin
        @test loggamma(Inf + 0im) === Inf + 0im
        @test loggamma(Inf - 0.0im) === Inf - 0.0im
        @test loggamma(Inf + 1im) === Inf + Inf*im
        @test loggamma(Inf - 1im) === Inf - Inf*im
        @test loggamma(-Inf + 0.0im) === -Inf - Inf*im
        @test loggamma(-Inf - 0.0im) === -Inf + Inf*im
        @test loggamma(Inf*im) === -Inf + Inf*im
        @test loggamma(-Inf*im) === -Inf - Inf*im
        @test loggamma(Inf + Inf*im) === loggamma(NaN + 0im) === loggamma(NaN*im) === NaN + NaN*im
    end

    @testset "Other float types" begin
        struct NotAFloat <: AbstractFloat
        end
        let x = one(Float16)
            @test gamma(x) ≈ one(Float16)
            @test gamma(x) isa Float16
            @test loggamma(x) ≈ zero(Float16)
            @test loggamma(x) isa Float16
        end
        @test_throws MethodError gamma(NotAFloat())
        @test_throws MethodError logabsgamma(NotAFloat())
        @test_throws MethodError loggamma(NotAFloat())
    end
end

@testset "zeta" begin
    @test zeta(0) ≈ -0.5
    @test zeta(2) ≈ pi^2/6
    @test zeta(Complex{Float32}(2)) ≈ zeta(2)
    @test zeta(4) ≈ pi^4/90
    @test zeta(1,Float16(2.)) ≈ zeta(1,2.)
    @test zeta(1.,Float16(2.)) ≈ zeta(1,2.)
    @test zeta(Float16(1.),Float16(2.)) ≈ zeta(1,2.)
    @test isnan(zeta(NaN))
    @test isnan(zeta(1.0e0))
    @test isnan(zeta(1.0f0))
    @test isnan(zeta(complex(0,Inf)))
    @test isnan(zeta(complex(-Inf,0)))
end

#(compared to Wolfram Alpha)
@testset "digamma, trigamma, polygamma & zeta" begin
    for x in -10.2:0.3456:50
        @test 1e-12 > relerr(digamma(x+0im), digamma(x))
    end
    @test digamma(7+0im) ≅ 1.872784335098467139393487909917597568957840664060076401194232
    @test digamma(7im) ≅ 1.94761433458434866917623737015561385331974500663251349960124 + 1.642224898223468048051567761191050945700191089100087841536im
    @test digamma(-3.2+0.1im) ≅ 4.65022505497781398615943030397508454861261537905047116427511+2.32676364843128349629415011622322040021960602904363963042380im
    @test trigamma(8+0im) ≅ 0.133137014694031425134546685920401606452509991909746283540546
    @test trigamma(8im) ≅ -0.0078125000000000000029194973110119898029284994355721719150 - 0.12467345030312762782439017882063360876391046513966063947im
    @test trigamma(-3.2+0.1im) ≅ 15.2073506449733631753218003030676132587307964766963426965699+15.7081038855113567966903832015076316497656334265029416039199im
    @test polygamma(2, 8.1+0im) ≅ -0.01723882695611191078960494454602091934457319791968308929600
    @test polygamma(30, 8.1+2im) ≅ -2722.8895150799704384107961215752996280795801958784600407589+6935.8508929338093162407666304759101854270641674671634631058im
    @test polygamma(3, 2.1+1im) ≅ 0.00083328137020421819513475400319288216246978855356531898998-0.27776110819632285785222411186352713789967528250214937861im
    @test 1e-11 > relerr(polygamma(3, -4.2 + 2im),-0.0037752884324358856340054736472407163991189965406070325067-0.018937868838708874282432870292420046797798431078848805822im)
    @test polygamma(13, 5.2 - 2im) ≅ 0.08087519202975913804697004241042171828113370070289754772448-0.2300264043021038366901951197725318713469156789541415899307im
    @test 1e-11 > relerr(polygamma(123, -47.2 + 0im), 5.7111648667225422758966364116222590509254011308116701029e291)
    @test zeta(4.1+0.3im, -3.2+0.1im) ≅ -281.34474134962502296077659347175501181994490498591796647 + 286.55601240093672668066037366170168712249413003222992205im
    @test zeta(4.1+0.3im, 3.2+0.1im) ≅ 0.0121197525131633219465301571139288562254218365173899270675-0.00687228692565614267981577154948499247518236888933925740902im
    @test zeta(4.1, 3.2+0.1im) ≅ 0.0137637451187986846516125754047084829556100290057521276517-0.00152194599531628234517456529686769063828217532350810111482im
    @test 1e-12 > relerr(zeta(1.0001, -4.5e2+3.2im), 10003.765660925877260544923069342257387254716966134385170 - 0.31956240712464746491659767831985629577542514145649468090im)
    @test zeta(3.1,-4.2) ≅ zeta(3.1,-4.2+0im) ≅ 149.7591329008219102939352965761913107060718455168339040295
    @test 1e-15 > relerr(zeta(3.1+0im,-4.2), zeta(3.1,-4.2+0im))
    @test zeta(3.1,4.2) ≅ 0.029938344862645948405021260567725078588893266227472565010234
    @test zeta(27, 3.1) ≅ 5.413318813037879056337862215066960774064332961282599376e-14
    @test zeta(27, 2) ≅ 7.4507117898354294919810041706041194547190318825658299932e-9
    @test 1e-12 > relerr(zeta(27, -105.3), 1.3113726525492708826840989036205762823329453315093955e14)
    @test polygamma(4, -3.1+Inf*im) == polygamma(4, 3.1+Inf*im) == 0
    @test polygamma(4, -0.0) == Inf == -polygamma(4, +0.0)
    @test zeta(4, +0.0) == zeta(4, -0.0) ≅ pi^4 / 90
    @test zeta(5, +0.0) == zeta(5, -0.0) ≅ 1.036927755143369926331365486457034168057080919501912811974
    @test zeta(Inf, 1.) == 1
    @test zeta(Inf, 2.) == 0
    @test isnan(zeta(NaN, 1.))
    @test isa([digamma(x) for x in [1.0]], Vector{Float64})
    @test isa([trigamma(x) for x in [1.0]], Vector{Float64})
    @test isa([polygamma(3,x) for x in [1.0]], Vector{Float64})
    @test zeta(2 + 1im, -1.1) ≅ zeta(2 + 1im, -1.1+0im) ≅ -64.580137707692178058665068045847533319237536295165484548 + 73.992688148809018073371913557697318846844796582012921247im
    @test polygamma(3,5) ≈ polygamma(3,5.)

    @test zeta(-3.0, 7.0) ≅ -52919/120
    @test zeta(-3.0, -7.0) ≅ 94081/120
    @test zeta(-3.1, 7.2) ≅ -587.457736596403704429103489502662574345388906240906317350719
    @test zeta(-3.1, -7.2) ≅ 1042.167459863862249173444363794330893294733001752715542569576
    @test zeta(-3.1, 7.0) ≅ -518.431785723446831868686653718848680989961581500352503093748
    @test zeta(-3.1, -7.0) ≅ 935.1284612957581823462429983411337864448020149908884596048161
    @test zeta(-3.1-0.1im, 7.2) ≅ -579.29752287650299181119859268736614017824853865655709516268 - 96.551907752211554484321948972741033127192063648337407683877im
    @test zeta(-3.1-0.1im, -7.2) ≅ 1025.17607931184231774568797674684390615895201417983173984531 + 185.732454778663400767583204948796029540252923367115805842138im
    @test zeta(-3.1-0.1im, 7.2 + 0.1im) ≅ -571.66133526455569807299410569274606007165253039948889085762 - 131.86744836357808604785415199791875369679879576524477540653im
    @test zeta(-3.1-0.1im, -7.2 + 0.1im) ≅ 1035.35760409421020754141207226034191979220047873089445768189 + 130.905870774271320475424492384335798304480814695778053731045im
    @test zeta(-3.1-0.1im, -7.0 + 0.1im) ≅ 929.546530292101383210555114424269079830017210969572819344670 + 113.646687807533854478778193456684618838875194573742062527301im
    @test zeta(-3.1, 7.2 + 0.1im) ≅ -586.61801005507638387063781112254388285799318636946559637115 - 36.148831292706044180986261734913443701649622026758378669700im
    @test zeta(-3.1, -7.2 + 0.1im) ≅ 1041.04241628770682295952302478199641560576378326778432301623 - 55.7154858634145071137760301929537184886497572057171143541058im
    @test zeta(-13.4, 4.1) ≅ -3.860040842156185186414774125656116135638705768861917e6
    @test zeta(3.2, -4) ≅ 2.317164896026427640718298719837102378116771112128525719078
    @test zeta(3.2, 0) ≅ 1.166773370984467020452550350896512026772734054324169010977
    @test zeta(-3.2+0.1im, 0.0) ≅ zeta(-3.2+0.1im, 0.0+0im) ≅ 0.0070547946138977701155565365569392198424378109226519905493 + 0.00076891821792430587745535285452496914239014050562476729610im
    @test zeta(-3.2, 0.0) ≅ zeta(-3.2, 0.0+0im) ≅ 0.007011972077091051091698102914884052994997144629191121056378

    @test 1e-14 > relerr(eta(1+1e-9), 0.693147180719814213126976796937244130533478392539154928250926)
    @test 1e-14 > relerr(eta(1+5e-3), 0.693945708117842473436705502427198307157819636785324430166786)
    @test 1e-13 > relerr(eta(1+7.1e-3), 0.694280602623782381522315484518617968911346216413679911124758)
    @test 1e-13 > relerr(eta(1+8.1e-3), 0.694439974969407464789106040237272613286958025383030083792151)
    @test 1e-13 > relerr(eta(1 - 2.1e-3 + 2e-3 * im), 0.69281144248566007063525513903467244218447562492555491581+0.00032001240133205689782368277733081683574922990400416791019im)
    @test 1e-13 > relerr(eta(1 + 5e-3 + 5e-3 * im), 0.69394652468453741050544512825906295778565788963009705146+0.00079771059614865948716292388790427833787298296229354721960im)
    @test 1e-12 > relerr(zeta(1e-3+1e-3im), -0.5009189365276307665899456585255302329444338284981610162-0.0009209468912269622649423786878087494828441941303691216750im)
    @test 1e-13 > relerr(zeta(1e-4 + 2e-4im), -0.5000918637469642920007659467492165281457662206388959645-0.0001838278317660822408234942825686513084009527096442173056im)

    # Issue #7169:
    @test 1e-13  > relerr(zeta(0 + 99.69im), 4.67192766128949471267133846066040655597942700322077493021802+3.89448062985266025394674304029984849370377607524207984092848im)
    @test 1e-12 > relerr(zeta(3 + 99.69im), 1.09996958148566565003471336713642736202442134876588828500-0.00948220959478852115901654819402390826992494044787958181148im)
    @test 1e-13  > relerr(zeta(-3 + 99.69im), 10332.6267578711852982128675093428012860119184786399673520976+13212.8641740351391796168658602382583730208014957452167440726im)
    @test 1e-13 > relerr(zeta(2 + 99.69im, 1.3), 0.41617652544777996034143623540420694985469543821307918291931-0.74199610821536326325073784018327392143031681111201859489991im)

    # issue #128
    @test 1e-13 > relerr(zeta(.4 + 453.0im), 5.595631794716693 - 4.994584420588448im)
    @test 1e-10 > relerr(zeta(.4 + 4053.0im), -0.1248993234383550+0.9195498409364987im)
    @test 1e-13 > relerr(zeta(.4 + 12.01im), 1.0233184799021265846512208845-0.8008078492939259287905322251im)
    @test zeta(.4 + 12.01im) == conj(zeta(.4 - 12.01im))
end

@testset "logabsbinomial" begin
    @test logabsbinomial(10, -1) == (-Inf, 0.0)
    @test logabsbinomial(10, 11) == (-Inf, 0.0)
    @test logabsbinomial(10,  0) == ( 0.0, 1.0)
    @test logabsbinomial(10, 10) == ( 0.0, 1.0)

    @test logabsbinomial(10,  1)[1]   ≈ log(10)
    @test logabsbinomial(10,  1)[2]   == 1.0
    @test logabsbinomial(-6, 10)[1]   ≈ log(binomial(-6, 10))
    @test logabsbinomial(-6, 10)[2]   == 1.0
    @test logabsbinomial(-6, 11)[1]   ≈ log(abs(binomial(-6, 11)))
    @test logabsbinomial(-6, 11)[2]   == -1.0
    @test first.(logabsbinomial.(200, 0:200)) ≈ log.(binomial.(BigInt(200), (0:200)))
end

@testset "beta, lbeta" begin
    @test beta(3/2,7/2) ≈ 5π/128
    @test beta(3,5) ≈ 1/105
    @test logbeta(5,4) ≈ log(beta(5,4))
    @test beta(5,4) ≈ beta(4,5)
    @test beta(-2.0,2.0) ≈ 0.5
    @test logabsbeta(-2.0,2.0)[1] ≈ -0.69314718055994529
    @test beta(1e8,0.5) ≈ 0.00017724538531210809
    @test logabsbeta(1e8,0.5)[1] ≈ -8.637975427801484
    @test beta(-1/2, 3) ≈ beta(-1/2 + 0im, 3 + 0im) ≈ -16/3
    @test logabsbeta(-1/2, 3)[1] ≈ log(16/3)
    @test beta(Float32(5),Float32(4)) == beta(Float32(4),Float32(5))
    @test beta(3,5) ≈ beta(3+0im,5+0im)
    @test(beta(3.2+0.1im,5.3+0.3im) ≈ exp(logbeta(3.2+0.1im,5.3+0.3im)) ≈
          0.00634645247782269506319336871208405439180447035257028310080 -
          0.00169495384841964531409376316336552555952269360134349446910im)

    @test beta(big(1.0),big(1.2)) ≈ beta(1.0,1.2) rtol=4*eps()
end


text:
module BackendBase
export AbstractBigArrayBackend, get_info, get_scale_name
abstract type AbstractBigArrayBackend <: AbstractDict{String, Any} end

function get_info end
function get_scale_name end

end # module