text:
const INT_PER_PAGE = div(ccall(:jl_getpagesize, Clong, ()), sizeof(Int))
# we use this to mark pages not in use, it must never be written to.
const NULL_INT_PAGE = Vector{Int}()

mutable struct SparseIntSet
    packed ::Vector{Int}
    reverse::Vector{Vector{Int}}
    counters::Vector{Int}  # counts the number of real elements in each page of reverse.
end

SparseIntSet() = SparseIntSet(Int[], Vector{Int}[], Int[])

SparseIntSet(indices) = union!(SparseIntSet(), indices)

eltype(::Type{SparseIntSet}) = Int

empty(::SparseIntSet) = SparseIntSet()

function empty!(s::SparseIntSet)
    empty!(s.packed)
    empty!(s.reverse)
    empty!(s.counters)
    return s
end

isempty(s::SparseIntSet) = isempty(s.packed)

copy(s::SparseIntSet) = copy!(SparseIntSet(), s)

function copy!(to::SparseIntSet, from::SparseIntSet)
    to.packed = copy(from.packed)
    #we want to keep the null pages === NULL_INT_PAGE
    resize!(to.reverse, length(from.reverse))
    for i in eachindex(from.reverse)
        page = from.reverse[i]
        if page === NULL_INT_PAGE
            to.reverse[i] = NULL_INT_PAGE
        else
            to.reverse[i] = copy(from.reverse[i])
        end
    end
    to.counters = copy(from.counters)
    return to
end

function pageid_offset(s::SparseIntSet, i)
    pageid = div(i - 1, INT_PER_PAGE) + 1
    return pageid, (i - 1) & (INT_PER_PAGE - 1) + 1
end

function in(i, s::SparseIntSet)
    pageid, offset = pageid_offset(s, i)
    if pageid > length(s.reverse)
        return false
    else
        page = @inbounds s.reverse[pageid]
        return page !== NULL_INT_PAGE &&  @inbounds page[offset] != 0
    end
end

length(s::SparseIntSet) = length(s.packed)

@inline function push!(s::SparseIntSet, i::Integer)
    i <= 0 && throw(DomainError("Only positive Ints allowed."))

    pageid, offset = pageid_offset(s, i)
    pages = s.reverse
    plen = length(pages)

    if pageid > plen
        # Create new null pages up to pageid and fresh (zero-filled) one at pageid
        sizehint!(pages, pageid)
        sizehint!(s.counters, pageid)
        for i in 1:pageid - plen - 1
            push!(pages, NULL_INT_PAGE)
            push!(s.counters, 0)
        end
        push!(pages, zeros(Int, INT_PER_PAGE))
        push!(s.counters, 0)
    elseif pages[pageid] === NULL_INT_PAGE
        #assign a page to previous null page
        pages[pageid] = zeros(Int, INT_PER_PAGE)
    end
    page = pages[pageid]
    if page[offset] == 0
        @inbounds page[offset] = length(s) + 1
        @inbounds s.counters[pageid] += 1
        push!(s.packed, i)
        return s
    end
    return s
end

@inline function push!(s::SparseIntSet, is::Integer...)
    for i in is
        push!(s, i)
    end
    return s
end

@inline Base.@propagate_inbounds function pop!(s::SparseIntSet)
    if isempty(s)
        throw(ArgumentError("Cannot pop an empty set."))
    end
    id = pop!(s.packed)
    pageid, offset = pageid_offset(s, id)
    @inbounds s.reverse[pageid][offset] = 0
    @inbounds s.counters[pageid] -= 1
    cleanup!(s, pageid)
    return id
end

@inline Base.@propagate_inbounds function pop!(s::SparseIntSet, id::Integer)
    id < 0 && throw(ArgumentError("Int to pop needs to be positive."))

    @boundscheck if !in(id, s)
        throw(BoundsError(s, id))
    end
    @inbounds begin
        packed_endid = s.packed[end]
        from_page, from_offset = pageid_offset(s, id)
        to_page, to_offset = pageid_offset(s, packed_endid)

        packed_id = s.reverse[from_page][from_offset]
        s.packed[packed_id] = packed_endid
        s.reverse[to_page][to_offset] = s.reverse[from_page][from_offset]
        s.reverse[from_page][from_offset] = 0
        s.counters[from_page] -= 1
        pop!(s.packed)
    end
    cleanup!(s, from_page)
    return id
end

@inline function cleanup!(s::SparseIntSet, pageid::Int)
    if s.counters[pageid] == 0
        s.reverse[pageid] = NULL_INT_PAGE
    end
end

@inline function pop!(s::SparseIntSet, id::Integer, default)
    id < 0 && throw(ArgumentError("Int to pop needs to be positive."))
    return in(id, s) ? (@inbounds pop!(s, id)) : default
end
popfirst!(s::SparseIntSet) = pop!(s, first(s))

@inline iterate(set::SparseIntSet, args...) = iterate(set.packed, args...)

last(s::SparseIntSet) = isempty(s) ? throw(ArgumentError("Empty set has no last element.")) : last(s.packed)

union(s::SparseIntSet, ns) = union!(copy(s), ns)
function union!(s::SparseIntSet, ns)
    for n in ns
        push!(s, n)
    end
    return s
end

intersect(s1::SparseIntSet) = copy(s1)
intersect(s1::SparseIntSet, ss...) = intersect(s1, intersect(ss...))
function intersect(s1::SparseIntSet, ns)
    s = SparseIntSet()
    for n in ns
        n in s1 && push!(s, n)
    end
    return s
end

intersect!(s1::SparseIntSet, ss...) = intersect!(s1, intersect(ss...))

#Is there a more performant way to do this?
intersect!(s1::SparseIntSet, ns) = copy!(s1, intersect(s1, ns))

setdiff(s::SparseIntSet, ns) = setdiff!(copy(s), ns)
function setdiff!(s::SparseIntSet, ns)
    for n in ns
        pop!(s, n, nothing)
    end
    return s
end

function ==(s1::SparseIntSet, s2::SparseIntSet)
    length(s1) != length(s2) && return false
    return all(in(s1), s2)
end

issubset(a::SparseIntSet, b::SparseIntSet) = isequal(a, intersect(a, b))

<(a::SparseIntSet, b::SparseIntSet) = ( a<=b ) && !isequal(a, b)
<=(a::SparseIntSet, b::SparseIntSet) = issubset(a, b)

function findfirst_packed_id(i, s::SparseIntSet)
    pageid, offset = pageid_offset(s, i)
    if pageid > length(s.counters) || s.counters[pageid] == 0
        return 0
    end
    @inbounds id = s.reverse[pageid][offset]
    return id
end

collect(s::SparseIntSet) = copy(s.packed)

struct ZippedSparseIntSetIterator{VT,IT}
    valid_sets::VT
    shortest_set::SparseIntSet
    excluded_sets::IT
    function ZippedSparseIntSetIterator(valid_sets::SparseIntSet...; exclude::NTuple{N, SparseIntSet}=()) where{N}
        shortest = valid_sets[findmin(map(length, valid_sets))[2]]
        new{typeof(valid_sets), NTuple{N, SparseIntSet}}(valid_sets, shortest, exclude)
    end
end

Base.zip(s::SparseIntSet...;kwargs...) = ZippedSparseIntSetIterator(s...;kwargs...)

length(it::ZippedSparseIntSetIterator) = length(it.shortest_set)

# we know it is not in_excluded, as there are no excluded
@inline in_excluded(id, it::ZippedSparseIntSetIterator{VT,Tuple{}}) where {VT} = false

@inline function in_excluded(id, it)
    for e in it.excluded_sets
        if id in e
            return true
        end
    end
    return false
end

@inline function in_valid(id, it)
    for e in it.valid_sets
        if !in(id, e)
            return false
        end
    end
    return true
end

@inline function iterate(it::ZippedSparseIntSetIterator, state=1)
    iterator_length = length(it)
    if state > iterator_length
        return nothing
    end
    for i in state:iterator_length
        @inbounds id = it.shortest_set.packed[i]
        if in_valid(id, it) && !in_excluded(id, it)
            return map(x->findfirst_packed_id(id, x), it.valid_sets), i + 1
        end
    end
    return nothing
end


text:
import Base.LinAlg.BlasInt
using Base.LinAlg.BLAS

@testset "----- Testing Lapack -----" begin

  A = [1.0 2 3; 4 5 6; 8 8 9]
  A2 = copy(A)
  ipiv = zeros(BlasInt, 3)
  info = getrf!(A, ipiv)
  @test ( info )== 0
  b = [1.0, 2, 3]
  b2 = copy(b)
  info = getrs2!('N', A, ipiv, b)
  @test ( info )== 0

  x = copy(b)
  x2 = A2\b2
  @test isapprox( norm(x - x2), 0.0) atol=1e-13

  # test getrs2!
  # test A*x against P*L*U*x

  N = 30
  for i=1:10
    A = rand(N, N)
    x = rand(N)
    b = A*x
    b2 = zeros(b)
    ipiv = zeros(BlasInt, N)
    getrf!(A,ipiv)
    copy!(b2, x)
    trmv!('U', 'N', 'N', A, b2)
    trmv!('L', 'N', 'U', A, b2)
    applyIpiv!(ipiv, b2)

    @test isapprox( norm(b2 - b), 0.0) atol=1e-12
  end

  # test suitesparse solve
  A = [1.0 2 3; 4 5 6; 8 8 9]
  As = sparse(A)
  b = Float64[1, 2, 3]

  x = A\b
  lu = factorize(As)
  x2 = zeros(x)
  solve_suitesparse(lu, b, UMFPACK_A, x2)

  @test isapprox( norm(x2 - x), 0.0) atol=1e-13

end

ext: jl


text:
using LightGraphs
using SimpleWeightedGraphs
using Test

testdir = dirname(@__FILE__)

testgraphs(g) = [g, SimpleWeightedGraph{UInt8,Float64}(g), SimpleWeightedGraph{Int16,Float32}(g)]
testdigraphs(g) = [g, SimpleWeightedDiGraph{UInt8,Float64}(g), SimpleWeightedDiGraph{Int16,Float32}(g)]

testsimplegraphs(g) = [g, LightGraphs.SimpleGraph{UInt8}(g), LightGraphs.SimpleGraph{Int16}(g)]
testsimpledigraphs(g) = [g, LightGraphs.SimpleDiGraph{UInt8}(g), LightGraphs.SimpleDiGraph{Int16}(g)]

tests = [
    "simpleweightededge",
    "simpleweightedgraph",
    "overrides",
    "persistence",
    "connectivity"
]

@testset "SimpleWeightedGraphs" begin
    for t in tests
        tp = joinpath(testdir, "$(t).jl")
        include(tp)
    end
end


text:
# Treasury uses Act/Act (ICMA) day count convention
# unclear how this translates to constant maturities
# I use months = 30 days, years = 365 days
const DAYSTOMATURITY = [30,60,90,180,365,730,1095,1825,2555,3650,7300,10950]
const DAYSTOMATURITYREAL = [1825,2555,3650,7300,10950]

const TTM = (;Iterators.zip(DailyTreasuryYieldCurve.COLNAMES[3:end-1],DAYSTOMATURITY)...)
const TTMREAL = (;Iterators.zip(DailyTreasuryYieldCurve.COLNAMESREAL[3:end],DAYSTOMATURITYREAL)...)

"""
    RateInterpolator

A container for the series of daily yield curve interpolators.

Fields: `dates` and `interpolators`, both `Vector`.

To interpolate, just call it with the value to interpolate (days to maturity) and the date of the curve.
```
(ri::RateInterpolator)(d::Real,dt::Date)
```

Currently supports only linear interpolation and extrapolation. See [`createRateInterpolator`](@ref) for construction.
"""
struct RateInterpolator{T}
    dates::Vector{Date}
    interpolators::Vector{T}
    realrates::Bool
end

function (ri::RateInterpolator)(d::Real,dt::Date)
    idx = searchsortedlast(ri.dates,dt)
    return ri.interpolators[idx](d)
end

"""
    createRateInterpolator(df;realrates::Bool=false)

Preferred method to construct a [`RateInterpolator`](@ref), just pass the `df` that you get from [`getyieldcurves`](@ref).
"""
function createRateInterpolator(df0;realrates::Bool=false)
    loc = [!all(ismissing(x) for x in r) for r in eachrow(select(df0,Not(:date)))] # 0 for rows with only missing data
    df = df0[loc,:]
    dates = unique(df.date)

    dfs = DataFrames.stack(df,Not(:date); variable_eltype=Symbol)
    if realrates
        dfs.DTM = [TTMREAL[k] for k in dfs.variable]
    else
        dfs.DTM = [TTM[k] for k in dfs.variable]
    end

    gd = groupby(dfs,:date)
    rateinterps = [buildsingleinterpolator(g) for g in gd]

    @assert length(dates) == length(rateinterps) "Number of dates is different than number of daily interpolators."

    return RateInterpolator(dates,rateinterps,realrates)
end

"""
    buildsingleinterpolator(df::AbstractDataFrame)

Builds a single (one day) interpolation of the yield curve. Used in `createRateInterpolator`.
"""
function buildsingleinterpolator(df::AbstractDataFrame)
    dfdm = dropmissing(df)
    x = dfdm.DTM
    y = dfdm.value
    itp = LinearInterpolation(x,y,extrapolation_bc=Line())
    return itp
end

function Base.show(io::IO, ri::RateInterpolator)
    itpT = Interpolations.itptype(first(ri.interpolators))
    begdt = first(ri.dates)
    enddt = last(ri.dates)
    ndays = length(ri.dates)
    kind = ri.realrates ? "Real" : "Nominal"
    print(io,"$kind RateInterpolator{$itpT interpolators} starting $begdt and ending $enddt ($ndays days)")
end


text:
setprecision(100)
T = BigFloat
# T = Float64

PI = T(π)
πi = PI * 1im
πi2 = 2πi

using QuadGK
using OffsetArrays
using SpecialFunctions

roots = [
     14.134725142
     21.022039639
     25.010857580
     30.424876126
     32.935061588
     37.586178159
     40.918719012
     43.327073281
     48.005150881
     49.773832478
     52.970321478
     56.446247697
     59.347044003
     60.831778525
     65.112544048
     67.079810529
     69.546401711
     72.067157674
     75.704690699
     77.144840069
     79.337375020
     82.910380854
     84.735492981
     87.425274613
     88.809111208
     92.491899271
     94.651344041
     95.870634228
     98.831194218
    101.317851006
    103.725538040
    105.446623052
    107.168611184
    111.029535543
    111.874659177
    114.320220915
    116.226680321
    118.790782866
    121.370125002
    122.946829294
    124.256818554
    127.516683880
    129.578704200
    131.087688531
    133.497737203
    134.756509753
    138.116042055
    139.736208952
    141.123707404
    143.111845808
    146.000982487
    147.422765343
    150.053520421
    150.925257612
    153.024693811
    156.112909294
    157.597591818
    158.849988171
    161.188964138
    163.030709687
    165.537069188
    167.184439978
    169.094515416
    169.911976479
    173.411536520
    174.754191523
    176.441434298
    178.377407776
    179.916484020
    182.207078484
    184.874467848
    185.598783678
    187.228922584
    189.416158656
    192.026656361
    193.079726604
    195.265396680
    196.876481841
    198.015309676
    201.264751944
    202.493594514
    204.189671803
    205.394697202
    207.906258888
    209.576509717
    211.690862595
    213.347919360
    214.547044783
    216.169538508
    219.067596349
    220.714918839
    221.430705555
    224.007000255
    224.983324670
    227.421444280
    229.337413306
    231.250188700
    231.987235253
    233.693404179
    236.524229666
    237.769820481
    239.555477573
    241.049157796
    242.823271934
    244.070898497
    247.136990075
    248.101990060
    249.573689645
    251.014947795
    253.069986748
    255.306256455
    256.380713694
    258.610439492
    259.874406990
    260.805084505
    263.573893905
    265.557851839
    266.614973782
    267.921915083
    269.970449024
    271.494055642
    273.459609188
    275.587492649
    276.452049503
    278.250743530
    279.229250928
    282.465114765
    283.211185733
    284.835963981
    286.667445363
    287.911920501
    289.579854929
    291.846291329
    293.558434139
    294.965369619
    295.573254879
    297.979277062
    299.840326054
    301.649325462
    302.696749590
    304.864371341
    305.728912602
    307.219496128
    310.109463147
    311.165141530
    312.427801181
    313.985285731
    315.475616089
    317.734805942
    318.853104256
    321.160134309
    322.144558672
    323.466969558
    324.862866052
    327.443901262
    329.033071680
    329.953239728
    331.474467583
    333.645378525
    334.211354833
    336.841850428
    338.339992851
    339.858216725
    341.042261111
    342.054877510
    344.661702940
    346.347870566
    347.272677584
    349.316260871
    350.408419349
    351.878649025
    353.488900489
    356.017574977
    357.151302252
    357.952685102
    359.743754953
    361.289361696
    363.331330579
    364.736024114
    366.212710288
    367.993575482
    368.968438096
    370.050919212
    373.061928372
    373.864873911
    375.825912767
    376.324092231
    378.436680250
    379.872975347
    381.484468617
    383.443529450
    384.956116815
    385.861300846
    387.222890222
    388.846128354
    391.456083564
    392.245083340
    393.427743844
    395.582870011
    396.381854223
    397.918736210
    399.985119876
    401.839228601
    402.861917764
    404.236441800
    405.134387460
    407.581460387
    408.947245502
    410.513869193
    411.972267804
    413.262736070
    415.018809755
    415.455214996
    418.387705790
    419.861364818
    420.643827625
    422.076710059
    423.716579627
    425.069882494
    427.208825084
    428.127914077
    430.328745431
    431.301306931
    432.138641735
    433.889218481
    436.161006433
    437.581698168
    438.621738656
    439.918442214
    441.683199201
    442.904546303
    444.319336278
    446.860622696
    447.441704194
    449.148545685
    450.126945780
    451.403308445
    453.986737807
    454.974683769
    456.328426689
    457.903893064
    459.513415281
    460.087944422
    462.065367275
    464.057286911
    465.671539211
    466.570286931
    467.439046210
    469.536004559
    470.773655478
    472.799174662
    473.835232345
    475.600339369
    476.769015237
    478.075263767
    478.942181535
    481.830339376
    482.834782791
    483.851427212
    485.539148129
    486.528718262
    488.380567090
    489.661761578
    491.398821594
    493.314441582
    493.957997805
    495.358828822
    496.429696216
    498.580782430
    500.309084942
    501.604446965
    502.276270327
    504.499773313
    505.415231742
    506.464152710
    508.800700336
    510.264227944
    511.562289700
    512.623144531
    513.668985555
    515.435057167
    517.589668572
    518.234223148
    520.106310412
    521.525193449
    522.456696178
    523.960530892
    525.077385687
    527.903641601
    528.406213852
    529.806226319
    530.866917884
    532.688183028
    533.779630754
    535.664314076
    537.069759083
    538.428526176
    540.213166376
    540.631390247
    541.847437121
    544.323890101
    545.636833249
    547.010912058
    547.931613364
    549.497567563
    550.970010039
    552.049572201
    553.764972119
    555.792020562
    556.899476407
    557.564659172
    559.316237029
    560.240807497
    562.559207616
    564.160879111
    564.506055938
    566.698787683
    567.731757901
    568.923955180
    570.051114782
    572.419984132
    573.614610527
    575.093886014
    575.807247141
    577.039003472
    579.098834672
    580.136959362
    581.946576266
    583.236088219
    584.561705903
    585.984563205
    586.742771891
    588.139663266
    590.660397517
    591.725858065
    592.571358300
    593.974714682
    595.728153697
    596.362768328
    598.493077346
    599.545640364
    601.602136736
    602.579167886
    603.625618904
    604.616218494
    606.383460422
    608.413217311
    609.389575155
    610.839162938
    611.774209621
    613.599778676
    614.646237872
    615.538563369
    618.112831366
    619.184482598
    620.272893672
    621.709294528
    622.375002740
    624.269900018
    626.019283428
    627.268396851
    628.325862359
    630.473887438
    630.805780927
    632.225141167
    633.546858252
    635.523800311
    637.397193160
    637.925513981
    638.927938267
    640.694794669
    641.945499666
    643.278883781
    644.990578230
    646.348191596
    647.761753004
    648.786400889
    650.197519345
    650.668683891
    653.649571605
    654.301920586
    655.709463022
    656.964084599
    658.175614419
    659.663845973
    660.716732595
    662.296586431
    664.244604652
    665.342763096
    666.515147704
    667.148494895
    668.975848820
    670.323585206
    672.458183584
    673.043578286
    674.355897810
    676.139674364
    677.230180669
    677.800444746
    679.742197883
    681.894991533
    682.602735020
    684.013549814
    684.972629862
    686.163223588
    687.961543185
    689.368941362
    690.474735032
    692.451684416
    693.176970061
    694.533908700
    695.726335921
    696.626069900
    699.132095476
    700.296739132
    701.301742955
    702.227343146
    704.033839296
    705.125813955
    706.184654800
    708.269070885
    709.229588570
    711.130274180
    711.900289914
    712.749383470
    714.082771821
    716.112396454
    717.482569703
    718.742786545
    719.697100988
    721.351162219
    722.277504976
    723.845821045
    724.562613890
    727.056403230
    728.405481589
    728.758749796
    730.416482123
    731.417354919
    732.818052714
    734.789643252
    735.765459209
    737.052928912
    738.580421171
    739.909523674
    740.573807447
    741.757335573
    743.895013142
    745.344989551
    746.499305899
    747.674563624
    748.242754465
    750.655950362
    750.966381067
    752.887621567
    754.322370472
    755.839308976
    756.768248440
    758.101729246
    758.900238225
    760.282366984
    762.700033250
    763.593066173
    764.307522724
    766.087540100
    767.218472156
    768.281461807
    769.693407253
    771.070839314
    772.961617566
    774.117744628
    775.047847097
    775.999711963
    777.299748530
    779.157076949
    780.348925004
    782.137664391
    782.597943946
    784.288822612
    785.739089701
    786.461147451
    787.468463816
    790.059092364
    790.831620468
    792.427707609
    792.888652563
    794.483791870
    795.606596156
    797.263470038
    798.707570166
    799.654336211
    801.604246463
    802.541984878
    803.243096204
    804.762239113
    805.861635667
    808.151814936
    809.197783363
    810.081804886
    811.184358847
    812.771108389
    814.045913608
    814.870539626
    816.727737714
    818.380668866
    819.204642171
    820.721898444
    821.713454133
    822.197757493
    824.526293872
    826.039287377
    826.905810954
    828.340174300
    829.437010968
    830.895884053
    831.799777659
    833.003640909
    834.651915148
    836.693576188
    837.347335060
    838.249021993
    839.465394810
    841.036389829
    842.041354207
    844.166196607
    844.805993976
    846.194769928
    847.971717640
    848.489281181
    849.862274349
    850.645448466
    853.163112583
    854.095511720
    855.286710244
    856.484117491
    857.310740603
    858.904026466
    860.410670896
    861.171098213
    863.189719772
    864.340823930
    865.594664327
    866.423739904
    867.693122612
    868.670494229
    870.846902326
    872.188750822
    873.098978971
    873.908389235
    875.985285109
    876.600825833
    877.654698341
    879.380951970
    880.834648848
    882.386696627
    883.430331839
    884.198743115
    885.272304480
    886.852801963
    888.475566674
    889.735294294
    890.813132113
    892.386433260
    893.119117567
    894.886292321
    895.397919675
    896.632251556
    899.221522668
    899.858884608
    900.849739861
    902.243207587
    903.099674443
    904.702902722
    905.829940758
    907.656729469
    908.333543645
    910.186334057
    911.234951486
    912.331045600
    912.823999247
    914.730096958
    916.355000809
    917.825377570
    918.836535244
    919.448344440
    921.156395507
    922.500629307
    923.285719802
    924.773483933
    926.551552785
    927.850858986
    928.663659329
    929.874092851
    931.009211337
    931.852740746
    934.385306837
    934.995424864
    936.228649379
    937.532925712
    939.024300899
    939.660940615
    941.156999642
    942.052341643
    944.188035810
    945.333562503
    946.765842205
    947.079183096
    948.346646255
    950.151612685
    951.033248734
    952.727988620
    954.129719270
    954.829308938
    956.675479343
    957.510052596
    958.414593390
    959.459168807
    961.669572474
    963.182086671
    963.567040192
    965.055579624
    966.110754818
    967.371153766
    968.636301906
    970.125610557
    971.071491486
    973.185361294
    973.873078993
    974.774635066
    976.178502421
    976.917202117
    978.766671535
    980.578000640
    981.288615302
    982.396485169
    983.575076006
    985.186928656
    986.130515110
    986.756008408
    988.992622371
    990.223917804
    991.374294148
    992.728696337
    993.214580957
    994.404590571
    996.205336164
    997.511934752
    998.827547137
    999.791571557
]

function get_primes(l, r)
    a = ones(Bool, r - l + 1)
    for i in 2:Int64(floor(√r))
        li = max((l - 1) ÷ i + 1, 2)
        ri = r ÷ i
        for j in li : ri
            a[i * j - l + 1] = 0
        end
    end

    ret = []
    for i in max(l, 2):r
        if a[i - l + 1]
            push!(ret, i)
        end
    end
    ret
end

function brentq(f, xa, xb, xtol, rtol, iter)
    xpre, xcur, xblk, fpre, fcur, fblk, spre, scur = xa, xb, 0, f(xa), f(xb), 0, 0, 0

    fpre * fcur > 0 && return 0
    fpre == 0 && return xpre
    fcur == 0 && return xcur

    for _ in 1:iter
        if fpre * fcur < 0
            xblk = xpre
            fblk = fpre
            spre = xcur - xpre
            scur = spre
        end

        if abs(fblk) < abs(fcur)
            xpre = xcur
            xcur = xblk
            xblk = xpre

            fpre = fcur
            fcur = fblk
            fblk = fpre
        end

        Δ = (xtol + rtol * abs(xcur)) / 2
        sbis = (xblk - xcur) / 2
        (fcur == 0 || abs(sbis) < Δ) && return xcur

        if abs(spre) < Δ && abs(fcur) < abs(fpre)
            if xpre == xblk
                stry = -fcur * (xcur - xpre) / (fcur - fpre)
            else
                dpre = (fpre - fcur) / (xpre - xcur);
                dblk = (fblk - fcur) / (xblk - xcur);
                stry = -fcur * (fblk * dblk - fpre * dpre) / (dblk * dpre * (fblk - fpre))
            end

            if 2 * abs(stry) < min(abs(spre), 3 * abs(sbis) - Δ)
                spre = scur
                scur = stry
            else
                spre, scur = sbis, sbis
            end
        else
            spre, scur = sbis, sbis
        end

        xpre = xcur
        fpre = fcur

        if abs(scur) > Δ
            xcur += scur
        elseif sbis > 0
            xcur += Δ
        else
            xcur -= Δ
        end

        fcur = f(xcur)
    end

    xcur
end

function Δ(x, λ, ϵ)
    Φ(ρ) = erfc(ρ / √2) / 2
    ϕ(u, x, λ) = Φ(log(u / x) / λ)
    E_plus(u, x, λ) = x * exp(λ^2 / 2) * Φ(log(u / x) / λ - λ) - u * Φ(log(u / x) / λ)
    E_minus(u, x, λ) = u * Φ(-log(u / x) / λ) - x * exp(λ^2 / 2) * Φ(λ - log(u / x) / λ)

    ϵ = ϵ / 2
    x1 = floor(x)
    if E_minus(x1, x, λ) > ϵ
        g(u) = E_minus(u, x, λ) - 0.9 * ϵ
        u = brentq(g, 2, x, 1e-8, 1e-8, 100)
        x1 = floor(u)
    end

    x2 = ceil(x)
    if E_plus(x2, x, λ) > ϵ
        f(u) = E_plus(u, x, λ) - 0.9 * ϵ
        u = brentq(f, x * 2, x, 1e-8, 1e-8, 100)
        x2 = ceil(u)
    end
    println("[Δ] x1 = $x1, x2 = $x2, diff = $(x2 - x1), est diff = $(√(log(λ * x / ϵ)) * λ * x)")
    println("Eminus $(E_minus(x1, x, λ))")
    println("Eplus $(E_plus(x2, x, λ))")

    N = (x2 - x1 + 1) + ceil(√x2)
    ϵ_1 = ϵ * 2^ceil(-log2(N))

    ret = 0
    for p in get_primes(Int64(x1), Int64(x2))
        ret -= ϕ(T(p), x, λ)
        if p ≤ x
            ret += 1
        end
    end
    for p in get_primes(2, Int64(floor(√x2)))
        m = 2
        while p^m < x2
            if p^m < x1
                ret -= T(1 / m)
            else
                ret -= ϕ(p^m, x, λ) / m
            end
            m += 1
        end
    end

    ret
end

function myint(x, λ, r0, r1)
    ϕ̂(s) = x^s / s * exp(λ^2 * s^2 / 2)
    s0 = 1/2 + r0 * 1im

    N = 30

    # approx exp(-λ^2 h^2)
    c1 = OffsetArray(zeros(Complex{Float64}, N + 1), 0:N)
    c1[0] = 1
    for i in 1:N ÷ 2
        c1[i * 2] = c1[(i - 1) * 2] / -i * λ^2 / 2
    end

    # approx 1/(1+ih/s0)
    c2 = OffsetArray(zeros(Complex{Float64}, N + 1), 0:N)
    c2[0] = 1
    for i in 1:N
        c2[i] = c2[i - 1] * -1im / s0
    end

    coefs = OffsetArray(zeros(Complex{Float64}, N * 2 + 1), 0:2N)
    for i in 0:2:N
        for j in 0:N
            coefs[i + j] += c1[i] * c2[j]
        end
    end

    c = 1im * (s0 * λ^2 + log(x))
    println(c, " ", s0, λ)
    dh = r1 - r0
    # redt = ϕ̂(s0) * \sum_i=0^{2N} \int_0^{dh} exp(c*h) h^i d h
    factorial = OffsetArray(zeros(0:2N), 0:2N)
    factorial[0] = 1
    for i in 1:2N
        factorial[i] = factorial[i - 1] * i
    end

    integral = 0
    for n in 0:2N
        w = 0
        for i in 0:n
            w += (-1)^(n - i) * factorial[n] / factorial[i] / c^(n - i + 1) * dh^i
        end
        w *= exp(c * dh)
        w -= (-1)^n * factorial[n] / c^(n + 1)

        quad = quadgk((h) -> imag(exp(c * h) * h^n), 0, dh)
#         println("n = $n, w = $w, quad = $quad")
        integral += w * coefs[n]
    end

    base = ϕ̂(s0)

    ret = imag(base * integral)
    quad = quadgk((h) -> imag(ϕ̂(1/2 + h * 1im)), r0, r1)
#     println("ret = $ret, quad = $quad")

    return ret
end

function platt(x)
    err_remain_int(λ) = exp(λ^2 / 2) / 2π / x / λ * (5π * √2 + 2 / λ)

    ϵ = 0.5
#     λ = brentq((λ) -> err_remain_int(λ) - ϵ * 0.4, 0, 2, 1e-9, 1e-9, 100)
    λ = 0.01
    est_λ = √(1 / π / x / (ϵ * 0.4))

    println("λ = $λ, err_remain_int = $(err_remain_int(λ)), est λ = $est_λ")
    @assert err_remain_int(λ) < ϵ / 2

    ϕ̂(s) = x^s / s * exp(λ^2 * s^2 / 2)
    B(σ, T) = exp(λ^2 * (1 - T^2) / 2) * (x^σ / T / log(x) + 1 / λ^2 / T^2 / x)
    ans = 0
    integral = quadgk((h) -> imag(ϕ̂(1 + h * 1im)), 0, +Inf)[1] - log(2)

    whole_int = 0

    for i in 1:length(roots) - 1
        r0 = roots[i]
        r1 = roots[i + 1]
        err = B(1/2, r0)
        if err < 1 / Float64(x)^3
            break
        end
        v = myint(x, λ, r0, r1)
        whole_int += v
        integral -= v * i * 2  # count two because two symmetric roots
        println("r = $r0, v = $v, err = $err")
    end
    Delta = Δ(x, λ, ϵ * 0.4)
    ans = integral + Delta
    println("integral = $integral, Delta = $Delta, whole_int = $whole_int, fixed = $(quadgk((h) -> imag(ϕ̂(1/2 + h * 1im)), roots[1], +Inf)[1])")
    println(ans)
end

platt(100000)

max_stars_repo_name: roosephu/apc

max_issues_repo_name: roosephu/apc

max_forks_repo_name: roosephu/apc


text:
using Genie, Genie.Renderer.Html, Stipple, StipplePlotly

Genie.config.log_requests = false

xrange = 0.0:(2π/1000):2π

dxmeasured = 0.05
dymeasured = 0.1

xexperiment = (0.0:(2π/10):2π) .+ 3 .* dxmeasured .* (rand(Float64,11) .- 0.5)
dx = dxmeasured .* ones(Float64, size(xexperiment))
yexperiment = sin.(xexperiment) .+ 3 .* dymeasured .* (rand(Float64,size(xexperiment)) .- 0.5)
dy = dymeasured .* ones(Float64, size(yexperiment))
# create outlier:
yexperiment[6] = 0.5

pd_line(name, xar) = PlotData(
  x = xar,
  y = sin.(xar),
  plot = StipplePlotly.Charts.PLOT_TYPE_SCATTER,
  mode = "lines",
  name = name
)

pd_scatter(name, xar, dx, yar, dy) = PlotData(
  x = xar,
  error_x = ErrorBar(dx),
  y = yar,
  error_y = ErrorBar(dy),
  plot = StipplePlotly.Charts.PLOT_TYPE_SCATTER,
  mode = "markers",
  name = name
)

pl() = PlotLayout(
  plot_bgcolor = "#FFFFFF",
  title = PlotLayoutTitle(text="Wave", font=Font(24)),
  legend = PlotLayoutLegend(bgcolor = "rgb(212,212,212)", font=Font(6)),
  hovermode = "closest",
  showlegend = true,
  xaxis = [PlotLayoutAxis(xy="x", index=1,
    title = "time (s)",
    ticks = "outside",
    tickfont = Font(size=24, color="#FF00FF"),
    showline = true,
    zeroline = false,
    mirror = true
  )],
  yaxis = [PlotLayoutAxis(xy="y", index=1,
    showline = true,
    zeroline = false,
    title = "displacement (mm)",
    ticks = "outside",
    mirror = true
  )],
  annotations = [PlotAnnotation(visible=true, x=xexperiment[6], y=yexperiment[6], text="possible outlier")]
)

Base.@kwdef mutable struct Model <: ReactiveModel
  data::R{Vector{PlotData}} = [pd_line("Sinus", xrange), pd_scatter("Experiment", xexperiment, dx, yexperiment, dy)]
  layout::R{PlotLayout} = pl()
  config::R{PlotConfig} = PlotConfig()
end

model = Stipple.init(Model())

function ui()
  page(
    vm(model), class="container", [
        plot(:data, layout = :layout, config = :config)
    ]
  ) |> html
end

route("/", ui)

up()


text:
#### Solution types

export CompactSolution, VerboseSolution

abstract type SolutionType end

"""
    VerboseSolution <: SolutionType

Used to return a [`VerboseSolutionResults`](@ref)

!!! note
    VerboseSolution is designed for debugging purposes
    only, and is not GPU-friendly.
"""
struct VerboseSolution <: SolutionType end

abstract type AbstractSolutionResults{AbstractFloat} end

"""
    VerboseSolutionResults{AT} <: AbstractSolutionResults{AT}

Result returned from `find_zero` when
`VerboseSolution` is passed as the `soltype`.
"""
struct VerboseSolutionResults{AT} <: AbstractSolutionResults{AT}
    "solution ``x^*`` of the root of the equation ``f(x^*) = 0``"
    root::AT
    "indicates convergence"
    converged::Bool
    "error of the root of the equation ``f(x^*) = 0``"
    err::AT
    "number of iterations performed"
    iter_performed::Int
    "solution per iteration"
    root_history::Vector{AT}
    "error of the root of the equation ``f(x^*) = 0`` per iteration"
    err_history::Vector{AT}
end
SolutionResults(soltype::VerboseSolution, args...) =
    VerboseSolutionResults(args...)

"""
    CompactSolution <: SolutionType

Used to return a [`CompactSolutionResults`](@ref)
"""
struct CompactSolution <: SolutionType end

"""
    CompactSolutionResults{AT} <: AbstractSolutionResults{AT}

Result returned from `find_zero` when
`CompactSolution` is passed as the `soltype`.
"""
struct CompactSolutionResults{AT} <: AbstractSolutionResults{AT}
    "solution ``x^*`` of the root of the equation ``f(x^*) = 0``"
    root::AT
    "indicates convergence"
    converged::Bool
end
SolutionResults(soltype::CompactSolution, root, converged, args...) =
    CompactSolutionResults(root, converged)

init_history(::VerboseSolution, ::Type{AT}) where {AT <: AbstractArray} = AT[]
init_history(::CompactSolution, ::Type{AT}) where {AT <: AbstractArray} =
    nothing

function push_history!(
    history::Vector{AT},
    x::AT,
    ::VerboseSolution,
) where {AT <: AbstractArray}
    push!(history, x)
end
function push_history!(
    history::Nothing,
    x::AT,
    ::CompactSolution,
) where {AT <: AbstractArray}
    nothing
end


text:
#path = "C:\\Users\\wilhe\\Dropbox\\My PC (DESKTOP-P6322LG)\\Desktop\\GELU4_ANN.csv"
#path = "C:\\Users\\wilhe\\Dropbox\\My PC (DESKTOP-P6322LG)\\Desktop\\TANH2_ANN.csv"
path = "C:\\Users\\wilhe\\Dropbox\\My PC (DESKTOP-P6322LG)\\Desktop\\SILU2_ANN.csv"
nn_data_raw = CSV.read(path, DataFrame)
size_df = size(nn_data_raw)

input_len = 2
hidden_len = 6
output_len = 1
depth = 2
g = swish1

# extracts vector of W and b that
W = Matrix{Float64}[]
b = Vector{Float64}[]

# extract first hidden layer info
H1frame = filter(row -> row.Layer === "H1", nn_data_raw)
H1frame_W = H1frame[:, [Symbol("W"*string(i)) for i=1:input_len]]
H1frame_b = H1frame[:, [:b]]
Welement = Matrix{Float64}(H1frame_W)
belement = Matrix{Float64}(H1frame_b)
push!(W, Welement)
push!(b, Float64[belement[i] for i = 1:size(belement,1)])

# extract info for 2 to last hidden layer
for i = 2:depth
    Hiframe = filter(row -> row.Layer === "H$i", nn_data_raw)
    Hiframe_W = Hiframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
    Hiframe_b = Hiframe[:, [:b]]
    Welement2 = Matrix{Float64}(Hiframe_W)
    belement2 = Matrix{Float64}(Hiframe_b)
    push!(W, Welement2)
    push!(b, Float64[belement2[i] for i = 1:size(belement2,1)])
end

# extract info for output layer
Hoframe = filter(row -> row.Layer === "Output", nn_data_raw)
Hoframe_W = Hoframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
Hoframe_b = Hoframe[:, [:b]]
Welement = Matrix{Float64}(Hoframe_W)
belement = Matrix{Float64}(Hoframe_b)
push!(W, Welement)
push!(b, Float64[belement[i] for i = 1:size(belement,1)])

#mutable struct ANN1 <: Function
#end

xmin = zeros(2)
xmax = [40.0; 9.1]
ymin = 0.0
ymax = 0.638966464124928

#=
function ANN1(x)
    xs = (x .- xmin)./(xmax .- xmin)

    yi1 = W[1]*xs .+ b[1]
    y1 = swish1.(yi1)
    #@show yi1
    #@show y1

    yi2 = W[2]*y1 .+ b[2]
    y2 = swish1.(yi2)

    yi3 = W[3]*y2 .+ b[3]
    y3 = swish1.(yi3)

    yi4 = W[4]*y3 .+ b[4]
    y4 = swish1.(yi4)

    yi5 = W[5]*y4 .+ b[5]
    y5 = sigmoid.(yi5)

    y5[1]*(ymax - ymin) + ymin
end
=#

function ANN1(x)

    #@show x
    xs = (x .- xmin)./(xmax .- xmin)

    yi1 = W[1]*xs .+ b[1]
    y1 = swish1.(yi1)
    #@show yi1
    #@show y1

    yi2 = W[2]*y1 .+ b[2]
    y2 = swish1.(yi2)
    #@show yi2
    #@show y2

    yi3 = W[3]*y2 .+ b[3]
    y3 = sigmoid.(yi3)
    #@show yi4
    #@show y4

    #@show y4[1]*(ymax - ymin) + ymin
    return y3[1]*(ymax - ymin) + ymin
end

#=
# builds a nonallocating multilayer perceptron (vanilla ANN) with static storage types
calc_block = quote end

# used to revert minmax scaling
xmin = zeros(2)
xmax = [40.0; 9.1]

# unpacks input variables (CHECKED)
for i = 1:input_len
    sym = Symbol("Input_"*string(i))
    #push!(calc_block.args, :($sym = (x[$i] - $(xmin[i]))/($(xmax[i]) - $(xmin[i]))))
    push!(calc_block.args, :($sym = x[$i]))
end

# input to first hidden layer (CHECKED)
for i = 1:hidden_len
    sym = Symbol("H1_"*string(i))
    sum_expr = :($(b[1][i]))
    for j = 1:input_len
        prior_sym =  Symbol("Input_"*string(j))
        sum_expr = :($sum_expr + $(W[1][i,j])*$prior_sym)
    end
    push!(calc_block.args, :($sym = f($sum_expr)))
end

# 2 to last hidden_layers
for k = 2:depth
    @show k
    for i = 1:hidden_len
        sym = Symbol("H$(k)_"*string(i))
        sum_expr = :($(b[1][i]))
        for j = 1:input_len
            prior_sym =  Symbol("H$(k-1)_"*string(j))
            sum_expr = :($sum_expr + $(W[k][i,j])*$prior_sym)
        end
        push!(calc_block.args, :($sym = f($sum_expr)))
    end
end

output_tuple = :((Output_1,))
for i = 1:output_len
    sym = Symbol("Output_"*string(i))
    if i > 1
        push!(output_tuple, sym)
    end
    sum_expr = :($(b[end][i]))
    for j = 1:hidden_len
        prior_sym = Symbol("H$(depth)_"*string(j))
        sum_expr = :($sum_expr + $(W[end][i,j])*$prior_sym)
    end
    push!(calc_block.args, :($sym = $sum_expr))
end

@generated function ANN(x, f)
    return quote
        $calc_block
        return $output_tuple
    end
end

ANN1(x) = ANN(x, g)[1]

function ANN_GLM1(x)
    ANN1(x) + 0.0020*x[1] + 0.0245*x[2] + 0.4163
end
=#

#=
path = "C:\\Users\\wilhe\\Dropbox\\My PC (DESKTOP-P6322LG)\\Desktop\\GELU1_ANN.csv"
nn_data_raw = CSV.read(path, DataFrame)
size_df = size(nn_data_raw)

input_len = 2
hidden_len = 5
output_len = 1
depth = 1
g = swish1

# extracts vector of W and b that
W = Matrix{Float64}[]
b = Vector{Float64}[]

# extract first hidden layer info
H1frame = filter(row -> row.Layer === "H1", nn_data_raw)
H1frame_W = H1frame[:, [Symbol("W"*string(i)) for i=1:input_len]]
H1frame_b = H1frame[:, [:b]]
Welement = Matrix{Float64}(H1frame_W)
belement = Matrix{Float64}(H1frame_b)
push!(W, Welement)
push!(b, Float64[belement[i] for i = 1:size(belement,1)])

# extract info for 2 to last hidden layer
for i = 2:depth
    Hiframe = filter(row -> row.Layer === "H$i", nn_data_raw)
    Hiframe_W = Hiframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
    Hiframe_b = Hiframe[:, [:b]]
    Welement2 = Matrix{Float64}(Hiframe_W)
    belement2 = Matrix{Float64}(Hiframe_b)
    push!(W, Welement2)
    push!(b, Float64[belement2[i] for i = 1:size(belement2,1)])
end

# extract info for output layer
Hoframe = filter(row -> row.Layer === "Output", nn_data_raw)
Hoframe_W = Hoframe[:, [Symbol("W"*string(i)) for i=1:hidden_len]]
Hoframe_b = Hoframe[:, [:b]]
Welement = Matrix{Float64}(Hoframe_W)
belement = Matrix{Float64}(Hoframe_b)
push!(W, Welement)
push!(b, Float64[belement[i] for i = 1:size(belement,1)])

# builds a nonallocating multilayer perceptron (vanilla ANN) with static storage types
calc_block = quote end

# used to revert minmax scaling
xmin = zeros(2)
xmax = [40.0; 9.1]

# unpacks input variables (CHECKED)
for i = 1:input_len
    sym = Symbol("Input_"*string(i))
    push!(calc_block.args, :($sym = (x[$i] - $(xmin[i]))/($(xmax[i]) - $(xmin[i]))))
end

# input to first hidden layer (CHECKED)
for i = 1:hidden_len
    sym = Symbol("H1_"*string(i))
    sum_expr = :($(b[1][i]))
    for j = 1:input_len
        prior_sym =  Symbol("Input_"*string(j))
        sum_expr = :($sum_expr + $(W[1][i,j])*$prior_sym)
    end
    push!(calc_block.args, :($sym = f($sum_expr)))
end

# 2 to last hidden_layers
for k = 2:depth
    @show k
    for i = 1:hidden_len
        sym = Symbol("H$(k)_"*string(i))
        sum_expr = :($(b[1][i]))
        for j = 1:input_len
            prior_sym =  Symbol("H$(k-1)_"*string(j))
            sum_expr = :($sum_expr + $(W[k][i,j])*$prior_sym)
        end
        push!(calc_block.args, :($sym = f($sum_expr)))
    end
end

output_tuple = :((Output_1,))
for i = 1:output_len
    sym = Symbol("Output_"*string(i))
    if i > 1
        push!(output_tuple, sym)
    end
    sum_expr = :($(b[end][i]))
    for j = 1:hidden_len
        prior_sym = Symbol("H$(depth)_"*string(j))
        sum_expr = :($sum_expr + $(W[end][i,j])*$prior_sym)
    end
    push!(calc_block.args, :($sym = $sum_expr))
end

@generated function ANN(x, f)
    return quote
        $calc_block
        return $output_tuple
    end
end

ANN1(x) = ANN(x, g)[1]

function ANN_GLM1(x)
    ANN1(x) + 0.0020*x[1] + 0.0245*x[2] + 0.4163
end
=#

max_stars_repo_name: PSORLab/RobustHybridModels

max_issues_repo_name: PSORLab/RobustHybridModels

max_forks_repo_name: PSORLab/RobustHybridModels


text:
__precompile__()
module SFGTools

using DataFrames
# using MAT

include("sfspectrum.jl")
include("spectraio.jl")
include("processing.jl")

const N_PIXEL = 512
const VIS_WAVELENGTH = 512.4  # According to service protocol of May 2018

export  SFSpectrum,
        average,
        grab,
        list_spectra,
        load_spectra,
        fieldcorrection!,
        get_attribute,
        get_metadata,
        is_attribute,
        rm_events!,
        get_wavelength,
        get_wavenumber,
        get_ir_wavelength,
        get_ir_wavenumber,
        get_pump_delay,
        dlpos2t,
        get_variables,
        pixelshift,
        pixelshift!,
        findpixelshift,
        savejson,
        save_data,
        failed_sifs,
        check_sif_files,
        get_exposure_time
        # save_mat
end


text:
using Gadfly, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("lattice", "singer"), x=:VoicePart, y=:Height, Geom.boxplot)


text:
using Test, PolyChaos
import LinearAlgebra: norm, dot

##########################################################
# MEASURE
##########################################################
@testset "Error handling for Measure" begin
    @test_throws DomainError JacobiMeasure(-1, 2.)
    @test_throws DomainError JacobiMeasure(-1., -2.)
    @test_throws DomainError JacobiMeasure(1, -2)

    @test_throws DomainError Beta01Measure(0., 2.)
    @test_throws DomainError Beta01Measure(0., -2.)
    @test_throws DomainError Beta01Measure(1, -2)

    @test_throws MethodError Measure("something")

    @test_throws DomainError Measure("test",t->t,(3,-1),true)
end

##########################################################
# ORTHOPOLY
##########################################################
# @test_throws ErrorException OrthoPoly("idonotexist",3)

quads = [ fejer, fejer2, clenshaw_curtis ]
discs = [ lanczos, stieltjes ]
degs  = 1:5:20
m = Uniform01Measure()
tol = 1e-6
@testset "General constructor for Gaussian OrthoPoly" begin
    for quad in quads, disc in discs, deg in degs
        op = Uniform01OrthoPoly(deg)
        op1 = OrthoPoly("myuniform", deg, m; quadrature=quad, discretization = disc, Nquad = 100*deg)
        @test isapprox(norm(coeffs(op)-coeffs(op1),Inf),0; atol = tol)
        op2 = OrthoPoly("myuniform",deg, t->1, (0.,1.), true; quadrature=quad, discretization = disc, Nquad = 100*deg)
        @test isapprox(norm(coeffs(op)-coeffs(op2),Inf),0; atol = tol)
    end
end

##########################################################
# QUAD
##########################################################
Ns = 10:5:30
m = Uniform01Measure()
myf(t) = 6t^5
result = 1.
tol = 1e-8
@testset "Last-resort-constructor for Quad using quadgp" begin
    for quad in quads, N in Ns
        op = OrthoPoly("myOP", N, m)
        NW = nw(Quad(N,m;quadrature=quad))
        @test isapprox(dot(myf.(NW[:,1]),NW[:,2]), result; atol=tol)
        NW = nw(Quad(N, t->1.,(0,1); quadrature=quad))
        @test isapprox(dot(myf.(NW[:,1]),NW[:,2]), result; atol=tol)

        NW = nw(Quad(N,op))
        @test isapprox(dot(myf.(NW[:,1]),NW[:,2]), result; atol=tol)
        # NW = nw(Quad(N, m.w, op.α, op.β, m.dom, m.symmetric))
        # @test isapprox(dot(myf.(NW[:,1]),NW[:,2]), result; atol=tol)
    end
end

d = 5

ops = [ HermiteOrthoPoly(d),
        Uniform01OrthoPoly(d),
        Beta01OrthoPoly(d, 3.2, 5.34) ]
mop = MultiOrthoPoly(ops,d)

is = 0:d
n = length(ops)

x = [1., 2., 3.]

@testset "Evaluation of multivariate basis" begin
    for ind in Iterators.product([collect(0:d) for i = 1:n]...)
        @test isapprox(prod( map(i->evaluate(ind[i],x[i],ops[i]),1:n) ) -  evaluate([ind...],x,mop)[1],0;atol=tol)
    end
end

##########################################################
# TENSOR
##########################################################
deg, M = 5, 1:4
opq = LegendreOrthoPoly(deg; Nrec=3*deg)
mop = MultiOrthoPoly([GaussOrthoPoly(deg; Nrec=3*deg), LogisticOrthoPoly(deg; Nrec=3*deg), Uniform01OrthoPoly(deg; Nrec=3*deg)], deg)

# @test_throws AssertionError Tensor(0,mop)
@testset "Check tensor" begin
    for m in M, op in [opq, mop]
        tensor = Tensor(m,op)
        for ind in Iterators.product([collect(0:deg) for i=1:m]...)
            @test isapprox(tensor.get(collect(ind)), computeSP(collect(ind), op); atol = tol)
        end
    end
end

# test that error is thrown in case a quadrature rule is missing
opq = LegendreOrthoPoly(deg; Nrec=3*deg, addQuadrature=false)
mop = MultiOrthoPoly([GaussOrthoPoly(deg; Nrec=3*deg, addQuadrature=false), LogisticOrthoPoly(deg; Nrec=3*deg), Uniform01OrthoPoly(deg; Nrec=3*deg)], deg)

@test_throws InconsistencyError Tensor(3,opq)
@test_throws InconsistencyError Tensor(3,mop)


text:
# “One size fits all” analytic solutions to the Grad–Shafranov equation
# Phys. Plasmas 17, 032502 (2010); https://doi.org/10.1063/1.3328818
#
# Ideal MHD Freidberg Chapter 6.6.1

@fastmath function _solovev_polynomials(x, y, ::Val{N} = Val(12); logx=log(x)) where N

    x2 = x^2
    x4 = x2^2
    x6 = x2^3
    y2 = y^2
    y4 = y2^2
    y6 = y2^3

    # Eq. 8 & Eq. 27
    ψ₁  = 1
    ψ₂  = x2
    ψ₃  = y2 - x2*logx
    ψ₄  = x4 - 4*x2*y2
    ψ₅  = 2*y4 - 9*y2*x2 + 3*x4*logx - 12*x2*y2*logx
    ψ₆  = x6 - 12*x4*y2 + 8*x2*y4
    ψ₇  = 8*y6 - 140*y4*x2 + 75*y2*x4 - 15*x6*logx + 180*x4*y2*logx - 120*x2*y4*logx

    N == 7 && return @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇]

    x3 = x^3
    y3 = y^3
    y5 = y^5

    ψ₈  = y
    ψ₉  = y*x2
    ψ₁₀ = y3 - 3*y*x2*logx
    ψ₁₁ = 3*y*x4 - 4*y3*x2
    ψ₁₂ = 8*y5 - 45*y*x4 - 80*y3*x2*logx + 60*y*x4*logx

    ψ_polys = @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇, ψ₈, ψ₉, ψ₁₀, ψ₁₁, ψ₁₂]

    return ψ_polys
end

@fastmath function _solovev_polynomials_Dx(x, y, ::Val{N} = Val(12); logx=log(x)) where N

    x2 = x^2
    x3 = x^3
    x4 = x2^2
    x5 = x^5
    x6 = x2^3
    y2 = y^2
    y4 = y2^2

    # Derivatives from Mathematica
    ψ₁  = zero(x)
    ψ₂  = 2*x
    ψ₃  = -(2*x*logx+ x)
    ψ₄  = 4*x3 - 8*x*y2
    ψ₅  = 3*x3 - 30*x*y2 + 12*x3*logx - 24*x*y2*logx
    ψ₆  = 6*x5 - 48*x3*y2 + 16*x*y4
    ψ₇  = -5*x*(3*x4 - 96*x2*y2 + 80*y4 + 6*(3*x4 - 24*x2*y2 + 8*y4)*logx)

    N == 7 && return @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇]

    y3 = y^3

    ψ₈  = zero(x)
    ψ₉  = 2*y*x
    ψ₁₀ = -3*x*y*(1 + 2*logx)
    ψ₁₁ = 12*x3*y - 8*x*y3
    ψ₁₂ = 40*x*y*(-3*x2 - 2*y2 + (6*x2 - 4*y2)*logx)

    ψ_polys = @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇, ψ₈, ψ₉, ψ₁₀, ψ₁₁, ψ₁₂]

    return ψ_polys
end

@fastmath function _solovev_polynomials_Dxx(x, y, ::Val{N} = Val(12); logx=log(x)) where N

    x2 = x^2
    x3 = x^3
    x4 = x2^2
    x6 = x2^3
    y2 = y^2
    y4 = y2^2

    # Derivatives from Mathematica
    ψ₁  = zero(x)
    ψ₂  = 2*one(x)
    ψ₃  = -2*logx - 3
    ψ₄  = 12*x2 - 8*y2
    ψ₅  = 3*(7*x2 - 18*y2 + 4*(3*x2 - 2*y2)*logx)
    ψ₆  = 2*(15*x4 - 72*x2*y2 + 8*y4)
    ψ₇  = -5*(33*x4 - 432*x2*y2 + 128*y4 + 6*(15*x4 - 72*x2*y2 + 8*y4)*logx)

    N == 7 && return @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇]

    y3 = y^3

    ψ₈  = zero(x)
    ψ₉  = 2*y
    ψ₁₀ = -3*y*(3 + 2*logx)
    ψ₁₁ = 36*x2*y - 8*y3
    ψ₁₂ = -40*y*(3*x2 + 6*y2 - 18*x2*logx + 4*y2*logx)

    ψ_polys = @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇, ψ₈, ψ₉, ψ₁₀, ψ₁₁, ψ₁₂]

    return ψ_polys
end

@fastmath function _solovev_polynomials_Dy(x, y, ::Val{N} = Val(12); logx=log(x)) where N

    x2 = x^2
    x3 = x^3
    x4 = x2^2
    x6 = x2^3
    y2 = y^2
    y3 = y^3
    y4 = y2^2

    # Derivatives from Mathematica
    ψ₁  = zero(y)
    ψ₂  = zero(y)
    ψ₃  = 2*y
    ψ₄  = -8*x2*y
    ψ₅  = 2*y*(-9*x2 + 4*y2 - 12*x2*logx)
    ψ₆  = -24*x4*y + 32*x2*y3
    ψ₇  = 2*y*(75*x4 - 280*x2*y2 + 24*y4 + 60*(3*x4 - 4*x2*y2)*logx)

    N == 7 && return @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇]

    ψ₈  = one(y)
    ψ₉  = x2
    ψ₁₀ = 3*(y2 - x2*logx)
    ψ₁₁ = 3*(x4 - 4*x2*y2)
    ψ₁₂ = 5*(-9*x4 + 8*y4 + 12*(x4 - 4*x2*y2)*logx)

    ψ_polys = @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇, ψ₈, ψ₉, ψ₁₀, ψ₁₁, ψ₁₂]

    return ψ_polys
end

@fastmath function _solovev_polynomials_Dyy(x, y, ::Val{N} = Val(12); logx=log(x)) where N

    x2 = x^2
    x3 = x^3
    x4 = x2^2
    x6 = x2^3
    y2 = y^2
    y3 = y^3
    y4 = y2^2

    # Derivatives from Mathematica
    ψ₁  = zero(y)
    ψ₂  = zero(y)
    ψ₃  = 2*one(y)
    ψ₄  = -8*x2
    ψ₅  = -6*(3*x2 - 4*y2 + 4*x2*logx)
    ψ₆  = -24*x2*(x2 - 4*y2)
    ψ₇  = 30*(5*x4 - 56*x2*y2 + 8*y4 + 12*(x4 - 4*x2*y2)*logx)

    N == 7 && return @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇]

    ψ₈  = zero(y)
    ψ₉  = zero(y)
    ψ₁₀ = 6*y
    ψ₁₁ = -24*x2*y
    ψ₁₂ = 160*y*(y2 - 3*x2*logx)

    ψ_polys = @SVector [ψ₁, ψ₂, ψ₃, ψ₄, ψ₅, ψ₆, ψ₇, ψ₈, ψ₉, ψ₁₀, ψ₁₁, ψ₁₂]

    return ψ_polys
end

@fastmath function _solovev_psi_p(x,y,A; logx=log(x))
    # Derivatives from Mathematica
    x2 = x^2
    x4 = x2*x2
    return x4/8 + A*(x2*logx/2 - x4/8)
end

@fastmath function _solovev_psi_p_Dx(x,y,A; logx=log(x))
    # Derivatives from Mathematica
    x2 = x^2
    return (x*(A + x2 - A*x2 + 2*A*logx))/2
end

@fastmath function _solovev_psi_p_Dxx(x,y,A; logx=log(x))
    # Derivatives from Mathematica
    x2 = x^2
    return (3*(A + x2 - A*x2))/2 + A*logx
end

function _solovev_psi_p_Dy(x,y,A)
    return zero(y)
end

function _solovev_psi_p_Dyy(x,y,A)
    return zero(y)
end

function _solovev_psi(x,y,A,coeffs::SVector{N};logx=log(x)) where N

    ψ_p = _solovev_psi_p(x,y,A; logx=logx)
    ψ_polys = _solovev_polynomials(x,y,Val(N),logx=logx)

    ψ_h = zero(x)
    @inbounds for i=1:N
        ψ_h += coeffs[i]*ψ_polys[i]
    end

    return ψ_p + ψ_h
end

function _solovev_psi_Dx(x,y,A,coeffs::SVector{N};logx=log(x)) where N

    ψ_p = _solovev_psi_p_Dx(x,y,A,logx=logx)
    ψ_polys = _solovev_polynomials_Dx(x,y,Val(N),logx=logx)

    ψ_h = zero(x)
    @inbounds for i=1:N
        ψ_h += coeffs[i]*ψ_polys[i]
    end

    return ψ_p + ψ_h
end

function _solovev_psi_Dxx(x,y,A,coeffs::SVector{N};logx=log(x)) where N

    ψ_p = _solovev_psi_p_Dxx(x,y,A,logx=logx)
    ψ_polys = _solovev_polynomials_Dxx(x,y,Val(N),logx=logx)

    ψ_h = zero(x)
    @inbounds for i=1:N
        ψ_h += coeffs[i]*ψ_polys[i]
    end

    return ψ_p + ψ_h
end

function _solovev_psi_Dy(x,y,A,coeffs::SVector{N};logx=log(x)) where N

    ψ_p = _solovev_psi_p_Dy(x,y,A)
    ψ_polys = _solovev_polynomials_Dy(x,y,Val(N),logx=logx)

    ψ_h = zero(y)
    @inbounds for i=1:N
        ψ_h += coeffs[i]*ψ_polys[i]
    end

    return ψ_p + ψ_h
end

function _solovev_psi_Dyy(x,y,A,coeffs::SVector{N};logx=log(x)) where N

    ψ_p = _solovev_psi_p_Dyy(x,y,A)
    ψ_polys = _solovev_polynomials_Dyy(x,y,Val(N),logx=logx)

    ψ_h = zero(x)
    @inbounds for i=1:N
        ψ_h += coeffs[i]*ψ_polys[i]
    end

    return ψ_p + ψ_h
end

"""
SolovevEquilibrium Structure

Defines the equilibrium that satisfy Δ⋆ψ(x,y) = α + (1 - α)x^2\\
where (r,z) = (R0 x,R0 y). F(dF/dψ) = -A, μ₀dp/dψ = -C, α = A/(A + C R0²), ψ₀ = R0²(A + C R0²)

Fields:\\
`cocos` - COCOS\\
`B0` - Toroidal magnetic field magnitude on axis [T]\\
`R0` - Major Radius [m]\\
`epsilon` - inverse aspect ratio a/R0\\
`delta` - triangularity\\
`kappa` - elongation/ellipticity\\
`alpha` - constant relating beta regime (α)\\
`qstar` - Kink safety factor\\
`psi0` - Poloidal flux normalization (ψ₀)\\
`beta_p` - Poloidal beta\\
`beta_t` - Toroidal beta\\
`c` - Coefficients for Solov'ev polynomials\\
`diverted` - If true then equilibrium has one or more x-points\\
`symmetric` - If true then equilibrium is up-down symmetric\\
`sigma` - Sign of dot(B,J)
"""
struct SolovevEquilibrium{T,N} <: AbstractEquilibrium
    cocos::COCOS
    B0::T
    R0::T
    epsilon::T
    delta::T
    kappa::T
    alpha::T
    qstar::T
    psi0::T
    beta_p::T
    beta_t::T
    c::SVector{N,T}
    diverted::Bool
    symmetric::Bool
    sigma_B0::Int
    sigma_Ip::Int
end

"""
    solovev(B0, R0, ϵ, δ, κ, q⋆; B0_dir=1, Ip_dir=1, diverted=false, xpoint=nothing, symmetric = true)

Creates a SolovevEquilibrium Structure

Arguments:\\
`B0` - Toroidal magnetic field magnitude on axis [T]\\
`R0` - Major Radius [m]\\
`ϵ` - inverse aspect ratio a/R0\\
`δ` - triangularity\\
`κ` - elongation/ellipticity\\
`α` - constant relating beta regime (α)\\
`q⋆` - Kink safety factor\\

Keyword Arguments:\\
`B0_dir` - Direction of the toroidal field (+1 = CCW, -1 = CW)\\
`Ip_dir` - Direction of the plasma current (+1 = CCW, -1 = CW)\\
`diverted` - If true then equilibrium has one or more x-points\\
`xpoint` - If diverted = true, then xpoint is set to (R0*(1-1.1*δ*ϵ), -R0*1.1*κ*ϵ) else nothing\\
`symmetric` - Is equilibrium up-down symmetric\\
"""
function solovev(B0, R0, ϵ, δ, κ, α, qstar;
                 B0_dir = 1, Ip_dir = 1,
                 diverted::Bool = false,
                 xpoint::Union{NTuple{2},Nothing} = (diverted ? (R0*(1-1.1*δ*ϵ),-R0*1.1*κ*ϵ) : nothing),
                 symmetric::Bool = (xpoint === nothing))

    if δ > sin(1)
        @warn "Equilibrium is not convex. δ > sin(1)"
    end

    diverted = xpoint !== nothing

    #Eq. 11
    δ₀ = asin(δ)
    N₁ = -(1 + δ₀)^2/(ϵ*κ^2) # [d²x/dy²]_(τ=0)
    N₂ =  (1 - δ₀)^2/(ϵ*κ^2) # [d²x/dy²]_(τ=pi)
    N₃ = -κ/(ϵ*cos(δ₀)^2)    # [d²y/dx²]_(τ=pi/2)

    Ψ_h = zeros(typeof(α),12,12)
    Ψ_p = zeros(typeof(α),12)
    # Boundary Conditions: (ψ_p - bc) + ψ_h'*c = 0
    if symmetric
        if xpoint === nothing
            # Eq. 10
            # Outer equatorial point
            x, y = 1 + ϵ, 0
            Ψ_p[1] = _solovev_psi_p(x, y, α)
            Ψ_h[:,1] .= _solovev_polynomials(x,y)
            # Inner equatorial point
            x, y = 1 - ϵ, 0
            Ψ_p[2] = _solovev_psi_p(x, y, α)
            Ψ_h[:,2] .= _solovev_polynomials(x, y)
            # High point
            x, y = 1 - δ*ϵ, κ*ϵ
            Ψ_p[3] = _solovev_psi_p(x, y, α)
            Ψ_h[:,3] .= _solovev_polynomials(x, y)
            # High point maximum
            x, y = 1 - δ*ϵ, κ*ϵ
            Ψ_p[4] = _solovev_psi_p_Dx(x, y, α)
            Ψ_h[:,4] .= _solovev_polynomials_Dx(x, y)
            # Outer equatorial point curvature
            x, y = 1 + ϵ, 0
            Ψ_p[5] = _solovev_psi_p_Dyy(x, y, α) + N₁*_solovev_psi_p_Dx(x, y, α)
            Ψ_h[:,5] .= _solovev_polynomials_Dyy(x, y) .+ N₁*_solovev_polynomials_Dx(x, y)
            # Inner equatorial point curvature
            x, y = 1 - ϵ, 0
            Ψ_p[6] = _solovev_psi_p_Dyy(x, y, α) + N₂*_solovev_psi_p_Dx(x, y, α)
            Ψ_h[:,6] .= _solovev_polynomials_Dyy(x, y) .+ N₂*_solovev_polynomials_Dx(x, y)
            # High point curvature
            x, y = 1 - δ*ϵ, κ*ϵ
            Ψ_p[7] = _solovev_psi_p_Dxx(x, y, α) + N₃*_solovev_psi_p_Dy(x, y, α)
            Ψ_h[:,7] .= _solovev_polynomials_Dxx(x, y) .+ N₃*_solovev_polynomials_Dy(x, y)

            c = SVector{7}(Ψ_h[1:7,1:7]'\(-Ψ_p[1:7]))
        else
            xsep, ysep = xpoint[1]/R0, xpoint[2]/R0 # normalize xpoint

            # Eq. 12
            # Outer equatorial point
            x, y = 1 + ϵ, 0
            Ψ_p[1] = _solovev_psi_p(x, y, α)
            Ψ_h[:,1] .= _solovev_polynomials(x,y)
            # Inner equatorial point
            x, y = 1 - ϵ, 0
            Ψ_p[2] = _solovev_psi_p(x, y, α)
            Ψ_h[:,2] .= _solovev_polynomials(x, y)
            # High point
            x, y = xsep, ysep
            Ψ_p[3] = _solovev_psi_p(x, y, α)
            Ψ_h[:,3] .= _solovev_polynomials(x, y)
            # B_norm = 0 at high point
            x, y = xsep, ysep
            Ψ_p[4] = _solovev_psi_p_Dx(x, y, α)
            Ψ_h[:,4] .= _solovev_polynomials_Dx(x, y)
            # B_tang = 0 at high point
            x, y = xsep, ysep
            Ψ_p[5] = _solovev_psi_p_Dy(x, y, α)
            Ψ_h[:,5] .= _solovev_polynomials_Dy(x, y)
            # Outer equatorial point curvature
            x, y = 1 + ϵ, 0
            Ψ_p[6] = _solovev_psi_p_Dyy(x, y, α) + N₁*_solovev_psi_p_Dx(x, y, α)
            Ψ_h[:,6] .= _solovev_polynomials_Dyy(x, y) .+ N₁*_solovev_polynomials_Dx(x, y)
            # Inner equatorial point curvature
            x, y = 1 - ϵ, 0
            Ψ_p[7] = _solovev_psi_p_Dyy(x, y, α) + N₂*_solovev_psi_p_Dx(x, y, α)
            Ψ_h[:,7] .= _solovev_polynomials_Dyy(x, y) .+ N₂*_solovev_polynomials_Dx(x, y)

            c = SVector{7}(Ψ_h[1:7,1:7]'\(-Ψ_p[1:7]))
        end
    else
        xsep, ysep = xpoint[1]/R0, xpoint[2]/R0 # normalize xpoint
        if ysep > 0
            throw(ArgumentError("X-point should be below the midplane"))
        end

        # Eq. 28
        # Outer equatorial point
        x, y = 1 + ϵ, 0
        Ψ_p[1] = _solovev_psi_p(x, y, α)
        Ψ_h[:,1] .= _solovev_polynomials(x,y)
        # Inner equatorial point
        x, y = 1 - ϵ, 0
        Ψ_p[2] = _solovev_psi_p(x, y, α)
        Ψ_h[:,2] .= _solovev_polynomials(x, y)
        # Upper high point
        x, y = 1 - δ*ϵ, κ*ϵ
        Ψ_p[3] = _solovev_psi_p(x, y, α)
        Ψ_h[:,3] .= _solovev_polynomials(x, y)
        # Lower X-point
        x, y = xsep, ysep
        Ψ_p[4] = _solovev_psi_p(x, y, α)
        Ψ_h[:,4] .= _solovev_polynomials(x, y)
        # Outer equatorial point up-down symmetry
        x, y = 1 + ϵ, 0
        Ψ_p[5] = _solovev_psi_p_Dy(x,y, α)
        Ψ_h[:,5] .= _solovev_polynomials_Dy(x,y)
        # Inner equatorial point up-down symmetry
        x, y = 1 - ϵ, 0
        Ψ_p[6] = _solovev_psi_p_Dy(x, y, α)
        Ψ_h[:,6] .= _solovev_polynomials_Dy(x, y)
        # Upper high point maximum
        x, y = 1 - δ*ϵ, κ*ϵ
        Ψ_p[7] = _solovev_psi_p_Dx(x, y, α)
        Ψ_h[:,7] .= _solovev_polynomials_Dx(x, y)
        # B_y = 0 at X-point
        x, y = xsep, ysep
        Ψ_p[8] = _solovev_psi_p_Dx(x, y, α)
        Ψ_h[:,8] .= _solovev_polynomials_Dx(x, y)
        # B_x = 0 at X-point
        x, y = xsep, ysep
        Ψ_p[9] = _solovev_psi_p_Dy(x, y, α)
        Ψ_h[:,9] .= _solovev_polynomials_Dy(x, y)
        # Outer equatorial point curvature
        x, y = 1 + ϵ, 0
        Ψ_p[10] = _solovev_psi_p_Dyy(x, y, α) + N₁*_solovev_psi_p_Dx(x, y, α)
        Ψ_h[:,10] .= _solovev_polynomials_Dyy(x, y) .+ N₁*_solovev_polynomials_Dx(x, y)
        # Inner equatorial point curvature
        x, y = 1 - ϵ, 0
        Ψ_p[11] = _solovev_psi_p_Dyy(x, y, α) + N₂*_solovev_psi_p_Dx(x, y, α)
        Ψ_h[:,11] .= _solovev_polynomials_Dyy(x, y) .+ N₂*_solovev_polynomials_Dx(x, y)
        # High point curvature
        x, y = 1 - δ*ϵ, κ*ϵ
        Ψ_p[12] = _solovev_psi_p_Dxx(x, y, α) + N₃*_solovev_psi_p_Dy(x, y, α)
        Ψ_h[:,12] .= _solovev_polynomials_Dxx(x, y) .+ N₃*_solovev_polynomials_Dy(x, y)

        c = SVector{12}(Ψ_h'\(-Ψ_p))
    end

    τ = range(0,2pi,length=100)
    x = (1 .+ ϵ .* cos.(τ .+ δ₀*sin.(τ)))
    y = (ϵ*κ*sin.(τ))
    bdry = PlasmaBoundary(collect(zip(x,y)))

    #Eq. 6.158 in Ideal MHD
    V = volume(bdry,dx=0.005,dy=0.005)
    K₁ = area_average(bdry, (x,y) -> (α + (1-α)*x^2)/x,dx=0.005,dy=0.005)*area(bdry)
    K₂ = volume_average(bdry,(x,y) -> -_solovev_psi(x,y,α,c),dx=0.005,dy=0.005)*V/(2pi)
    K₃ = V/(2pi)

    #Eq. 19 from paper (error in β_t in Ideal MHD (extra ^2))
    Cp = circumference(bdry)
    a = ϵ*R0
    ψ0 = -Cp*a*R0*B0/(qstar*K₁*Ip_dir)
    β_p = 2*(1-α)*Cp^2 * K₂/(K₁^2*K₃)
    β_t = ϵ^2*β_p/(qstar^2)

    cc = 3
    return SolovevEquilibrium(cocos(cc),promote(B0,R0,ϵ,δ,κ,α,qstar,ψ0,β_p,β_t)...,c, diverted, symmetric, B0_dir, Ip_dir)
end

function (S::SolovevEquilibrium)(r,z)
    x = r/S.R0
    y = z/S.R0
    logx = log(x)
    return S.psi0*_solovev_psi(x,y,S.alpha,S.c,logx=logx)
end

function Base.show(io::IO, S::SolovevEquilibrium)
    print(io, "SolovevEquilibrium\n")
    print(io, "  B0 = $(S.sigma_B0*S.B0) [T]\n")
    print(io, "  R0 = $(S.R0) [m]\n")
    print(io, "  ϵ  = $(S.epsilon)\n")
    print(io, "  δ  = $(S.delta)\n")
    print(io, "  κ  = $(S.kappa)\n")
    print(io, "  α  = $(S.alpha)\n")
    print(io, "  q⋆ = $(S.qstar)\n")
    print(io, "  βp = $(S.beta_p)\n")
    print(io, "  βt = $(S.beta_t)\n")
    print(io, "  σ  = $(S.sigma_B0*S.sigma_Ip)\n")
    print(io, "  diverted  = $(S.diverted)\n")
    print(io, "  symmetric = $(S.symmetric)")
end

Base.broadcastable(S::SolovevEquilibrium) = (S,)

function cocos(S::SolovevEquilibrium)
    return S.cocos
end

function B0Ip_sign(S::SolovevEquilibrium)
    return S.sigma_B0*S.sigma_Ip
end

function limits(S::SolovevEquilibrium)
    xlims = (S.R0*(1-1.2*S.epsilon), S.R0*(1 + 1.2*S.epsilon))
    ylims = (S.R0*(-1.2*S.epsilon*S.kappa), S.R0*(1.2*S.epsilon*S.kappa))
    return xlims, ylims
end

function boundary(S::SolovevEquilibrium; n=100, kwargs...)
    if S.diverted
        return boundary(S,0.0; kwargs...)
    else
        τ = range(0,2pi,length=n)
        δ₀ = asin(S.delta)
        x = S.R0*(1 .+ S.epsilon .* cos.(τ .+ δ₀*sin.(τ)))
        y = S.R0*(S.epsilon*S.kappa*sin.(τ))
        return PlasmaBoundary(collect(zip(x,y)))
    end
end

function psi_gradient(S::SolovevEquilibrium,r,z)
    x, y = r/S.R0, z/S.R0
    logx=log(x)
    return (S.psi0/S.R0)*SVector{2}(_solovev_psi_Dx(x,y,S.alpha,S.c; logx=logx), _solovev_psi_Dy(x,y,S.alpha,S.c;logx=logx))
end

_solovev_magnetic_axis = Dict{SolovevEquilibrium,NTuple{2}}()
function clear_cache(S)
    delete!(_solovev_magnetic_axis,S)
end

function magnetic_axis(S::SolovevEquilibrium; x0 = (S.R0, zero(S.R0)))
    if S in keys(_solovev_magnetic_axis)
        return _solovev_magnetic_axis[S]
    else
        sigma_psi = sign(S(x0[1],x0[2]))
        res = Optim.optimize(x->-sigma_psi*S(x[1],x[2]), x -> -sigma_psi*psi_gradient(S,x[1],x[2]), collect(x0), inplace=false)
        axis = (res.minimizer[1],res.minimizer[2])
        _solovev_magnetic_axis[S] = axis
    end
    return axis
end

function psi_limits(S::SolovevEquilibrium)
    psimag = S(magnetic_axis(S)...)
    psibry = zero(psimag)
    return (psimag, psibry)
end

function pressure_gradient(S::SolovevEquilibrium)
    C = S.psi0*(1-S.alpha)/(S.R0^4)
    return -C/mu0
end

function pressure_gradient(S::SolovevEquilibrium, psi)
    return pressure_gradient(S) + psi*0
end

function pressure(S::SolovevEquilibrium, psi; p0=zero(psi))
    C = S.psi0*(1-S.alpha)/(S.R0^4)
    dpdpsi = -C/mu0
    p = dpdpsi*psi + p0
    return p
end

function poloidal_current(S::SolovevEquilibrium,psi)
    A = S.psi0*S.alpha/(S.R0^2)
    dF2dpsi = -2*A
    F2 = abs(dF2dpsi*psi + S.R0^2*S.B0^2)
    return S.sigma_B0*sqrt(F2)
end

function poloidal_current_gradient(S::SolovevEquilibrium,psi)
    A = S.psi0*S.alpha/(S.R0^2)
    return -A/poloidal_current(S,psi)
end

electric_potential(S::SolovevEquilibrium,psi) = zero(psi)
electric_potential_gradient(S::SolovevEquilibrium,psi) = zero(psi)

# === Special Cases ===

function curlB(S::SolovevEquilibrium, r, z)
    # curlB = mu0*J
    psi = S(r,z)
    gval = poloidal_current(S,psi)
    grad_psi = psi_gradient(S,r,z)

    gp = poloidal_current_gradient(S,psi)
    pp = pressure_gradient(S)

    cc = cocos(S)

    invr = inv(r)
    cocos_factor = cc.sigma_RpZ*invr
    cBr = -cocos_factor*gp*grad_psi[2]
    cBz =  cocos_factor*gp*grad_psi[1]

    cocos_factor = -cc.sigma_Bp*((2pi)^cc.exp_Bp)
    cBt = cocos_factor*(r*pp*mu0 + gval*gp*invr)

    return SVector{3}(cylindrical_cocos(cc,cBr,cBt,cBz))
end

function curlB(S::SolovevEquilibrium, x, y, z)
    r = hypot(x,y)
    phi = atan(y,x)
    cB = curlB(S,r,z)
    sp, cp = sincos(phi)
    cB_xyz = SVector{3}(cB[1]*cp - cB[2]*sp, cB[1]*sp + cB[2]*cp, cB[3])
    return cB_xyz
end


text:
module myTesting
using SpinLattice
using PyCall, RecursiveArrayTools
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport mkl
mkl.set_num_threads(1)
@pyimport madness as oop
@pyimport scipy.sparse as pysparse
@pyimport scipy.sparse.linalg as slinalg

L1d_1 = ((15,), (-0.41,-0.41,0.82), Interaction(nearest_neighbours), (0.,0.,0.))
L1d_2 = ((15,), (1.0,1.0,1.0), Interaction(nearest_neighbours), (0.,0.,0.))
L1d_3 = ((15,), (0.707,0.707,0.0), Interaction(nearest_neighbours), (0.,0.,0.))
L1d_4 = ((15,), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.,0.,0.))
L1d_5 = ((15,), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.207,0.,0.))
L1d_6 = ((15,), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.207,0.,0.100))
L1d_7 = ((15,), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.207,0.,0.207))

L2d_1 = ((9,9), (-0.41,-0.41,0.82), Interaction(nearest_neighbours), (0.,0.,0.))
L2d_2 = ((9,9), (1.0,1.0,1.0), Interaction(nearest_neighbours), (0.,0.,0.))
L2d_3 = ((9,9), (0.707,0.707,0.0), Interaction(nearest_neighbours), (0.,0.,0.))
L2d_4 = ((9,9), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.,0.,0.))
L2d_5 = ((9,9), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.207,0.,0.))
L2d_6 = ((9,9), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.207,0.,0.100))
L2d_7 = ((9,9), (0.518, 0.830, 0.207), Interaction(nearest_neighbours), (0.207,0.,0.207))

L3d_1 = ((9,9,9), (-0.41,-0.41,0.82), Interaction(nearest_neighbours), (0.,0.,0.))
Lfid0x0x1 = ((9,9,9), (-0.5,-0.5,1.0), Interaction(MDI(0,0,1)), (0.,0.,0.), [(1,1,i) for i in 1:9], [(1,1,i) for i in 1:9])
Lfid0x1x1 = ((9,9,9), (-0.5,-0.5,1.0), Interaction(MDI(0,1,1)), (0.,0.,0.), [(i,1,1) for i in 1:9], [(i,1,1) for i in 1:9])
Lfid1x1x1 = ((9,9,9), (-0.5,-0.5,1.0), Interaction(MDI(1,1,1)), (0.,0.,0.), [(i,i,i) for i in 1:9], [(i,i,i) for i in 1:9])

L1d = [L1d_1, L1d_2, L1d_3, L1d_4, L1d_5, L1d_6, L1d_7]
L2d = [L2d_1, L2d_2, L2d_3, L2d_4, L2d_5, L2d_6, L2d_7]
L3d = [L3d_1, Lfid0x0x1, Lfid0x1x1, Lfid1x1x1]

d1block = (5,)
d2block = (3,3)
d3block = (:cross, sort([(2,2,2),(2,2,1),(2,1,2),(1,2,2),(2,2,3),(2,3,2),(3,2,2)]), [(2,2,2)])
d3block_py = ([Tuple(reverse(collect(ind) .- 1)) for ind in d3block[2]], [Tuple(reverse(collect(ind) .- 1)) for ind in d3block[3]])
get_block(::Type{Val{1}}) = d1block
get_block(::Type{Val{2}}) = d2block
get_block(::Type{Val{3}}) = d3block
get_block_py(x::Type{Val{D}}) where D = get_block(x)
get_block_py(::Type{Val{3}}) = d3block_py

jlmat2pymat(S::Union{SparseMatrixCSC, SpinLattice.SharedSparseMatrixCSC}) = pysparse.csc_matrix((S.nzval, S.rowval .- 1, S.colptr .- 1), shape=size(S))
function py_initialize(M::Exact, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}; kwargs...) where {D}
    sp = oop.UniformLattice(reverse(dims), hs, Js; kwargs...)
    return sp
end
function py_initialize(M::Exact, dims::NTuple{2,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}; kwargs...)
    sp = oop.UniformLattice(reverse((4,4)), hs, Js; kwargs...)
    return sp
end
function py_initialize(M::Clustered, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}; kwargs...) where {D}
    dims_r = reverse(dims)
    block_dims = reverse(get_block_py(Val{D}))
    blocks = Tuple(div.(collect(dims_r), collect(block_dims)))
    sp = oop.ClusteredLattice(block_dims, blocks, hs, Js; kwargs...)
    return sp
end
function py_initialize(M::PureClassical, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}; kwargs...) where {D}
    sp = oop.ClassicalLattice(reverse(dims), hs, Js; kwargs...)
    return sp
end
function py_initialize(M::Hybrid, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}; kwargs...) where {D}
    dims_r = reverse(dims)
    block_dims = reverse(get_block_py(Val{D}))
    sp = oop.HybridLattice(block_dims, dims_r, hs, Js; kwargs...)
    return sp
end
function py_initialize(M::Hybrid, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}; kwargs...)
    q_spins, middle_spins = d3block_py
    q_spins_r = [reverse(x) for x in q_spins]
    sp = oop.HybridLatticeIR(q_spins_r, reverse(dims), middle_spins, hs, Js; kwargs...)
    return sp
end
function py_initialize(M::Hybrid, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{MDI{3}}, hs::NTuple{3,Float64}; kwargs...)
    q_spins, middle_spins = d3block_py
    q_spins_r = [reverse(x) for x in q_spins]
    sp = oop.FIDLattice(3, reverse(dims), q_spins_r, middle_spins, direction = reverse(I.func.dir); kwargs...)
    return sp
end
function py_initialize(M::Hybrid, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{MDI{3}}, hs::NTuple{3,Float64}, q_spins::Vector{NTuple{3,Int64}}, middle_spins::Vector{NTuple{3,Int64}}; kwargs...)
    q_spins_r = [Tuple(reverse(collect(ind) .- 1)) for ind in q_spins]
    middle_spins_r = [Tuple(reverse(collect(ind) .- 1)) for ind in middle_spins]
    sp = oop.FIDLattice(3, reverse(dims), q_spins_r, middle_spins_r, direction = reverse(I.func.dir); kwargs...)
    return sp
end
function py_initialize(M::PureClassical, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{MDI{3}}, hs::NTuple{3,Float64}; kwargs...)
    sp = oop.FIDLatticeClassical(3, reverse(dims), direction = reverse(I.func.dir); kwargs...)
    return sp
end
@inline function py_initialize(M::PureClassical, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, q_spins::Vector{NTuple{3,Int64}}, middle_spins::Vector{NTuple{3,Int64}}; kwargs...) where {D,F}
    py_initialize(M::PureClassical, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}; kwargs...)
end

@inline default_links(L::Lattice, M::Model) where Model<:Union{Exact, PureClassical} = (:all,:all)
@inline default_links(L::Lattice{D}, M::Model) where {D, Model<:Union{Clustered,Hybrid}} = (:all,:central)
@inline default_links(L::Lattice{3}, M::Hybrid) = (:all, SpinArray(L, [(2,2,2)], :central))

function julia_initialize(M::Model, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, B::Bool=true) where {D, Model<:Union{Exact,PureClassical}} where F
    L = Lattice(dims, Js, I, hs, Int64, B)
    A = build_Approximation(M, L)
    RHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, default_links(L,M))
    return A, RHS, OSET, rules
end
@inline function julia_initialize(M::PureClassical, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, q_spins::Vector{NTuple{3,Int64}}, middle_spins::Vector{NTuple{3,Int64}}, B::Bool=true) where {D, F}
    julia_initialize(M, dims, Js, I, hs, Int64, B)
end
function julia_initialize(M::Exact, dims::NTuple{2,Int64}, Js::NTuple{3,Float64}, I::Interaction{typeof(nearest_neighbours)}, hs::NTuple{3,Float64}, B::Bool=true)
    L = Lattice((4,4), Js, I, hs, Int64, B)
    A = build_Approximation(M, L)
    RHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, default_links(L,M))
    return A, RHS, OSET, rules
end

function julia_initialize(M::Model, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, B::Bool=true) where {D, Model<:Union{Clustered,Hybrid}, F}
    L = Lattice(dims, Js, I, hs, Int64, B)
    A = build_Approximation(M, L, get_block(Val{D}))
    RHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, default_links(L,M))
    return A, RHS, OSET, rules
end

function julia_initialize(M::Hybrid, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, B::Bool=true) where F
    L = Lattice(dims, Js, I, hs, Int64, B)
    name, q_spins, middle_spins = d3block
    sarray = SpinArray(L, q_spins, name)
    A = build_Approximation(M, L, sarray)
    RHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, default_links(L,M))
    return A, RHS, OSET, rules
end
function julia_initialize(M::Hybrid, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, q_spins::Vector{NTuple{3,Int64}}, middle_spins::Vector{NTuple{3,Int64}}, B::Bool=true) where F
    L = Lattice(dims, Js, I, hs, Int64, B)
    name = :column
    sarray = SpinArray(L, q_spins, name)
    A = build_Approximation(M, L, sarray)
    RHS = build_RHS_function(A)
    rules, OSET = set_CorrFuncs(A, (:all, SpinArray(L, middle_spins, :central)))
    return A, RHS, OSET, rules
end

function test_julia_init()
    Models = [Exact(), PureClassical(), Clustered(), Hybrid()]
    for model in Models
        for lat in L1d
            try
                julia_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
    for model in Models
        for lat in L2d
            try
                julia_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
    for model in Models[[2,4]]
        for lat in L3d
            try
                julia_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
end
function test_py_init()
    Models = [Exact(), PureClassical(), Clustered(), Hybrid()]
    for model in Models
        for lat in L1d
            try
                py_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
    for model in Models
        for lat in L2d
            try
                py_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
    for model in Models[[2,4]]
        for lat in L3d
            try
                py_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
end
function test_py_init_3d()
    Models = [Exact(), PureClassical(), Clustered(), Hybrid()]
    for model in Models[[2,4]]
        for lat in L3d
            try
                py_initialize(model, lat...)
                print("all green\n\n")
            catch
                print("$model   $lat\n\n")
            end
        end
    end
end

const toler = 1e-14

function convert_state(A::ExactApprox, sp)
    v = Vector{Vector{Complex{Float64}}}()
    push!(v, copy(sp[:psi]))
    for i = 1:size(sp[:psi_aid],1)
        push!(v, sp[:psi_aid][i,:])
    end
    return VectorOfArray(v)
end

function test_operators(M::Exact, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, B::Bool=true) where {D, F}
    A,RHS,OSET, rules = julia_initialize(M, dims, Js, I, hs, B)
    sp = py_initialize(M, dims, Js, I, hs; cb=B)
    H_j = jlmat2pymat(RHS.H)
    Ms_j = Dict()
    for Obs in OSET.Observables
        axis = axis_dict[Obs.axis]-1
        Ms_j[axis] = jlmat2pymat(Obs.QO.O)
    end
    sp[:initState]()
    u0 = copy(sp[:psi])
    du = copy(u0)
    tj = -time()
    RHS(0.0,u0,du)
    tj += time()
    tp = -time()
    du2 = sp[:rhs](sp[:psi])
    tp += time()
    print("$((dims, Js, I, hs))\n")
    @printf("Check Hamiltonian:   %e\n", slinalg.norm(sp[:H]-H_j)/slinalg.norm(sp[:H]))
    @printf("Check tot magnetiz: %e   %e   %e\n", [slinalg.norm(sp[:Ms][i]-Ms_j[i])/slinalg.norm(sp[:Ms][i]) for i in 0:2]...)
    @printf("Check RHS: %e\n", norm(du.*(2.0^-7)-du2)/norm(du2))
    @printf("Check time: %10e", tj/tp)
    print("\n\n")

end

function convert_state(A::ClusteredApprox, sp)
    cluster_num = size(sp[:psi],1)
    u0 = VectorOfArray([copy(sp[:psi][i,:]) for i in 1:cluster_num])
    return u0
end

function test_operators(M::Clustered, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}) where {D, F}
    d = Dict([(:all,"Global"), (:central,"Local")])
    A,RHS,OSET, rules = julia_initialize(M, dims, Js, I, hs)
    sp = py_initialize(M, dims, Js, I, hs)
    H_j = jlmat2pymat(RHS.H)
    ops_j = Dict()
    for Obs in OSET.Observables
        index = (string(Obs.axis), d[Obs.name])
        ops_j[index] = jlmat2pymat(Obs.QO.O)
    end
    num_of_espins = RHS.num_of_espins
    taus_j = Matrix(num_of_espins, 3)
    for espin in 1:num_of_espins for sigma in 1:3
        taus_j[espin, sigma] = jlmat2pymat(RHS.edgeSpinOperators[sigma, espin])
    end end
    sp[:initState]()
    cluster_num = size(sp[:psi],1)
    u0 = VectorOfArray([copy(sp[:psi][i,:]) for i in 1:cluster_num])
    du = VectorOfArray([copy(sp[:psi][i,:]) for i in 1:cluster_num])
    tj = -time()
    RHS(0.0,u0,du)
    tj += time()
    tp = -time()
    du2 = sp[:rhs](sp[:psi])
    tp += time()

    print("$((dims, Js, I, hs))\n")
    @printf("Check Hamiltonian:   %e\n", slinalg.norm(sp[:H]-H_j)/slinalg.norm(sp[:H]))
    print("Check spin operators:\n")
    for espin in 1:num_of_espins
        for sigma in 1:3
            @printf("%10e ", slinalg.norm(sp[:taus][espin, sigma]-taus_j[espin, sigma])/slinalg.norm(sp[:taus][espin, sigma]))
        end
        print("\n")
    end
    ###################
    links = sp[:links]
    noe = num_of_espins
    cl_num = cluster_num
    len = noe*cl_num
    IMp = zeros(len, len)
    for link in links
        p1 = sub2ind((noe,cl_num), link[1][2]+1, link[1][1]+1)
        p2 = sub2ind((noe,cl_num), link[2][2]+1, link[2][1]+1)
        IMp[p1,p2] += 1.0
    end
    IMp .= .-IMp.*0.25.*sqrt(sp[:N])
    @printf("Check intrablock matrix: %10e\n", norm(IMp-RHS.IM)/norm(IMp))
    ####################
    print("Check observables:\n")
    print("Local: ")
    for axis in ["x","y","x"]
        index = (axis,"Local")
        @printf("%10e ", slinalg.norm(sp[:ops][index] - ops_j[index])/slinalg.norm(sp[:ops][index]))
    end
    print("\n")
    print("Global: ")
    for axis in ["x","y","x"]
        index = (axis,"Global")
        @printf("%10e ", slinalg.norm(sp[:ops][index] - ops_j[index])/slinalg.norm(sp[:ops][index]))
    end
    print("\n")
    ###################
    print("Check rhs: ")
    for i in 1:cluster_num
        @printf("%10e ", norm(du[:,i].*2.0^-7-du2[i,:])/norm(du2[i,:]))
    end
    print("\n")
    @printf("Check time: %10e", tj/tp)
    print("\n\n")
end

function convert_state(A::HybridApprox, sp)
    u0 = ArrayPartition(copy(sp[:psi]), VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3]))
    return u0
end

function test_operators(M::Hybrid, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}) where {D, F}
    d = Dict([(:all,"Global"), (:central,"Local")])
    A,RHS,OSET, rules = julia_initialize(M, dims, Js, I, hs)
    sp = py_initialize(M, dims, Js, I, hs)
    H_j = jlmat2pymat(RHS.H)
    ops_j = Dict()
    for Obs in OSET.Observables
        index = (axis_dict[Obs.axis]-1, d[Obs.name])
        ops_j[index] = (jlmat2pymat(Obs.QO.O),Obs)
    end
    num_of_espins = size(RHS.edgeSpinOperators,1)
    taus_j = Matrix(num_of_espins, 3)
    for espin in 1:num_of_espins for sigma in 1:3
        taus_j[espin, sigma] = jlmat2pymat(RHS.edgeSpinOperators[espin, sigma])
    end end
    sp[:initState]()
    sp[:updateCorrelator]()
    u0 = ArrayPartition(copy(sp[:psi]), VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3]))
    du = ArrayPartition(copy(sp[:psi]), VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3]))
    tj = -time()
    RHS(0.0,u0,du)
    tj += time()
    tp = -time()
    du2 = sp[:rhs](sp[:psi], sp[:clState])
    tp += time()

    print("$((dims, Js, I, hs))\n")
    @printf("Check Hamiltonian:   %e\n", slinalg.norm(sp[:H]-H_j)/slinalg.norm(sp[:H]))
    print("Check spin operators:\n")
    for espin in 1:num_of_espins
        for sigma in 1:3
            @printf("%10e ", slinalg.norm(sp[:taus][espin, sigma]-taus_j[espin, sigma])/slinalg.norm(sp[:taus][espin, sigma]))
        end
        print("\n")
    end
    ###################
    print("Check observable ops:\n")
    print("Local: ")
    for axis in [0,1,2]
        index = (axis,"Local")
        @printf("%10e ", slinalg.norm(sp[:ops][index] - ops_j[index][1])/slinalg.norm(sp[:ops][index]))
    end
    print("\n")
    print("Global: ")
    for axis in [0,1,2]
        index = (axis,"Global")
        @printf("%10e ", slinalg.norm(sp[:ops][index] - ops_j[index][1])/slinalg.norm(sp[:ops][index]))
    end
    print("\n")
    ###################
    print("Check observable vals:\n")
    print("Local: ")
    for axis in [0,1,2]
        index = (axis,"Local")
        @printf("%10e ", norm(sp[:means][index][end] - ops_j[index][2](u0))/norm(sp[:means][index][end]))
    end
    print("\n")
    print("Global: ")
    for axis in [0,1,2]
        index = (axis,"Global")
        @printf("%10e ", norm(sp[:means][index][end] - ops_j[index][2](u0))/norm(sp[:means][index][end]))
    end
    print("\n")
    ###################
    IM_j = jlmat2pymat(RHS.IM)
    q2cl_j = jlmat2pymat(RHS.q2cl)
    cl2q_j = jlmat2pymat(RHS.cl2q)
    q2cl_p = sp[:Q2CL]*0.5*sqrt(get_Dh(A))
    cl2q_p = sp[:CL2Q]*-0.5
    @printf("Check IM: %10e\n", slinalg.norm(sp[:FM]-IM_j)/slinalg.norm(sp[:FM]))
    @printf("Check Q2CL: %10e\n", slinalg.norm(q2cl_p-q2cl_j)/slinalg.norm(q2cl_p))
    @printf("Check CL2Q: %10e\n", slinalg.norm(cl2q_p-cl2q_j)/slinalg.norm(cl2q_p))
    ###################
    print("Check rhs:\n")
    @printf("%10e\n", norm(du2[1] - du.x[1].*2.0^-7)/norm(du2[1]))
    for sigma in 1:3
        @printf("%10e ", norm(du2[2][sigma,:]-du.x[2][:,sigma].*2.0^-7)/norm(du2[2][sigma,:]))
    end
    print("\n")
    @printf("Check time: %10e", tj/tp)
    print("\n\n")
end

function test_operators(M::Hybrid, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, q_spins::Vector{NTuple{3,Int64}}, middle_spins::Vector{NTuple{3,Int64}}) where {F}
    d = Dict([(:all,"Global"), (:central,"Local")])
    A,RHS,OSET, rules = julia_initialize(M, dims, Js, I, hs, q_spins, middle_spins)
    sp = py_initialize(M, dims, Js, I, hs, q_spins, middle_spins)
    H_j = jlmat2pymat(RHS.H)
    ops_j = Dict()
    for Obs in OSET.Observables
        if Obs.axis == :x
            index = d[Obs.name]
            ops_j[index] = (jlmat2pymat(Obs.QO.O),Obs)
        end
    end
    num_of_espins = size(RHS.edgeSpinOperators,1)
    taus_j = Matrix(num_of_espins, 3)
    for espin in 1:num_of_espins for sigma in 1:3
        taus_j[espin, sigma] = jlmat2pymat(RHS.edgeSpinOperators[espin, sigma])
    end end
    sp[:initState]()
    sp[:updateCorrelator]()
    u0 = ArrayPartition(copy(sp[:psi]), VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3]))
    du = ArrayPartition(copy(sp[:psi]), VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3]))
    tj = -time()
    RHS(0.0,u0,du)
    tj += time()
    tp = -time()
    du2 = sp[:rhs](sp[:psi], sp[:clState])
    tp += time()

    print("$((dims, Js, I, hs))\n")
    @printf("Check Hamiltonian:   %e\n", slinalg.norm(sp[:H]-H_j)/slinalg.norm(sp[:H]))
    print("Check spin operators:\n")
    for espin in 1:num_of_espins
        for sigma in 1:3
            @printf("%10e ", slinalg.norm(sp[:taus][espin, sigma]-taus_j[espin, sigma])/slinalg.norm(sp[:taus][espin, sigma]))
        end
        print("\n")
    end
    ###################
    print("Check observable ops:\n")
    index = "Local"
    @printf("Local: %10e\n", slinalg.norm(sp[:ops][index] - ops_j[index][1])/slinalg.norm(sp[:ops][index]))
    index = "Global"
    @printf("Global: %10e\n", slinalg.norm(sp[:ops][index] - ops_j[index][1])/slinalg.norm(sp[:ops][index]))
    ###################
    print("Check observable vals:\n")
    index = "Local"
    @printf("Local: %10e\n", norm(sp[:means][index][end] - ops_j[index][2](u0))/norm(sp[:means][index][end]))
    index = "Global"
    @printf("Global: %10e\n", norm(sp[:means][index][end] - ops_j[index][2](u0))/norm(sp[:means][index][end]))
    ###################
    IM_j = RHS.IM
    q2cl_j = RHS.q2cl
    cl2q_j = RHS.cl2q
    q2cl_p = sp[:Q2CL]
    cl2q_p = -sp[:CL2Q]
    @printf("Check IM: %10e\n", norm(sp[:FM]-IM_j)/norm(sp[:FM]))
    @printf("Check Q2CL: %10e\n", norm(q2cl_p-q2cl_j)/norm(q2cl_p))
    @printf("Check CL2Q: %10e\n", norm(cl2q_p-cl2q_j)/norm(cl2q_p))
    ###################
    print("Check rhs:\n")
    @printf("%10e\n", norm(du2[1] - du.x[1].*2.0^-7)/norm(du2[1]))
    for sigma in 1:3
        @printf("%10e ", norm(du2[2][sigma,:]-du.x[2][:,sigma].*2.0^-7)/norm(du2[2][sigma,:]))
    end
    print("\n")
    @printf("Check time: %10e", tj/tp)
    print("\n\n")
end

function convert_state(A::PureClassicalApprox, sp)
    u0 = VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3])
    return u0
end

function test_operators(M::PureClassical, dims::NTuple{D,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, B::Bool=true) where {D, F}
    A,RHS,OSET, rules = julia_initialize(M, dims, Js, I, hs, B)
    sp = py_initialize(M, dims, Js, I, hs; cb=B)
    ops_j = Dict()
    for Obs in OSET.Observables
        index = axis_dict[Obs.axis]-1
        ops_j[index] = Obs
    end
    sp[:initState]()
    sp[:updateCorrelator]()
    u0 = VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3])
    du = VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3])
    tj = -time()
    RHS(0.0,u0,du)
    tj += time()
    tp = -time()
    du2 = sp[:rhs](sp[:clState])
    tp += time()

    print("$((dims, Js, I, hs))\n")
    ###################
    print("Check observable vals: ")
    for axis in [0,1,2]
        @printf("%10e ", norm(sp[:means][axis][end] - ops_j[axis](u0))/norm(sp[:means][axis][end]))
    end
    print("\n")
    ###################
    IM_j = jlmat2pymat(RHS.IM)
    @printf("Check IM: %10e\n", slinalg.norm(sp[:FM]-IM_j)/slinalg.norm(sp[:FM]))
    ###################
    print("Check rhs: ")
    for sigma in 1:3
        @printf("%10e ", norm(du2[sigma,:]-du[:,sigma].*2.0^-7)/norm(du2[sigma,:]))
    end
    print("\n")
    @printf("Check time: %10e", tj/tp)
    print("\n\n")
end

function test_operators(M::PureClassical, dims::NTuple{3,Int64}, Js::NTuple{3,Float64}, I::Interaction{F}, hs::NTuple{3,Float64}, q_spins, middle_spins) where {F}
    A,RHS,OSET, rules = julia_initialize(M, dims, Js, I, hs)
    sp = py_initialize(M, dims, Js, I, hs)
    ops_j = OSET.Observables[1]
    for Obs in OSET.Observables
        if Obs.axis == :x
            ops_j = Obs
        end
    end
    sp[:initState]()
    sp[:updateCorrelator]()
    u0 = VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3])
    du = VectorOfArray([copy(sp[:clState][i,:]) for i in 1:3])
    tj = -time()
    RHS(0.0,u0,du)
    tj += time()
    tp = -time()
    du2 = sp[:rhs](sp[:clState])
    tp += time()

    print("$((dims, Js, I, hs))\n")
    ###################
    @printf("Check observable val: %10e\n", norm(sp[:means][end] - ops_j(u0))/norm(sp[:means][end]))
    ###################
    IM_j = RHS.IM
    @printf("Check IM: %10e\n", norm(sp[:FM]-IM_j)/norm(sp[:FM]))
    ###################
    print("Check rhs: ")
    for sigma in 1:3
        @printf("%10e ", norm(du2[sigma,:]-du[:,sigma].*2.0^-7)/norm(du2[sigma,:]))
    end
    print("\n")
    @printf("Check time: %10e", tj/tp)
    print("\n\n")
end

@inline store_state(A::ExactApprox, sp) = (deepcopy(sp[:psi]),deepcopy(sp[:psi_aid]))
@inline store_state(A::ClusteredApprox, sp) = (deepcopy(sp[:psi]),)
@inline store_state(A::HybridApprox, sp) = (deepcopy(sp[:psi]), deepcopy(sp[:clState]))
@inline store_state(A::PureClassicalApprox, sp) = (deepcopy(sp[:clState]),)

function test_propagation(M::AbstractModel, args...; Tmax::Float64 = 10.0, tstep::Float64 = 2.0^-7, delimiter = 10)
    A, ARHS, OSET, rules = julia_initialize(M, args...)
    sp = py_initialize(M, args...; :Tmax => Tmax, :tstep => tstep)
    sp[:initState]()
    u0 = convert_state(A, sp)
    state = store_state(A, sp)
    sp[:propagate](state...)
    ####
    l = length(OSET.Observables)
    tspan = (0.0,Tmax*delimiter)
    t_means = collect(0.0:tstep:(Tmax*delimiter))
    t_cors = collect(0.0:tstep:Tmax)
    saved_values = SpinLattice.SavedValues(t_means, [zeros(l) for i in eachindex(t_means)])
    prob = SpinLattice.ODEProblem(ARHS, u0, tspan)
    cb = SpinLattice.SavingCallback(OSET, saved_values, u0; saveat = t_means)
    #cb = SavingCallback(OSET, saved_values; save_everystep=true)
    cf_vals = CFVals(rules, t_cors)
    LP = LProblem(A, cb, prob, rules, cf_vals, Tmax, tstep, delimiter)
    ####
    cf_vals, integrator = simulate(LP, 1, SpinLattice.RK4(); adaptive = false, dt = LP.tstep)

    return LP, cf_vals, integrator, sp
end

function test_clustered()
    for L in L1d[1:4]
        test_operators(Clustered(), L...)
        gc()
    end
    for L in L2d[1:4]
        test_operators(Clustered(), L...)
        gc()
    end
end

function test_all()
    Models = [Exact(), PureClassical(), Clustered(), Hybrid()]
    for model in Models
        for L in L1d[1:4]
            print(model,"\n")
            test_operators(model, L...)
            gc()
        end
    end
    for model in Models
        for L in L2d[1:4]
            print(model,"\n")
            test_operators(model, L...)
            gc()
        end
    end
    for model in Models[2:4]
        for L in L3d[2:4]
            print(model,"\n")
            test_operators(model, L...)
            gc()
        end
    end
end

function test_boundary()
    Models = [Exact(), PureClassical()]
    cbs = [true, false]
    for model in Models
        for L in L1d[1:4]
            for cb in cbs
                print(model, " cb=", cb,"\n")
                test_operators(model, L...,cb)
                gc()
            end
        end
    end
end

function test_3d()
    Models = [Exact(), PureClassical(), Clustered(), Hybrid()]
    for model in (Models[4],)
        for L in L3d[2:4]
            print(model,"\n")
            test_operators(model, L...)
            gc()
        end
    end
end

end


text:
"""
    SEDiamond{N}(axes, [dims]; [r])

A (holy) trait type for the N-dimensional diamond shape structuring element. This is a
special case of [`SEMask`](@ref SEMask) that ImageMorphology algorithms
might provide optimized implementation.

It is recommended to use [`strel_diamond`](@ref) and [`strel_type`](@ref):

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> se = strel_diamond((3, 3)) # C4 connectivity
3×3 SEDiamondArray{2, 2, UnitRange{$Int}, 0} with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_type(se)
SEDiamond{2, 2, UnitRange{$Int}}((-1:1, -1:1), (1, 2), 1)

julia> se = centered(collect(se)) # converted to normal centered array
3×3 OffsetArray(::Matrix{Bool}, -1:1, -1:1) with eltype Bool with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_type(se)
SEMask{2}()
```
"""
struct SEDiamond{N,K,R<:AbstractUnitRange{Int}} <: MorphologySE{N}
    axes::NTuple{N,R}
    dims::Dims{K}
    r::Int # radius
    function SEDiamond{N,K,R}(axes::NTuple{N,R}, dims::Dims{K}, r) where {N,K,R<:AbstractUnitRange{Int}}
        if !all(r -> first(r) == -last(r), axes)
            throw(ArgumentError("axes must be symmetric along each dimension"))
        end
        _is_unique_tuple(dims) || throw(ArgumentError("dims should be unique"))
        N >= K || throw(ArgumentError("`axes` length should be at least $K"))
        if !all(i -> i <= N, dims)
            throw(ArgumentError("all `dims` values should be less than or equal to $N"))
        end
        return new{N,K,R}(axes, dims, r)
    end
end
function SEDiamond{N}(ax::NTuple{N,R}, dims=ntuple(identity, N); r=maximum(length.(ax)) ÷ 2) where {N,R}
    return SEDiamond{N,length(dims),R}(ax, dims, r)
end

# Helper array type to build a consistant array interface for SEs
"""
    SEDiamondArray(se::SEDiamond)

The instantiated array object of [`SEDiamond`](@ref).
"""
struct SEDiamondArray{N,K,R<:AbstractUnitRange{Int},S} <: MorphologySEArray{N}
    axes::NTuple{N,R}
    dims::Dims{K}
    r::Int # radius
    _rdims::Dims{S}
end
function SEDiamondArray(se::SEDiamond{N,K,R}) where {N,K,R}
    _rdims = _cal_rdims(Val(N), se.dims)
    return SEDiamondArray{N,K,R,length(_rdims)}(se.axes, se.dims, se.r, _rdims)
end
strel_type(A::SEDiamondArray{N}) where {N} = SEDiamond{N}(A.axes, A.dims; r=A.r)
function strel_size(se::SEDiamondArray)
    return ntuple(i -> in(i, se.dims) ? 1 + 2 * se.r : 1, strel_ndims(se))
end

@inline Base.axes(A::SEDiamondArray) = A.axes
@inline Base.size(A::SEDiamondArray) = map(length, axes(A))
@inline Base.IndexStyle(::SEDiamondArray) = IndexCartesian()
Base.@propagate_inbounds function Base.getindex(A::SEDiamondArray{N,K}, inds::Int...) where {N,K}
    # for remaining dimensions, check if it is at the center position
    ri = _tuple_getindex(inds, A._rdims)
    all(iszero, ri) || return false
    # for masked dimensions, compare if the city-block distance to center is within radius
    mi = _tuple_getindex(inds, A.dims)
    r = isempty(mi) ? 0 : sum(abs, mi)
    return ifelse(r > A.r, false, true)
end

_tuple_getindex(t::Tuple, inds::Dims) = ntuple(i -> t[inds[i]], length(inds))
function _cal_rdims(::Val{N}, dims::NTuple{K}) where {N,K}
    return Dims{N - K}(filter(i -> !in(i, dims), 1:N))
end
_is_unique_tuple(t::Tuple) = any(i -> t[i] in t[1:(i - 1)], 2:length(t)) ? false : true

"""
    strel_diamond(A::AbstractArray, [dims]; r=1)
    strel_diamond(size, [dims]; [r])

Construct the N-dimensional structuring element (SE) for a diamond shape.

If image `A` is provided, then the SE size will be `(2r+1, 2r+1, ...)` with default
half-size `r=1`. If `size` is provided, the default `r` will be `maximum(size)÷2`. The
default `dims` will be all dimensions, that is, `(1, 2, ..., length(size))`.

```jldoctest; setup=:(using ImageMorphology; using ImageMorphology.StructuringElements)
julia> img = rand(64, 64);

julia> strel_diamond(img) # default size for image input is (3, 3)
3×3 SEDiamondArray{2, 2, UnitRange{$Int}, 0} with indices -1:1×-1:1:
 0  1  0
 1  1  1
 0  1  0

julia> strel_diamond(img; r=2) # equivalent to `strel_diamond((5,5))`
5×5 SEDiamondArray{2, 2, UnitRange{$Int}, 0} with indices -2:2×-2:2:
 0  0  1  0  0
 0  1  1  1  0
 1  1  1  1  1
 0  1  1  1  0
 0  0  1  0  0

julia> strel_diamond(img, (1,)) # mask along dimension 1
3×1 SEDiamondArray{2, 1, UnitRange{$Int}, 1} with indices -1:1×0:0:
 1
 1
 1

julia> strel_diamond((3,3), (1,)) # 3×3 mask along dimension 1
3×3 SEDiamondArray{2, 1, UnitRange{$Int}, 1} with indices -1:1×-1:1:
 0  1  0
 0  1  0
 0  1  0
```

!!! note "specialization and performance"
    The diamond shape `SEDiamond` is a special type for which many morphology algorithms may
    provide much more efficient implementations. For this reason, if one tries to collect an
    `SEDiamondArray` into other array types (e.g. `Array{Bool}` via `collect`), then a
    significant performance drop is very likely to occur.

See also [`strel`](@ref) and [`strel_box`](@ref).
"""
function strel_diamond(img::AbstractArray{T,N}, dims=coords_spatial(img); r::Union{Nothing,Int}=nothing) where {T,N}
    dims = _to_dims(Val(N), dims)
    sz, r = if isnothing(r)
        ntuple(i -> !isempty(dims) && in(i, dims) ? 3 : 1, N), 1
    else
        ntuple(i -> !isempty(dims) && in(i, dims) ? 2r + 1 : 1, N), r
    end
    return strel_diamond(sz, dims; r)
end
function strel_diamond(sz::Dims{N}, dims=ntuple(identity, N); kw...) where {N}
    dims = _to_dims(Val(N), dims)
    all(isodd, sz) || throw(ArgumentError("size should be odd integers"))
    ax = map(r -> (-r):r, sz .÷ 2)
    return SEDiamondArray(SEDiamond{N}(ax, dims; kw...))
end

# Tuple(1) is not inferable before Julia 1.9
# we want to support Colon input
@inline _to_dims(::Val{N}, i::Int) where {N} = (i,)
@inline _to_dims(::Val{N}, dims::Dims) where {N} = dims
@inline _to_dims(::Val{N}, ::Union{Colon,Tuple{Colon}}) where {N} = ntuple(identity, N)
@inline _to_dims(::Val{N}, v) where {N} = Tuple(v) # fallback


text:
"""
    SimpleEdgeIter

The function [`edges`](@ref) returns a `SimpleEdgeIter` for `AbstractSimpleGraphs`.
The iterates are in lexicographical order, smallest first. The iterator is valid for
one pass over the edges, and is invalidated by changes to the graph.

# Examples
```jldoctest
julia> using LightGraphs

julia> g = PathGraph(3);

julia> es = edges(g)
SimpleEdgeIter 2

julia> e_it = iterate(es)
(Edge 1 => 2, SimpleEdgeIterState [2, 2])

julia> iterate(es, e_it[2])
(Edge 2 => 3, SimpleEdgeIterState [0, 1])
```
"""
struct SimpleEdgeIter{G} <: AbstractEdgeIter
    g::G
end

struct SimpleEdgeIterState{T<:Integer}
    s::T  # src vertex or zero if done
    di::Int # index into adj of dest vertex
end

eltype(::Type{SimpleEdgeIter{SimpleGraph{T}}}) where {T} = SimpleGraphEdge{T}
eltype(::Type{SimpleEdgeIter{SimpleDiGraph{T}}}) where {T} = SimpleDiGraphEdge{T}

function edge_start(g::AbstractSimpleGraph{T}) where T
    s = one(T)
    @inbounds while s <= nv(g)
        isempty(fadj(g, s)) || return SimpleEdgeIterState(s, 1)
        s += one(T)
    end
    return SimpleEdgeIterState(zero(T), 1)
end

function edge_next(g::AbstractSimpleGraph,
    state::SimpleEdgeIterState{T}) where {T <: Integer}
    s = state.s
    di = state.di
    e = SimpleEdge(s, fadj(g, s)[di])
    di += 1
    @inbounds while s <= nv(g)
        sadj = fadj(g, s)
        while di <= length(sadj)
            if is_directed(g) || s <= sadj[di]
                return e, SimpleEdgeIterState(s, di)
            end
            di += 1
        end
        s += one(T)
        di = 1
    end
    return e, SimpleEdgeIterState(zero(T), 1)
end

function iterate(eit::SimpleEdgeIter{G}) where {G<:AbstractSimpleGraph}
    state = edge_start(eit.g)
    return iterate(eit, state)
end

function iterate(eit::SimpleEdgeIter{G}, state::SimpleEdgeIterState{T}) where {T,G<:AbstractSimpleGraph{T}}
    state.s == zero(T) && return nothing
    return edge_next(eit.g, state)
end

length(eit::SimpleEdgeIter) = ne(eit.g)

function _isequal(e1::SimpleEdgeIter, e2)
    k = 0
    for e in e2
        has_edge(e1.g, e) || return false
        k += 1
    end
    return k == ne(e1.g)
end
==(e1::SimpleEdgeIter, e2::AbstractVector{SimpleEdge}) = _isequal(e1, e2)
==(e1::AbstractVector{SimpleEdge}, e2::SimpleEdgeIter) = _isequal(e2, e1)
==(e1::SimpleEdgeIter, e2::Set{SimpleEdge}) = _isequal(e1, e2)
==(e1::Set{SimpleEdge}, e2::SimpleEdgeIter) = _isequal(e2, e1)

function ==(e1::SimpleEdgeIter, e2::SimpleEdgeIter)
    g = e1.g
    h = e2.g
    ne(g) == ne(h) || return false
    m = min(nv(g), nv(h))
    for i in 1:m
        fadj(g, i) == fadj(h, i) || return false
    end
    nv(g) == nv(h) && return true
    for i in m+1:nv(g)
        isempty(fadj(g, i)) || return false
    end
    for i in m+1:nv(h)
        isempty(fadj(h, i)) || return false
    end
    return true
end

in(e, es::SimpleEdgeIter) = has_edge(es.g, e)

show(io::IO, eit::SimpleEdgeIter) = write(io, "SimpleEdgeIter $(ne(eit.g))")
show(io::IO, s::SimpleEdgeIterState) = write(io, "SimpleEdgeIterState [$(s.s), $(s.di)]")