text:
println("===============================================")
println("Testing poynting-thomson.jl")
println("===============================================")

function Fract_PT_G_reduce_SLS()
    dt = 0.001
    # NaN at 0.0 due to InverseLaplace so start at dt
    t = collect(dt:dt:10.0)

    fracpt = relaxmod(Fract_PT, t, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = relaxmod(SLS_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanSLS[i]), eachindex(t))
end
@test Fract_PT_G_reduce_SLS()

function FractSLS_PT_G_reduce_SLS()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracSLS = relaxmod(FractSLS_PT, t, [2.0, 1.0, 1.0, 0.5])
    stanSLS = relaxmod(SLS_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(t))
end
@test FractSLS_PT_G_reduce_SLS()

function Fract_PT_J_reduce_SLS(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracpt = creepcomp(Fract_PT, t, [2.0, 0.999, 1.0, 0.0, 0.5, 0.0])
    stanSLS = creepcomp(SLS_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanSLS[i], atol=tol), eachindex(t))
end
@test Fract_PT_J_reduce_SLS(tol)

function FractSLS_PT_J_reduce_SLS(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracSLS = creepcomp(FractSLS_PT, t, [2.0, 0.999, 1.0, 0.5])
    stanSLS = creepcomp(SLS_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i], atol=tol), eachindex(t))
end
@test FractSLS_PT_J_reduce_SLS(tol)

function Fract_PT_Gp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracpt = storagemod(Fract_PT, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = storagemod(SLS_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanSLS[i]), eachindex(chirp))
end
@test Fract_PT_Gp_reduce_SLS()

function FractSLS_PT_Gp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracSLS = storagemod(FractSLS_PT, chirp, [2.0, 1.0, 1.0, 0.5])
    stanSLS = storagemod(SLS_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(chirp))
end
@test FractSLS_PT_Gp_reduce_SLS()

function Fract_PT_Gpp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracpt = lossmod(Fract_PT, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 0.0])
    stanSLS = lossmod(SLS_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanSLS[i]), eachindex(chirp))
end
@test Fract_PT_Gpp_reduce_SLS()

function FractSLS_PT_Gpp_reduce_SLS()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracSLS = lossmod(FractSLS_PT, chirp, [2.0, 1.0, 1.0, 0.5])
    stanSLS = lossmod(SLS_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracSLS[i], stanSLS[i]), eachindex(chirp))
end
@test FractSLS_PT_Gpp_reduce_SLS()

function Fract_PT_G_reduce_Jeffreys(tol)
    dt = 0.001
    # NaN at 0.0 due to Gamma function so start at dt
    t = collect(dt:dt:10.0)

    fracpt = relaxmod(Fract_PT, t, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = relaxmod(Jeffreys_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanJef[i], atol=tol), eachindex(t))
end
@test Fract_PT_G_reduce_Jeffreys(tol)

function FractJeffreys_PT_G_reduce_Jeffreys()
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracJef = relaxmod(FractJeffreys_PT, t, [2.0, 1.0, 0.0, 0.5])
    stanJef = relaxmod(Jeffreys_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(t))
end
@test FractJeffreys_PT_G_reduce_Jeffreys()

function Fract_PT_J_reduce_Jeffreys(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracpt = creepcomp(Fract_PT, t, [2.0, 0.999, 1.0, 0.0, 0.5, 1.0])
    stanJef = creepcomp(Jeffreys_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanJef[i], atol=tol), eachindex(t))
end
@test Fract_PT_J_reduce_Jeffreys(tol)

function FractJeffreys_PT_J_reduce_Jeffreys(tol)
    dt = 0.001
    t = collect(0.0:dt:10.0)

    fracJef = creepcomp(FractJeffreys_PT, t, [2.0, 1.0, 0.001, 0.5])
    stanJef = creepcomp(Jeffreys_PT, t, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i], atol=tol), eachindex(t))
end
@test FractJeffreys_PT_J_reduce_Jeffreys(tol)

function Fract_PT_Gp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracpt = storagemod(Fract_PT, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = storagemod(Jeffreys_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanJef[i]), eachindex(chirp))
end
@test Fract_PT_Gp_reduce_Jeffreys()

function FractJeffreys_PT_Gp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracJef = storagemod(FractJeffreys_PT, chirp, [2.0, 1.0, 0.0, 0.5])
    stanJef = storagemod(Jeffreys_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(chirp))
end
@test FractJeffreys_PT_Gp_reduce_Jeffreys()

function Fract_PT_Gpp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracpt = lossmod(Fract_PT, chirp, [2.0, 1.0, 1.0, 0.0, 0.5, 1.0])
    stanJef = lossmod(Jeffreys_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracpt[i], stanJef[i]), eachindex(chirp))
end
@test Fract_PT_Gpp_reduce_Jeffreys()

function FractJeffreys_PT_Gpp_reduce_Jeffreys()
    ω_step = 0.01
    chirp = collect(0.0:ω_step:10000.0)

    fracJef = lossmod(FractJeffreys_PT, chirp, [2.0, 1.0, 0.0, 0.5])
    stanJef = lossmod(Jeffreys_PT, chirp, [2.0, 1.0, 0.5])

    all(i -> isapprox(fracJef[i], stanJef[i]), eachindex(chirp))
end
@test FractJeffreys_PT_Gpp_reduce_Jeffreys()


text:
using Base.Test
using QuantLib

cal = QuantLib.Time.TargetCalendar()
settlementDate = Date(2002, 2, 19)
todays_date = Date(2002, 2, 15)
set_eval_date!(settings, todays_date)

flat_rate = Quote(0.04875825)
rhTermStructure = FlatForwardTermStructure(settlementDate, cal, flat_rate, QuantLib.Time.Actual365())

# Define swap
fixedLegFrequency = QuantLib.Time.Annual()
fixedLegConvention = QuantLib.Time.Unadjusted()
floatingLegConvention = QuantLib.Time.ModifiedFollowing()
fixedLegDayCounter = QuantLib.Time.EuroThirty360()
floatingLegFrequency = QuantLib.Time.Semiannual()

swapType = Payer()
dummyFixedRate = 0.03
indexSixMonths = euribor_index(QuantLib.Time.TenorPeriod(Dates.Month(6)), rhTermStructure)

startDate = QuantLib.Time.advance(Dates.Year(1), cal, settlementDate, floatingLegConvention)
maturity = QuantLib.Time.advance(Dates.Year(5), cal, startDate, floatingLegConvention)

fixedSchedule = QuantLib.Time.Schedule(startDate, maturity, QuantLib.Time.TenorPeriod(fixedLegFrequency), fixedLegConvention, fixedLegConvention, QuantLib.Time.DateGenerationForwards(), false, cal)
floatSchedule = QuantLib.Time.Schedule(startDate, maturity, QuantLib.Time.TenorPeriod(floatingLegFrequency), floatingLegConvention, floatingLegConvention, QuantLib.Time.DateGenerationForwards(), false, cal)

swap = VanillaSwap(swapType, 1000.0, fixedSchedule, dummyFixedRate, fixedLegDayCounter, indexSixMonths, 0.0, floatSchedule, indexSixMonths.dc, DiscountingSwapEngine(rhTermStructure))

# Calculations
@test fair_rate(swap) == 0.050000005910246725

# Other tests
@test QuantLib.get_fixed_reset_dates(swap) == Date[Date(2003, 2, 19), Date(2004, 2, 19), Date(2005, 2, 19), Date(2006, 2, 19), Date(2007, 2, 19)]
@test QuantLib.get_fixed_pay_dates(swap) == Date[Date(2004, 2, 19), Date(2005, 2, 21), Date(2006, 2, 20), Date(2007, 2, 19), Date(2008, 2, 19)]
@test QuantLib.get_floating_reset_dates(swap) == Date[Date(2003, 2, 19), Date(2003, 8, 19), Date(2004, 2, 19), Date(2004, 8, 19), Date(2005, 2, 21), Date(2005, 8, 19), Date(2006, 2, 20), Date(2006, 8, 21), Date(2007, 2, 19), Date(2007, 8, 20)]
@test QuantLib.get_floating_pay_dates(swap) == Date[Date(2003, 8, 19), Date(2004, 2, 19), Date(2004, 8, 19), Date(2005, 2, 21), Date(2005, 8, 19), Date(2006, 2, 20), Date(2006, 8, 21), Date(2007, 2, 19), Date(2007, 8, 20), Date(2008, 2, 19)]
@test QuantLib.get_floating_accrual_times(swap) == Float64[0.5027777777777778,0.5111111111111111,0.5055555555555555,0.5166666666666667,0.49722222222222223,0.5138888888888888,0.5055555555555555,0.5055555555555555,0.5055555555555555,0.5083333333333333]
@test QuantLib.get_floating_spreads(swap) == Float64[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

## CDS ##

ext: jl


text:
using Test

#--- test meshes
m1a = Mesh( (xmin=-2*π,xmax=π,N=3) )
m1b = Mesh( [-2*π -π 0 ] )

m2a = Mesh( (xmin=-2*π,xmax=2π,N=2) )
m2b = Mesh( [-2*π 0 ] )

function the_mesh(m::Mesh)
    [m.N, m.xmin, m.xmax, m.dx, m.x, m.kmin, m.kmax, m.dk, m.k]
end

@testset "Tests on Mesh" begin
    @test the_mesh(m1a) ≈ [3,-2π,π,π,[-2π; -π; 0],-2/3,2/3,2/3,[0; 2/3; -2/3]]
    @test the_mesh(m1a) == the_mesh(m1b)

    @test the_mesh(m2a) ≈ [2,-2π,2π,2π,[-2π; 0],-1//2,0,1//2,[0; -1/2]]
    @test the_mesh(m2a) == the_mesh(m2b)

end

#--- test times
t1a = Times( 0:0.1:1; ns=3 )
t1b = Times( (dt=0.1 , T=1); ns=3 )

t2a = Times( 0:1//10:1 ; Ns=3 )
t2b = Times( (dt=0.1 , T=1); Ns=3 )

function the_times(t::Times)
    [t.Nc, t.Ns, t.ns, t.tfin, t.dt, t.tc, t.ts]
end

@testset "Tests on times" begin
    @test the_times(t1a) == [11,4,[3,3,3],1.,0.1,0:0.1:1,[0.0, 0.3, 0.6, 0.9]]
    @test the_times(t1a) == the_times(t1b)

    @test the_times(t2a) == [11,4,[3,4,3],1,0.1,0:0.1:1,[0.0, 0.3, 0.7, 1.0]]
    @test the_times(t2a) == the_times(t2b)
end

#--- test data
v1 = [1; 3; 5];v2 = [2; 4; 6];
M = [ v1 v2 ];
data=Data(M)
@testset "Tests on data" begin
    @test data.U==[[1 2 ; 3 4 ; 5 6 ]]
    @test data.datalength==length(v1)
    @test data.datasize==2
end

#--- test problems
@testset "Tests on problems" begin
    param = (L=10, N=20, ϵ = 0.1, μ = 0.1, T=1, dt=1/30)
    init = Random(param)
    solver=Euler_naive()
    #solver=RK4(model1)  # the tests will fail with this solver, since it uses pre-allocation
    # the tests also fail if we define  model=WWn(param;n=1), and use model to build the problems
    pb1  = Problem(WWn(param;n=1,mesh=Mesh(param)), init, param; solver=solver)
    pb1a = Problem(WWn(param;n=1), init, param; solver=solver)
    pb1b = Problem(WWn(param;n=1), init, param; solver=solver)
    pb2  = Problem(WWn(param;n=2), init, merge(param,(ns=2,)); solver=solver)
    pb2a = Problem(WWn(param;n=2), init, merge(param,(ns=2,)); solver=solver)
    pb2b = Problem(WWn(param;n=2), init, merge(param,(ns=2,)); solver=solver)

    solve!([pb1a pb2a pb1b pb2b];verbose=false);
    solve!(pb1;verbose=false);solve!(pb2;verbose=false);

    for i in 1:2
        @test solution(pb1)[i]==solution(pb1a)[i]
        @test solution(pb1)[i]==solution(pb1b)[i]
        @test solution(pb2)[i]==solution(pb2a)[i]
        @test solution(pb2)[i]==solution(pb2b)[i]
    end
end

#--- test init
@testset "Tests on init" begin
    param = (L=10, N=2^8)
    mesh=Mesh(param);x=mesh.x;
    f(x) = exp.(-x.^2);g(x) = x.*exp.(-(x.-1).^2);
    inits=InitialData[]
    push!(inits, Init(f , g) )
    push!(inits, Init( mesh, f(x) , g(x); fast = false ) )
    push!(inits, Init( x, f(x) , g(x); fast = false ) )
    init = inits[1]
    for i in 2:3
        @test init.η([π])≈inits[i].η([π])
        @test init.v([π])≈inits[i].v([π])
    end
end


text:
using ApproxFun, Plots
using LinearAlgebra, SpecialFunctions
# using HCubature

θin = 0.4; X = 1.0; n = 3
num_coefs = 1000;
Y0 = 0.0; Y1 = 2.0;

d = Domain(0..2)

K2D(X,Y) = cos((Y+X*im)*sin(θin) * n*dot(Y,X))*sin(sqrt(X^2+Y^2+ X*Y*im))
K(Y) = cos((Y+X*im)*sin(θin) * n*dot(Y,X))*hankelh1(n,sqrt(X^2+Y^2+ X*Y*im))
K2Dreal(X,Y) = real(cos((Y+X*im)*sin(θin) * n*dot(Y,X))*sin(sqrt(X^2+Y^2+ X*Y*im)))
K2Dimag(X,Y) = imag(cos((Y+X*im)*sin(θin) * n*dot(Y,X))*sin(sqrt(X^2+Y^2+ X*Y*im)))

freal = Fun(K2Dreal, d^2, num_coefs,num_coefs)
fimag = ProductFun(K2Dimag, d^2, num_coefs,num_coefs)

f = Fun(K, Domain(0..2), num_coefs)
freal = Fun(Kreal,Domain(0..2), num_coefs)
fimag = Fun(Kreal,Domain(0..2), num_coefs)

function Bfull(n::Int,X; Y0= sqrt(k^2*a12^2-X^2), Y1 =1000000 , θin = 0.0)
    K(Y) = cos(Y*sin(θin) + n*atan(Y,X))*hankelh1(n,sqrt(X^2+Y^2))
    return 2.0*(-1.0)^n*sum(Fun(K,Y0..Y1))
end

function B(n::Int,X; Y0= sqrt(k^2*a12^2-X^2), θin = 0.0)
    Y1 = max(2000.0*X, 4000.0)
    # assymptotically approximate the integral from Y1 to Inf
    Binf = (1.0+1.0im)*exp(im*Y1*(1.0 - sin(θin)))*
        (1.0 + (-1.0)^n*exp(2.0im*Y1*sin(θin))*(1.0 - sin(θin)) + sin(θin))/(sqrt(pi*Y1)*cos(θin)^2)

    return Binf + Bfull(n, X ;Y0=Y0, Y1=Y1, θin=θin)
end

θin = 0.1
a12 =1.0
n=1
k=1.
X = 0.2
K(Y) = cos(Y*sin(θin) + n*atan(Y,X))*hankelh1(n,sqrt(X^2+Y^2))

interval = sqrt(k^2*a12^2-X^2)..1000000
interval1 = interval.left..1000.7
interval2 = interval1.right..interval.right
Ys = Fun(identity,interval)
Y1s = Fun(identity,interval1)
Y2s = Fun(identity,interval2)
Ks = K.(Ys)

interval1 = 0.0..1.0
B_app =  Fun(x -> B(0,x),interval1)

sum(Ks) - sum(K.(Y1s)) - sum(K.(Y2s))
@time 2.0*(-1.0)^n*sum(Ks)
@time 2.0*(-1.0)^n.*hcubature(Y -> K(Y[1]), (interval.left,), (interval.right,))

Σ = DefiniteIntegral(Chebyshev(interval))
@time sum(Σ[i]*Ks.coefficients[i] for i in eachindex(Ks.coefficients))

# Solving fredholm https://github.com/JuliaApproximation/ApproxFun.jl/issues/570
d = Interval(0.0, 2.0)
K2 = Fun((x,y) -> exp(-(x-y)), d^2)
V = DefiniteIntegral(Chebyshev(d))[LowRankFun(K2,d^2)]
one = Fun(x->1.0, d)
u = (I - V)\one

# Test
x = Fun(d)
u(0.1)-cumsum(exp(-(0.1-x))*u)(2.0) # ≈ 1

# Solving simple fredholm http://juliaapproximation.github.io/ApproxFun.jl/latest/usage/equations.html#Linear-equations-1
interval = 0.0..10.0
α=1.0
x = Fun(interval)
Σ = DefiniteIntegral(Chebyshev(interval))

b = exp((im-α)*x)-1.0im/(im-α)
u = (I+1.0im*Σ)\b;

#Analytically the solution should be
norm(u.(0.0:0.1:5.0)-exp((im-α).*(0.0:0.1:5.0)))/norm(exp((im-α).*(0.0:0.1:5.0)))

@time sum(Σ[i]*Ks.coefficients[i] for i in eachindex(Ks.coefficients))

a =1.0
θin = 0.1
n=1
k=1.
X = 0.2

θin = 0.3

d = Interval(0.0, 10.0)
xs = Fun(d)
test_fun = exp(-xs)
S_fun = S(0,xs; θin = 0.3)
sum(S_fun*test_fun) ≈ 2*(exp(10 .* (im*cos(θin)-1))-1.0)/(cos(θin)*(im*cos(θin)-1.0))

Sop = DefiniteIntegral(Chebyshev(d))[S_fun]
Sop*test_fun ≈ 2*(exp(10 .* (im*cos(θin)-1))-1.0)/(cos(θin)*(im*cos(θin)-1.0))

x1 = Fun(Interval(0.0, 1.0))
@time f1 = Fun( x-> exp(-x), 0.0..1.0)
@time f3 = Fun( x-> (x>1.0) ? 0.0 : exp(-x), 0.0..10.0)
f2 = exp(-xs)

f4 = f1+f2

f1(0.3)+f2(0.3)
2*f4(0.3)

f1(1.3)+f2(1.3)
f4(1.3)

dd_fun = Fun((x,X) -> (abs(X-x) < 1.0) ? S_fun(0,X-x) :0.0+0.0im, d^2,10)

dd_fun = Fun((x,X) -> S_fun(X), d^2,10)

d1 = 0.0..1.0
B_fun = Fun(x -> integrate_B(0,x, sqrt(1.0 - x^2)), d1,40)#450)
S_fun = S(0,xs)

BS_lowfun = LowRankFun((x,X) -> (abs(X-x) < 1.0) ? B_fun(X-x) - S_fun(X-x) : 0.0+0.0im, d^2)
BS_lowfun = LowRankFun((x,X) -> (abs(X-x) < 1.0) ? B_fun(X-x) - S_fun(X-x) : 0.0+0.0im, d^2)
BS_fun = Fun((x,X) -> (abs(X-x) < 1.0) ? B_fun(X-x) - S_fun(X-x) : 0.0+0.0im, d^2,1000)

# r1 = 10.0*rand(10);
# r2 = 10.0*rand(10);
#
# m = mean( abs((abs(X-x) < 1.0) ? B_fun(X-x) - S_fun(X-x) : 0.0+0.0im)  for x in r1, X in r2)
# [ abs( BS_fun(x,X) - ((abs(X-x) < 1.0) ? B_fun(X-x) - S_fun(X-x) : 0.0+0.0im))  for x in r1, X in r2]
# [BS_fun(x,X) for x in r1, X in r2]

# Y0 = sqrt(k^a12^2 - (X-x)^2)
BS_op = DefiniteIntegral(Legendre(d))[
    BS_fun
]
# Y0 = sqrt(k^a12^2 - (X-x)^2)
BS_op = DefiniteIntegral(Chebyshev(d))[
    LowRankFun((x,X) -> (abs(X-x) < 1.0) ? B_fun(X-x) - S_fun(X-x) : 0.0+0.0im, d^2)
]

d = Interval(0.0, 2.0)
fun = Fun((x,X) -> (abs(X-x) < 1.0) ? exp(-abs(X-x)) : 0.0+0.0im, d^2,1000)
DefiniteIntegral(Chebyshev(d))[fun]
DefiniteIntegral(Legendre(d))[Fun((x,X) -> (abs(X-x) < 1.0) ? exp(-abs(X-x)) : 0.0+0.0im, d^2,1000)]
DefiniteIntegral(Chebyshev(d))[Fun((x,X) -> (abs(X-x) < 1.0) ? exp(-abs(X-x)) : 0.0+0.0im, d^2,1000)]
DefiniteIntegral(Chebyshev(d))[LowRankFun((x,X) -> (abs(X-x) < 1.0) ? exp(-abs(X-x)) : 0.0+0.0im, d^2)]
DefiniteIntegral(Chebyshev(d))[LowRankFun((x,X) -> exp(-abs(X-x)), d^2)]

f = LowRankFun((x,X) -> exp(-abs(X-x)))
f1 = Fun(x -> exp(-abs(x)),d,100)
f = Fun((x,X) -> exp(-abs(X-x)),d^2,100)

int_op = DefiniteIntegral(Legendre(d))[
LowRankFun((x,X) -> (abs(X-x) < 1.0) ? cos(pi*(X-x)/2.0): 0.0+0.0im, d^2)
]
BS_op = DefiniteIntegral(Legendre(d))[
Fun((x,X) -> cos(pi*(X-x)/2.0), d^2,100)
]

max_stars_repo_name: arturgower/EffectiveWaves


text:
using GAP_pkg_io
using Test

import GAP_pkg_io
println("TODO")


text:
# Watts.jl

###### watts #####
"""
    watts(v, k)

    Compute the Watts Poverty Index of a vector `v` at a specified absolute
    poverty line `k`.

# Examples
```julia
julia> using Inequality
julia> watts([8, 5, 1, 3, 5, 6, 7, 6, 3], 4)
0.217962056224828
```
"""
watts(v::AbstractVector{<:Real}, k::Real)::Float64 = length(v[v .< k]) == 0 ?  0. : sum(log.(k./v[v .< k]))/length(v)

###### weighted watts #####
"""
    watts(v, w, α)

Compute the Watts Poverty Index of a vector `v` at a specified absolute
poverty line `α`, using weights given by a weight vector `w`.

Weights must not be negative, missing or NaN. The weights and data vectors must have the same length.

# Examples
```julia
julia> using Inequality
julia> watts([8, 5, 1, 3, 5, 6, 7, 6, 3], collect(0.1:0.1:0.9), 4)
0.17552777833850716
```
"""
function watts(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, α::Real)::Float64

    checks_weights(v, w)

    length(v) == 0 ? (return 0) : nothing

    return sum( (log.(α./(v[v .< α]))) .* w[v .< α] ) / sum(w)

end

function watts(v::AbstractVector{<:Real}, w::AbstractWeights, α::Real)::Float64

    checks_weights(v, w)

    length(v) == 0 ? (return 0) : nothing

    return sum( (log.(α./(v[v .< α]))) .* w[v .< α] ) / sum(w)

end

"""
    wwatts(v, w, α)

Compute the Watts Poverty Index of `v` with weights `w`. See also [`watts`](@watts)
"""
wwatts(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, α::Real)::Float64 = watts(v, w, α)

wwatts(v::AbstractVector{<:Real}, w::AbstractWeights, α::Real)::Float64 = watts(v, w, α)


text:
"""
Julia SNN Implementation of AdExp Neuron.
[Adaptive_exponential_integrate and fire neuron](http://www.scholarpedia.org/article/Adaptive_exponential_integrate-and-fire_model)
Dr. Wulfram Gerstner
Romain Brette, Ecole Normale Supérieure, Paris, France
"""

@snn_kw struct ADEXParameter{FT=Float32}
    a::FT = 4.0
    b::FT = 0.0805
    cm::FT = 0.281
    v_rest::FT = -70.6
    tau_m::FT = 9.3667
    tau_w::FT = 144.0
    v_thresh::FT = -50.4
    delta_T::FT = 2.0
    v_reset::FT = -70.6
    spike_height::FT = 30
end
@snn_kw mutable struct AD{VFT=Vector{Float32},VBT=Vector{Bool}}
    #    param::ADEXParameter = ADEXParameter()
    #end
    #@with_kw mutable struct AD
    param::ADEXParameter = ADEXParameter(a,
                                        b,
                                        cm,
                                        v_rest,
                                        tau_m,
                                        tau_w,
                                        v_thresh,
                                        delta_T,
                                        v_reset,
                                        spike_height)
    N::SNNInt = 1
    cnt::SNNInt = 2
    v::Vector{SNNFloat} = fill(param.v_rest, N)
    w::Vector{SNNFloat} = zeros(N)
    fire::Vector{Bool} = zeros(Bool, N)
    I::Vector{SNNFloat} = zeros(N)
    spike_raster::Vector{SNNInt} = zeros(N)
    records::Dict = Dict()
end

function integrate!(p::AD, param::ADEXParameter, dt::SNNFloat)
    @unpack N, cnt, v, w, fire, I,spike_raster = p
    @unpack a,b,cm,v_rest,tau_m,tau_w,v_thresh,delta_T,v_reset,spike_height = param
    if spike_raster[cnt-1] == 1 || fire[1]
      v[1] = v_reset
      w[1] += b
    end
    dv  = (((v_rest-v[1]) +
            delta_T*exp((v[1] - v_thresh)/delta_T))/tau_m +
            (I[1] - w[1])/cm) *dt
    v[1] += dv
    w[1] += dt * (a*(v[1] - v_rest) - w[1])/tau_w * dt

    fire[1] = v[1] > v_thresh

    if v[1]>v_thresh
        fire[1] = 1 # v[1] > vPeak
        v[1] = spike_height
        spike_raster[cnt] = 1

    else
        spike_raster[cnt] = 0
    end

    cnt+=1
end


text:
include("utils.jl")

""" Nominal facility location model. """
function facility_model(c::Matrix, f::Vector)
    n, m = size(c)
    @assert length(f) == n
    model = Model(Gurobi.Optimizer)

    # VARIABLES
    @variable(model, x[1:n], Bin)      # Facility locations
    @variable(model, y[1:n, 1:m] >= 0) # Flow variables
    # CONSTRAINTS
    @constraint(model, [j=1:m], sum(y[i, j] for i = 1:n) >= d[j])      # demand constraint
    @constraint(model, [i=1:n], sum(y[i, j] for j = 1:m) <= s[i]*x[i]) # capacity constraint
    # OBJECTIVE
    @objective(model, Min,
            sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * y[i, j] for i=1:n, j=1:m))
    return model, x, y
end

model, x, y = facility_model(c, f)

optimize!(model)

plt = plot_solution(model, x, y)


text: include("foo.jl")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@testset "warn-pv" begin
    resource()
    d = F.PageVars(
    	"a" => 0.5 => (Real,),
    	"b" => "hello" => (String, Nothing))
    s = @capture_out F.set_vars!(d, ["a"=>"\"blah\""])
    @test occursin("Franklin Warning: in <unknown>", s)
    @test occursin("Page var 'a' (type(s): (Real,)) cannot", s)
    @test occursin("assignment will be ignored",  s)
end

@testset "warn-config" begin
    resource(); gotd()
    s = @capture_out F.process_config()
    @test occursin("Warning: in <config.md>", s)
    @test occursin("No 'config.md' file found", s)
end

@testset "warn-rss" begin
    resource(); gotd()
    empty!(F.RSS_DICT)
    F.set_var!(F.GLOBAL_VARS, "website_title", "Website title")
    F.set_var!(F.GLOBAL_VARS, "website_descr", "Website descr")
    F.set_var!(F.GLOBAL_VARS, "website_url", "https://github.com/tlienart/Franklin.jl/")
    F.def_LOCAL_VARS!()
    set_curpath("hey/ho.md")
    F.set_var!(F.LOCAL_VARS, "rss_title", "")
    s = @capture_out F.add_rss_item()
    @test occursin("Warning: in <input string>", s)
    @test occursin("An RSS description was found but without title", s)
    @test occursin("for page '/hey/ho/index.html'.", s)
end

@testset "warn-paginate" begin
    resource(); gotd()
    write(joinpath(td, "config.md"), "@def aa = 5")
    mkpath(joinpath(td, "_css"))
    mkpath(joinpath(td, "_layout"))
    write(joinpath(td, "_layout", "head.html"), "HEAD\n")
    write(joinpath(td, "_layout", "foot.html"), "\nFOOT\n")
    write(joinpath(td, "_layout", "page_foot.html"), "\nPG_FOOT\n")
    write(joinpath(td, "index.md"), "INDEX")
    write(joinpath(td, "foo.md"), raw"""
        @def a = ["
Item $i
" for i in 1:10]
        Some content
        {{paginate abc 4}}
        """)
    s = @capture_out serve(single=true)
    @test occursin("Franklin Warning: in <foo.md>", s)
    @test occursin("A h-function call '{{paginate ...}}' has some", s)
    @test occursin("The page variable 'abc' does not match", s)

    write(joinpath(td, "foo.md"), raw"""
        @def a = ["
Item $i
" for i in 1:10]
        Some content
        {{paginate a iehva}}
        """)
    s = @capture_out serve(single=true)
    @test occursin("Failed to parse", s)
    @test occursin("(given: 'iehva')", s)
    @test occursin("Setting to", s)

    write(joinpath(td, "foo.md"), raw"""
        @def a = ["
Item $i
" for i in 1:10]
        Some content
        {{paginate a -5}}
        """)
    s = @capture_out serve(single=true)
    @test occursin("Non-positive number", s)
    @test occursin("('-5' read from -5)", s)
    @test occursin("Setting to", s)

    write(joinpath(td, "foo.md"), raw"""
        @def a = ["
Item $i
" for i in 1:10]
        Some content
        ~~~
~~~
        {{paginate a 4}}
        ~~~
~~~
        ~~~
~~~
        {{paginate a 4}}
        ~~~
~~~
        """)
    s = @capture_out serve(single=true)
    @test occursin("Multiple calls to", s)
end


text:
module TestDataFrames

using Test
using DataFrames
using DataFramesMeta

df = DataFrame(A = 1:3, B = [2, 1, 2])

x = [2, 1, 0]

@test  @with(df, :A .+ 1)   ==  df[:A] .+ 1
@test  @with(df, :A .+ :B)  ==  df[:A] .+ df[:B]
@test  @with(df, :A .+ x)   ==  df[:A] .+ x

x = @with df begin
    res = 0.0
    for i in 1:length(:A)
        res += :A[i] * :B[i]
    end
    res
end
idx = :A
@test  @with(df, cols(idx) .+ :B)  ==  df[:A] .+ df[:B]
idx2 = :B
@test  @with(df, cols(idx) .+ cols(idx2))  ==  df[:A] .+ df[:B]

@test  x == sum(df[:A] .* df[:B])
@test  @with(df, df[:A .> 1, ^([:B, :A])]) == df[df[:A] .> 1, [:B, :A]]
@test  @with(df, DataFrame(a = :A * 2, b = :A .+ :B)) == DataFrame(a = df[:A] * 2, b = df[:A] .+ df[:B])

@test DataFramesMeta.where(df, 1) == df[1, :]

@test  @where(df, :A .> 1)          == df[df[:A] .> 1,:]
@test  @where(df, :B .> 1)          == df[df[:B] .> 1,:]
@test  @where(df, :A .> x)          == df[df[:A] .> x,:]
@test  @where(df, :B .> x)          == df[df[:B] .> x,:]
@test  @where(df, :A .> :B)         == df[df[:A] .> df[:B],:]
@test  @where(df, :A .> 1, :B .> 1) == df[map(&, df[:A] .> 1, df[:B] .> 1),:]
@test  @where(df, :A .> 1, :A .< 4, :B .> 1) == df[map(&, df[:A] .> 1, df[:A] .< 4, df[:B] .> 1),:]

@test DataFramesMeta.select(df, :A) == df[:A]

@test DataFramesMeta.orderby(df, df[[1, 3, 2], :]) == df[[1, 3, 2], :]

@test @byrow!(df, if :A > :B; :A = 0 end) == DataFrame(A = [1, 0, 0], B = [2, 1, 2])
@test  df == DataFrame(A = [1, 0, 0], B = [2, 1, 2])

df = DataFrame(A = 1:3, B = [2, 1, 2])  # Restore df
y = 0
@byrow!(df, if :A + :B == 3; global y += 1 end)
@test  y == 2

df = DataFrame(A = 1:3, B = [2, 1, 2])
df2 = @byrow! df begin
    @newcol colX::Array{Float64}
    @newcol colY::Array{Float64}
    :colX = :B == 2 ? pi * :A : :B
    if :A > 1
        :colY = :A * :B
    end
end

@test  df2[:colX] == [pi, 1.0, 3pi]
@test  df2[2, :colY] == 2

end # module


text:
# Install (Julia 1.0.2)
using Pkg
pkg"activate ."
pkg"instantiate"
# Gurobi (requires license)
#pkg"add Gurobi"
pkg"precompile"

max_stars_repo_name: martinbiel/mpc-reproducibility

max_issues_repo_name: martinbiel/mpc-reproducibility

max_forks_repo_name: martinbiel/mpc-reproducibility


text:
# title: neten.jl
# author: twab
# description: julia implementation of network enhancement algorithm
# license: MIT

module neten

## ---- imports

using RData
using LinearAlgebra

## ---- functions

function hello()
		# a function for testing
		println("Hello World!")
end

function data(rdata::String)
		# load package data
		# using RData
		path = joinpath(dirname(pathof(neten)),"data","$rdata.rda")
		return load(path)[rdata]
end

#function load_rdata(rdata::String, package::String)
#	# FIXME: why not have function data(String)?
#	# load any data object associated with an R package
#	# NOTE: this doesn't work as expected bc most package data is not saved
#	# internally as .rda
#	# construct path to R data.rda file in data/
#	file = joinpath("data","$rdata.rda")
#	# construct R command
#	cmd = "system.file('$file', package='$package')"
#	# parse output of R command
#	out = read(`Rscript -e $cmd`,String)
#	path = strip(split(out)[2],'\"')
#	return load(path)[rdata]
#end

function colSums(a::Array{Float64,2}, dropdim::Bool=true)
	# return column sums of an array
	col_sums = sum(a, dims=1)
	if dropdim
		col_sums = dropdims(col_sums, dims=1)
	end
	return col_sums
end

function rowSums(a::Array{Float64,2}, dropdim::Bool=true)
	# return row sums of an array
	row_sums = sum(a, dims=2)
	if dropdim
		row_sums = dropdims(row_sums, dims=2)
	end
	return row_sums
end

function NE_dn(Win::Array{Float64,2}, k::Int64=size(Win,1))
	# what is this function doing?
	# NE = network enhancment?, dn == ??
	#k = size(Win, 1)
	Wk = Win * k
	d = rowSums(abs.(Wk))
	D = diagm(1 ./ d)
	Wne = D * Wk
	return Wne
end

function dominateset(Win::Array{Float64,2}, k::Int64=minimum([20,size(Win,1)-1]))
	# dominating set is set of connected nodes
	W0 = abs.(Win)
	#k = minimum([k, size(W0, 1)-1])
	# sort each row in decreasing order
	# do this to each row!
	A = sort(W0, dims=2, rev=true)
	B = mapslices(x -> sortperm(x, rev=true),
		      W0; dims=2)
	res = A[:,1:k]
	x = collect(1:size(W0, 1))
	inds=repeat(x, outer = [1, k])
	loc = B[:, 1:k]
	P0 = zeros(size(W0,2), size(W0,1))
	idx = (vec(loc .- 1)) * size(W0,2) + vec(inds)
	P0[idx] = res
	P = (P0 + transpose(P0)) / 2
	return P .* sign.(Win)
end

function transitionFields(Win::Array{Float64,2})
	zeroindex = rowSums(Win).==0
	W1 = Win * size(Win,1)
	W2 = NE_dn(W1)
	# do this operation on each col / sqrt(sum(abs(col))
	W3 = mapslices(x -> x./sqrt(sum(abs.(x)) + eps()),
		       W2; dims=1)
	# compute cross product: x %*% t(x) (tcrossprod)
	W4 = W3 * transpose(W3)
	W4[zeroindex,:] .= 0
	W4[:,zeroindex] .= 0
	return W4
end

function enhance(Win::Array{Float64,2}, k::Int64=20, alpha::Float64=0.9, diffusion::Float64=2.0)
	# network enhancement diffusion process is given by:
	# W_{t+1} = αΓ x W_{t} x Γ + Γ(1 - α) (eq 2)
	# set diag == 0
	Win[diagind(Win)] .= 0;
	# remove unconnected nodes
	nonzero = colSums(Win) .> 0;
	W0 = Win[nonzero, nonzero];
	# do NE_dn
	Wne = NE_dn(W0);
	# mean? or some sort of normalization
	Wnorm = (Wne + transpose(Wne)) ./ 2;
	# calc dominating set
	P = dominateset(Wnorm, k);
	# working with P...
	I = diagm(ones(size(P,1))); # identiy matrix
	X = diagm(rowSums(abs.(P))); # diagonal = rowSums
	Pnew = P .+ I .+ X;
	# do transitionFields
	Pt = transitionFields(Pnew)
	# compute eigen vectors and values
	e = eigen(Pt)
	d = e.values[end:-1:1] .- eps() # rev order to match Rout
	U = -1 .* e.vectors[:, end:-1:1] # rev order and flip sign to match R output
	# generate diag matrix D from difusion process
	dd = (1-alpha) * d ./ (1 .- alpha * d.^diffusion) # <- diffusion process
	D = diagm(dd)
	W2 = U * D * transpose(U)
	W3 = W2 .* (1 .- diagm(ones(size(W2,1)))) ./ (1 .- diag(W2))
	# DD is is the diagonal matrix whose entries are the degree of the vertices
	DD = diagm(colSums(abs.(W0)))
	W4 = DD * W3
	# remove negative
	W4[W4 .< 0] .= 0
	W5 = (W4 + transpose(W4)) / 2
	Wout = zeros(size(Win))
	Wout[nonzero, nonzero]  = W5
	# set diag to zero
	#Wout[diagind(Wout)] .= 0
	return Wout
end

## ---- exports

export hello
export data
export enhance

end


text:
using TiledArrays
using Test

@testset "TiledArrays.jl" begin

@testset "BufferedArrays" begin
    using TiledArrays.BufferedArrays

    data = reshape(1:240, 4, 10, 6)
    arr = BufferedArray(data, (2, 5, 3))

    @test size(arr) == (4, 10, 6)
    @inferred size(arr)

    @test IndexStyle(typeof(arr)) === IndexTiled()

    @test arr[2, 2, 1] == 6

    I = rand.(axes(data))
    @test arr[I...] = data[I...]

    @test arr[92] == 92

    i = rand(eachindex(data))
    @test arr[i] == data[i]
end

end


text:
function MersenneTwister()
    return MersenneTwister(())
end

function MersenneTwister(arg0::Vector{jint})
    return MersenneTwister((Vector{jint},), arg0)
end

function MersenneTwister(arg0::jint)
    return MersenneTwister((jint,), arg0)
end

function MersenneTwister(arg0::jlong)
    return MersenneTwister((jlong,), arg0)
end

function next_int(obj::MersenneTwister)
    return jcall(obj, "nextInt", jint, ())
end

function set_seed(obj::MersenneTwister, arg0::Vector{jint})
    return jcall(obj, "setSeed", void, (Vector{jint},), arg0)
end

function set_seed(obj::MersenneTwister, arg0::jint)
    return jcall(obj, "setSeed", void, (jint,), arg0)
end


text:
__precompile__

module QWNNM

    using Images
    using Random
    using Statistics
    using LinearAlgebra
    using Quaternions
    using GenericSVD

    include("addNoise.jl")
    include("psnr.jl")
    include("denoising.jl")

end


text:
function collect_mers(::Type{M}, count_mode::CountMode, input::DatastoreBuffer{<:PairedReads}, range::AbstractRange = 1:length(input)) where {M<:AbstractMer}
    max_read_size = max_read_length(ReadDatastores.datastore(input))
    v = Vector{M}(undef, length(input) * (max_read_size - BioSequences.ksize(M) + 1))
    wi = firstindex(v)
    read_sequence = eltype(input)()
    @inbounds for i in range
        for mer in each(M, load_sequence!(input, i, read_sequence))
            v[wi] = count_mode(mer)
            wi = wi + 1
        end
    end
    return resize!(v, wi - 1)
end

"""
    serial_mem(::Type{M}, input::DatastoreBuffer{<:ReadDatastore}, count_mode::CountMode, range::AbstractRange = 1:length(input)) where {M<:AbstractMer}

MerCounting's simplest kmer counting method.

Build a sorted list (vector) of kmer counts (MerCount), serially and entirely in memory.

!!! warning
    This function is a serial and in memory `MerCount` list builder that can build a
    kmer count from a ReadsDatastore on its own (if you have memory and time),
    but it is also intended to be composed into other multi-process or multi-threaded
    kmer counting strategies.

    This method estimates roughly how many kmers will be generated by the reads
    specified by `range` in the dataset. It then pre-allocates an array to contain
    them. It then collects the kmers, sorts, them, and then collapses them into a
    list of counts sorted by the kmer.

    So if you want to count kmers and have the resources to throw at it, this is
    the simplest method, and possibly even the quickest given that simplicity.
"""
function serial_mem(::Type{M}, input::DatastoreBuffer{<:ReadDatastore}, count_mode::CountMode, range::AbstractRange = 1:length(input)) where {M<:AbstractMer}
    @info "[Process: $(myid())] Collecting kmers from $(length(range)) reads in read datastore $(name(ReadDatastores.datastore(input)))"
    all_mers = collect_mers(M, count_mode, input, range)
    @info "[Process: $(myid())] Collapsing collected kmers into counts"
    return collapse_into_counts(all_mers)
end

macro serial_mem(K::Int, count_mode::Symbol, filename::String, range::Expr)
    # Figure out type of read datastore
    dstype = ReadDatastores.deduce_datastore_type(filename)
    # Process the range argument
    @assert range.head === :call
    @assert range.args[1] === :(:)
    interval = range.args[2]:range.args[3]
    quote
        open($dstype, $filename) do ds
            serial_mem(DNAMer{$K}, $count_mode, ds, $interval)
        end
    end
end


text:
using PanelMatrices
using BenchmarkTools
using StaticArrays
using StaticNumbers
using LinearAlgebra
using Statistics

# Make sure we're benchmarking BLAS single-threaded
LinearAlgebra.BLAS.set_num_threads(1)

# In case it affects benchmarking
set_zero_subnormals(true)

const y = reshape(collect(1:100), 10, 10)
const x = PanelMatrix(y)
const v = @SMatrix [1 2 3 4; 5 6 7 8; 9 10 11 12; 13 14 15 16]

println("Benchmarking using Julia", VERSION)

let
    println("Setting each individual element")
    function f(x, s)
        for i ∈ eachindex(x)
            @inbounds x[i] = s
        end
    end
    @btime $f($x, 0)
end

let
    println("Setting each panel")
    function f(x, v)
        for i ∈ 1:3, j ∈ 1:3
            @inbounds PanelMatrices.set_full_panel!(x, v, i, j)
        end
    end
    @btime $f($x, $v)
end

for T in (Float32, Float64), M = 7:5:32, N = 7:5:32, K = 7:5:32
    if M*N*K <= 25*25*25
        println("\n$T $M x $N x $K mul!")

        A0 = reshape(collect(T(1):T(M*K)), M, K)
        B0 = reshape(collect(T(1):T(K*N)), K, N)
        C0 = A0*B0

        A = PanelMatrix(A0)
        B = PanelMatrix(B0)
        C = PanelMatrix{T}(undef, (M,N))

        mul!(C, A, B)
        @assert C ≈ C0

        tb = @benchmark mul!($C, $A, $B) samples=10 evals=1000
        println(2*M*N*K/minimum(tb.times), " Gflops")

        println("Julia arrays using $(LinearAlgebra.BLAS.vendor()) (single thread), for comparison")
        tbr = @benchmark mul!($C0, $A0, $B0) samples=100 evals=100
        println(2*M*N*K/minimum(tbr.times), " Gflops")
    end
end


text:
############################################################################################
"Mostly generated functions that allow to transform NamedTuples without much allocations. Not Exported."
#!NOTE: @generated does not seem to work properly with DocStringExtensions.jl, so documentation is done manually.
############################################################################################
"""
$(SIGNATURES)
Subset `nt` with arguments `s`

# Examples
```julia
julia> subset( (a = 1., b = 2., c = 3.), (:a, :b) )
(a = 1., b = 2.)
```

"""
@inline subset(nt::NamedTuple, s::Tuple{Vararg{Symbol}}) = NamedTuple{s}(nt)
@inline subset(nt::NamedTuple, s::Symbol) = NamedTuple{(s,)}(nt)

############################################################################################
#=
"""
$(SIGNATURES)
Convert Tuple tup to NamedTuple with all fields equal to val.

# Examples
```julia
```

"""
=#
"Convert Tuple `tup` to NamedTuple with all fields equal to `val`."
function Tuple_to_Namedtuple(tup::Tuple, val)
    return NamedTuple{tup}(map(x -> val, eachindex(tup)))
end

#=
"""
$(SIGNATURES)
Convert Tuple tup to NamedTuple with fields equal to val.

# Examples
```julia
```

"""
=#
"Convert Tuple `tup` to NamedTuple with fields equal to `val`."
function to_NamedTuple(tup::Tuple, val::AbstractVector)
    if length(tup) == 1
        #!NOTE Safe method even if Tuple is only a single symbol and val non-scalar
        return NamedTuple{tup}(map(x -> val, eachindex(tup)))
    else
        return NamedTuple{tup}(map(iter -> val[iter], eachindex(val)))
    end
end
############################################################################################
#=
"""
$(SIGNATURES)
Subset NamedTuple obj with keys of s without allocations.

# Examples
```julia
```

"""
=#
"Subset NamedTuple `obj` with keys of `s` without allocations."
@generated function subset(obj::NamedTuple, s::NamedTuple)
    #!NOTE: https://discourse.julialang.org/t/generated-function-iterate-over-function-argument/66087/10
    fnames = fieldnames(s)
    fvals = map(fnames) do fname
        Expr(:call, :getproperty, :obj, QuoteNode(fname))
    end
    fvals = Expr(:tuple, fvals...)
    return :(NamedTuple{$fnames}($fvals))
end

############################################################################################
"Convert strcut `x` to tuple."
@generated function to_Tuple_generated(x)
    #!NOTE: see https://discourse.julialang.org/t/get-fieldnames-and-values-of-struct-as-namedtuple/8991/3
    tup = Expr(:tuple)
    for i in 1:fieldcount(x)
        push!(tup.args, :(getfield(x, $i)))
    end
    return :($tup)
end

############################################################################################
#=
#=
"""
$(SIGNATURES)
Subset NamedTuple x given symbols sym without additional allocations.

# Examples
```julia
```

"""
=#
"Subset NamedTuple `x` given symbols `sym` without additional allocations."
@generated function to_Tuple_generated(x, sym...)
    tup = Expr(:tuple)
    for v in sym
        push!(tup.args, :(getfield(x, getval($v))))
    end
    return :($tup)
end
getval(::Type{Val{T}}) where {T} = T
=#

############################################################################################
#=
"""
$(SIGNATURES)
Convert struct x to NamedTuple without additional allocations.

# Examples
```julia
```

"""
=#
"Convert struct `x` to NamedTuple without additional allocations."
@generated function to_NamedTuple_generated(x)
    nt = Expr(:quote, generate_tuple(x))
    tup = Expr(:tuple)
    for i in 1:fieldcount(x)
        push!(tup.args, :(getfield(x, $i)))
    end
    return :($nt($tup))
end
function generate_tuple(container)
    return NamedTuple{
        (fieldnames(container)...,),
        Tuple{(fieldtype(container, i) for i in 1:fieldcount(container))...},
    }
end

############################################################################################
#=
"""
$(SIGNATURES)
Generate union of fields from 2 NamedTuples x and y.

# Examples
```julia
```

"""
=#
"Generate union of fields from 2 NamedTuples `x` and `y`."
@generated function keyunion(x::NamedTuple{Kx,Tx}, y::NamedTuple{Ky,Ty}) where {Kx,Tx,Ky,Ty}
    return union(Kx, Ky)
end

############################################################################################
"""
$(SIGNATURES)
Return symbol input as tuple.

# Examples
```julia
```

"""
function to_Tuple(sym::Symbol)
    return (sym, )
end
function to_Tuple(sym::Tuple)
    return sym
end

############################################################################################
# Export


text:
export Species, add!, renew!, belongs_to_species, find_species!, compute_fitness_mean!, compute_fitness_max!

mutable struct Species
    id::Int64
    age::Int64
    members::Array{NEATIndiv}
    previous_members::Array{NEATIndiv}
    dist_threshold::Float64
    fitness_val::Float64
end

"Creates a new species from a config dict and stores it"
function Species(cfg::Dict)
    cfg["species_max"] += 1 # New species ID
    s = Species(cfg["species_max"], 0, [], [], cfg["dist_threshold"], 0.0)  # Create empty species
    cfg["Species"][cfg["species_max"]] = s   # Store new species in cfg
    s
end

"Adds indiv to the current members of species s"
function add!(s::Species, indiv::NEATIndiv)
    push!(s.members, indiv)
end

"Empties the list of members if the species has any, and stores them as previous_members"
function renew!(s::Species)
    if length(s.members) > 0
        s.previous_members = s.members
        s.members = []
    end
end

"Returns true if indiv is close enough to s, else false"
function belongs_to_species(s::Species, indiv::NEATIndiv, cfg::Dict)
    if length(s.previous_members) > 0 # Compare first to the alive individuals
        rand_indiv = rand(s.previous_members)
    elseif length(s.members) > 0 # Else compare to last dead ones
        rand_indiv = rand(s.members)
    else # If both are empty, the species is empty
        return true
    end

    distance(indiv, rand_indiv, cfg) <= s.dist_threshold
end

"
Assigns indiv to the first species that matches.
If no species match, creates a new one.
"
function find_species!(indiv::NEATIndiv, cfg::Dict)
    for s in values(cfg["Species"])
        if belongs_to_species(s, indiv, cfg)
            add!(s, indiv)
            return true
        end
    end
    # If no species fits, create a new one
    new_s = Species(cfg)
    add!(new_s, indiv)
    return true
end

"Computes fitness for each member with explicit fitness sharing, returns mean fitness."
function compute_fitness_mean!(s::Species, fitness::Function)
    if length(s.members)==0
        s.fitness_val
        return 0
    end
    s_size = length(s.members) # Species size
    total_fit = 0
    for i in s.members
        f = fitness(i)
        i.fitness .=  f
        total_fit += f[1]/s_size # Compute total adjusted fitness in species
    end

    s.fitness_val = total_fit /s_size # Compute average adjusted fitness in species
    s.fitness_val
end

"Computes fitness for each member with explicit fitness sharing, returns max fitness."
function compute_fitness_max!(s::Species, fitness::Function)
    if length(s.members)==0
        s.fitness_val
        return 0
    end
    s_size = length(s.members) # Species size
    s.fitness_val = 0
    total_fit = 0
    for i in s.members
        f = fitness(i)
        i.fitness .=  f
        # Explicit fitness sharing computed for the species only
        s.fitness_val = maximum([s.fitness_val, f[1]/s_size])  # Compute max adjusted fitness in species
        total_fit += f[1]/s_size
    end
    s.fitness_val
end

"Creates n_children offsprings by crossover or mutation"
function reproduction!(
    s::Species,
    selection::Function,
    n_children::Int64,
    cfg::Dict,
)
    s.age += 1
    renew!(s)
    for i = 1:n_children
        if rand() < cfg["p_mutate_only"] || length(s.previous_members)==1
            # Mutate only
            p1 = selection(s.previous_members)
            child = mutate(p1, cfg)
        else
            # Crossover
            p1 = selection(s.previous_members)
            p2 = selection(s.previous_members)
            if p1 != p2
                child = crossover(p1, p2, cfg)
                child = mutate(child, cfg)
            else
                # If the same individual is chosen twice, mutate it
                child = mutate(p1, cfg)
            end
        end
        add!(s, child)
    end

end


text:
using SeisIO, SeisNoise, JLD2, PlotlyJS

# This script tests the functionality of Stretching.jl from Noise.jl
# This assumes data generated via genSynth.jl

function testStretching(finame::String, foname::String, InputDict::Dict, windowlenlist::Array{Float64,1}, winsteplist::Array{Float64,1}, type::String)
    valData = jldopen(finame)
    dvVlist = valData["info/dvVlist"]
    noiselist = valData["info/noiselist"]

    f_err = jldopen(foname, "a+")
    if type == "realData"
        f_err["info/dvV"] = dvVlist
        f_err["info/noise"] = noiselist
    end

    for dvV in dvVlist
        println("dv/v: $dvV")
        for noiselvl in noiselist
            for win_len in windowlenlist, win_step in windowsteplist
                signals = valData["$type/$dvV.$noiselvl"]
                ref = signals["ref"]
                cur = signals["cur"]

                time_axis = collect(InputDict["mintimelag"]:1/InputDict["fs"]:InputDict["mintimelag"]+InputDict["dtt_width"])

                # allow windowed stretching analysis
                if win_step == 0.0
                    minind = [1]
                    window_length_samples = length(time_axis)
                else
                    window_length_samples = convert(Int,win_len*InputDict["fs"])
                    window_step_samples = convert(Int, win_step*InputDict["fs"])
                    minind = 1:window_step_samples:length(ref) - window_length_samples
                end

                # number of windows
                N = length(minind)
                dv_list = zeros(N)
                err_list = zeros(N)

                # iterate over windows, report average dv/v
                for ii=1:N
                    window = collect(minind[ii]:minind[ii]+window_length_samples-1)
                    dv, cc, cdp, eps, err, C = stretching(ref,
                                                          cur,
                                                          time_axis,
                                                          window,
                                                          InputDict["freqmin"],
                                                          InputDict["freqmax"],
                                                          dvmin=InputDict["dvmin"],
                                                          dvmax=InputDict["dvmax"])
                    dv_list[ii] = dv
                    err_list[ii] = err
                end
                comp_dvV = sum(dv_list)/N
                comp_err = sum(err_list)/N

                f_err["$type/$dvV.$noiselvl"] = [comp_dvV, comp_err]
            end
        end
    end
    close(f_err)
    close(valData)
end

# MWCS input parameters
InputDict_real = Dict( "freqmin"    => 0.1,
                       "freqmax"    => 9.99,
                       "fs"         => 20.0,
                       "mintimelag" => -100.0,
                       "dtt_width"  => 200.0,
                       "dvmin"      => -0.04,
                       "dvmax"      => 0.04,
                       "ntrial"     => 100 )

InputDict_synth = Dict( "freqmin"    => 0.1,
                        "freqmax"    => 9.99,
                        "fs"         => 20.0,
                        "mintimelag" => 0.0,
                        "dtt_width"  => 200.0,
                        "dvmin"      => -0.04,
                        "dvmax"      => 0.04,
                        "ntrial"     => 100 )

finame = "verificationData.jld2"
foname = "verificationDataError_stretching.jld2"
windowlenlist = [200.0]
windowsteplist = [0.0]
type = "rickerConv"

for type in ["rickerConv", "dampedSinusoid", "realData"]
    if type == "realData"
        InputDict = InputDict_real
    else
        InputDict = InputDict_synth
    end
    testStretching(finame, foname, InputDict, windowlenlist, windowsteplist, type)
end

max_stars_repo_name: kura-okubo/SeisXcorrelation

max_issues_repo_name: kura-okubo/SeisXcorrelation

max_forks_repo_name: kura-okubo/SeisXcorrelation


text:
using DelaySSAToolkit, Catalyst
using Test
rn = @reaction_network begin
    1/(1+Y^2), 0 --> X
    1/(1+Y),   Y --> 0
end
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)
# states(rn)
u0 = [0,0]
de_chan0 = [[]]
tf = 400.
tspan = (0,tf)
τ = 20.
dprob = DiscreteProblem(jumpsys, u0, tspan)
# jumpsys.dep_graph

delay_trigger_affect! = function (integrator, rng)
  append!(integrator.de_chan[1], τ)
end
delay_trigger = Dict(1=>delay_trigger_affect!)
delay_complete = Dict(1=>[2=>1, 1=>-1], 2=>[2=>1])
delay_interrupt = Dict()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
alg = DelayDirectCR()
# alg = DelayMNRM()
djprob = DelayJumpProblem(jumpsys, dprob, alg, delayjumpset, de_chan0)

p = djprob.discrete_jump_aggregation;
integrator = DelaySSAToolkit.init(djprob, SSAStepper())
integrator.de_chan
delay_complete[1]
@test p.dep_gr == [[1],[1,2]]
@test DelaySSAToolkit.dep_gr_delay(p, integrator) == Dict(2=>[2],1=>[1,2])


text:
"""
Calculate the euler characteristic of solid, working on V, copEV, copFE matrices
"""
function test_eulero(V, copEV, copFE)
    if !isempty(V)
        vertex=size(V,1)
    else
        println("WARNING!!! V matrix is empty!!!")
    end
    if !isempty(copEV)
        edges=size(copEV,1)
    else
        println("WARNING!!! copEV matrix is empty!!!")
    end
    if !isempty(copFE)
        faces = size(copFE,1)
    else
        println("WARNING!!! copFE matrix is empty!!!")
    end
    bc=size(Lar.Arrangement.biconnected_components(copEV))[1]
    println("EULERIAN CHARACTERISTIC-->$(vertex-edges+faces) BICON COMPS==$(bc) VERTEX = $(vertex) EDGES = $(edges) FACES = $(faces)")
end
"""
check the euler characteristic of solid for a single face at time, working on V, copEV, copFE matrices
"""
function test_eulero_face4face(V, copEV, copFE)
    edges=Set{Int64}()
    vertices=Set{Int64}()
    counter=0
    dict=Dict()
    array=[]
    array2=[]
    ciEV=[]
    ciFE=[]
    if !isempty(V)
        vertex=size(V)[1]
    end
    if !isempty(copEV)
        ciEV=findall(!iszero,copEV);
    end
    if !isempty(copFE)
        ciFE=findall(!iszero,copFE);
    end
    for j in 1:size(copFE)[1]
        for i in ciFE
            if i[1]==j
                push!(array,[i[1],i[2]])
            end
        end
    end
    for j in 1:size(copEV)[1]
        for i in ciEV
            if i[1]==j
                push!(array2,[i[1],i[2]])
            end
        end
    end
    prev=0
    for x in array
        if prev!=x[1]
            if length(vertices)!==length(edges)
                println("WARNING !!! TEST FAILED")
            end
            vertices=Set{Int64}()
            edges=Set{Int64}()
        end
        for y in array2
            if x[2]==y[1]
                push!(edges,y[1])
                push!(vertices,y[2])
            end

        end
        prev=x[1]
    end
end


text:
with_repo_access() do test_repo, path

    context("test walk revlist") do
        walker = GitRevWalker(test_repo)
        push!(walker, Oid("9fd738e8f7967c078dceed8190330fc8648ee56a"))
        data = collect(walker)
        ids = sort(data)
        str = join(map(id -> hex(id)[1:5], ids), ".")
        @test str == "4a202.5b5b0.84960.9fd73"
    end

    context("test walk partial revlist") do
        walker = GitRevWalker(test_repo)
        id = Oid("8496071c1b46c854b31185ea97743be6a8774479")
        push!(walker, id)
        walk = collect(walker)
        @test id == Oid(walk[1])
        @test length(walk) == 1
    end

    context("test hide part of list") do
        walker = GitRevWalker(test_repo)
        push!(walker, Oid("9fd738e8f7967c078dceed8190330fc8648ee56a"))
        hide!(walker, Oid("5b5b025afb0b4c913b4c338a42934a3863bf3644"))
        @test length(collect(walker)) == 2
    end

    context("test reset walker") do
        walker = GitRevWalker(test_repo)
        id = Oid("8496071c1b46c854b31185ea97743be6a8774479")
        push!(walker, id)
        walk = collect(walker)
        @test id == Oid(walk[1])
        @test length(walk) == 1
        reset!(walker)
        walk = collect(walker)
        @test length(walk) == 0
    end

   context("test walker is iterable") do
        walker = GitRevWalker(test_repo)
        push!(walker, Oid("9fd738e8f7967c078dceed8190330fc8648ee56a"))
        vs = sort(collect(walker))
        str = join(map(a -> hex(a)[1:4], vs), ".")
        @test str == "4a20.5b5b.8496.9fd7"
   end
end

function do_sort(walker, sorting; rev=false)
    id = Oid("a4a7dce85cf63874e984719f4fdd239f5145052f")
    sortby!(walker, sorting, rev=rev)
    push!(walker, id)
    return collect(walker)
end

function revlist_with_sorting(walker, sorting; rev=false)
    data = do_sort(walker, sorting, rev=rev)
    return join(map(a -> hex(a)[1:5], data), ".")
end

function is_toposorted(list)
    tfunc = (c) -> begin
        pcs = parents(c)
        for p in pcs
            idx_commit = first(findin(list, [c]))
            idx_parent = first(findin(list, [p]))
            if idx_commit > idx_parent
                return false
            end
        end
        return true
    end
    all(tfunc, list)
end

with_repo_access() do test_repo, path
    context("test walk revlist") do
        walker = GitRevWalker(test_repo)
        t = revlist_with_sorting(walker, :date)
        @test t ==  "a4a7d.c4780.9fd73.4a202.5b5b0.84960"
    end

    context("test sort by topo") do
        walker = GitRevWalker(test_repo)
        sort_list = do_sort(walker, :topo)
        @test is_toposorted(sort_list)
    end

    context("test sort by date reversed") do
        walker = GitRevWalker(test_repo)
        t = revlist_with_sorting(walker, :date, rev=true)
        @test t == "84960.5b5b0.4a202.9fd73.c4780.a4a7d"
    end

    context("test sort by topo reversed") do
        walker = GitRevWalker(test_repo)
        sort_list = do_sort(walker, :topo, rev=true)
        reverse!(sort_list)
        @test is_toposorted(sort_list)
    end
end


text:
module ROCKS

using LinearAlgebra
using FreqTables
using NamedArrays

using Reexport
@reexport using StatsBase
@reexport using DataFrames
@reexport using Plots

export ranks                # equal count binning
export kstest               # two sample KS and location
export concordance          # Concordance with bounds function
export roc                  # ssynonym for Concordance

export BCDiag
export bcdiag
export ksplot, rocplot, biasplot, accuracyplot
export liftcurve, cumliftcurve
export liftable, cumliftable

include("ranks.jl")
include("kstest.jl")
include("roc.jl")
include("bcdiag.jl")

end # module


text:
#! /usr/bin/julia
#
#	csv_read.jl
#
#					Jul/05/2018
# --------------------------------------------------------------------
include("text_manipulate.jl")
#
# --------------------------------------------------------------------
println("*** 開始 ***")
fname_in = ARGS[1]
println(fname_in)
#
dict_aa = csv_to_dict_proc(fname_in)
dict_display_proc(dict_aa)
#
println("*** 終了 ***")
# --------------------------------------------------------------------

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
#=
nufft.jl
Non-uniform FFT (NUFFT), currently a wrapper around NFFT.jl
todo: open issues: small N, odd N, nufft!, adjoint!
2019-06-06, Jeff Fessler, University of Michigan
=#

export nufft_init, nufft_plots, nufft

#using MIRT: dtft_init, map_many
using NFFT
using Plots
using Random: seed!
using LinearAlgebra: norm

"""
`nufft_eltype(w::AbstractArray{<:Real})`
ensure NFFTPlan is Float32 or Float64
"""
function nufft_eltype(w::AbstractArray{<:Number})
	T = eltype(w)
	if T <: Integer || T == Float16
		return Float32 # require at least Float32
	elseif T ∉ (Float32,Float64)
		throw("unknown type $T")
	end
	return T
end

# the following convenience routine ensures correct type passed to nfft()
# see https://github.com/tknopp/NFFT.jl/pull/33
# todo: may be unnecessary with future version of nfft()
"""
`nufft_typer(T::DataType, x::AbstractArray{<:Real}; warn::Bool=true)`
type conversion wrapper for `nfft()`
"""
function nufft_typer(T::DataType, x::AbstractArray{<:Number}; warn::Bool=true)
	if eltype(x) == T
		return x
	end
	@warn("converting $(eltype(x)) to $T")
	return T.(x)
end

"""
`p = nufft_init(w, N; nfft_m=4, nfft_sigma=2.0, pi_error=true, n_shift=0)`

Setup 1D NUFFT,
for computing fast ``O(N \\log N)`` approximation to

``X[m] = sum_{n=0}^{N-1} x[n] exp(-i w[m] (n - n_shift)), m=1,…,M``

in
* `w::AbstractArray{<:Real}` `[M]` frequency locations (units radians/sample)
	+ `eltype(w)` determines the NFFTPlan type; so to save memory use Float32!
* `N::Int` signal length

option
* `nfft_m::Int` 		see NFFT.jl documentation; default 4
* `nfft_sigma::Real`	"", default 2.0
* `n_shift::Real`		often is N/2; default 0
* `pi_error::Bool`		throw error if ``|w| > π``, default `true`
   + Set to `false` only if you are very sure of what you are doing!
* `do_many::Bool`	support extended inputs via `map_many`? default `true`

out
* `p NamedTuple` with fields
	`nufft = x -> nufft(x), adjoint = y -> nufft_adj(y), A=LinearMap`

The default settings are such that for a 1D signal of length N=512,
the worst-case error is below 1e-5 which is probably adequate
for typical medical imaging applications.
To verify this statement, run `nufft_plot1()` and see plot.
"""
function nufft_init(w::AbstractArray{<:Real}, N::Int;
		n_shift::Real = 0,
		nfft_m::Int = 4,
		nfft_sigma::Real = 2.0,
		pi_error::Bool = true,
		do_many::Bool = true,
	)

	N < 6 && throw("NFFT may be erroneous for small N")
	isodd(N) && throw("NFFT erroneous for odd N")
	pi_error && any(abs.(w) .> π) &&
		throw(ArgumentError("|w| > π is likely an error"))

	T = nufft_eltype(w)
	CT = Complex{T}
	f = convert.(T, w/(2π)) # note: NFFTPlan must have correct type
	p = NFFTPlan(f, N, nfft_m, nfft_sigma) # create plan
	M = length(w)
	# extra phase here because NFFT always starts from -N/2
	phasor = CT.(cis.(-w * (N/2 - n_shift)))
	phasor_conj = conj.(phasor)
	forw1 = x -> nfft(p, nufft_typer(CT, x)) .* phasor
#	forw! = x,y -> nfft!(p, nufft_typer(CT, x)) .* phasor # todo
	back1 = y -> nfft_adjoint(p, nufft_typer(CT, y .* phasor_conj))

	A = LinearMap{CT}(x -> forw1(x), y -> back1(y), M, N) # no "many" here!

	if do_many
		forw = x -> map_many(forw1, x, (N,))
		back = y -> map_many(back1, y, (M,))
	else
		forw = forw1
		back = back1
	end

	return (nufft=forw, adjoint=back, A=A)
end

"""
`p = nufft_init(w, N; nfft_m=4, nfft_sigma=2.0, pi_error=true, n_shift=?)`

Setup multi-dimensional NUFFT,
for computing fast ``O(N \\log N)`` approximation to

``X[m] = sum_{n=0}^{N-1} x[n] exp(-i w[m,:] (n - n_shift)), m=1,…,M``

in
* `w::AbstractMatrix{<:Real}` `[M,D]` frequency locations (units radians/sample)
	+ `eltype(w)` determines the NFFTPlan type; so to save memory use Float32!
* `N::Dims` `[D]` signal dimensions

option
* `nfft_m::Int` 		see NFFT.jl documentation; default 4
* `nfft_sigma::Real`	"", default 2.0
* `n_shift::AbstractVector{<:Real}`	`[D]`	often is N/2; default zeros(D)
* `pi_error::Bool`		throw error if ``|w| > π``, default `true`
   + Set to `false` only if you are very sure of what you are doing!
* `do_many::Bool`	support extended inputs via `map_many`? default `true`

The default `do_many` option is designed for parallel MRI where the k-space
sampling pattern applies to every coil.
It may also be useful for dynamic MRI with repeated sampling patterns.
The coil and/or time dimensions must come after the spatial dimensions.

out
* `p NamedTuple` with fields
	`nufft = x -> nufft(x), adjoint = y -> nufft_adj(y), A=LinearMap`
	(The LinearMap does not support `do_many`.)
"""
function nufft_init(w::AbstractMatrix{<:Real},
		N::Dims;
		n_shift::AbstractVector{<:Real} = zeros(Int, length(N)),
		nfft_m::Int = 4,
		nfft_sigma::Real = 2.0,
		pi_error::Bool = true,
		do_many::Bool = true,
	)

	any(N .< 6) && throw("NFFT may be erroneous for small N")
	any(isodd.(N)) && throw("NFFT erroneous for odd N")
	pi_error && any(abs.(w) .> π) &&
		throw(ArgumentError("|w| > π is likely an error"))

	M,D = size(w)
	length(N) != D && throw(DimensionMismatch("length(N) vs D=$D"))
	length(n_shift) != D && throw(DimensionMismatch("length(n_shift) vs D=$D"))

	T = nufft_eltype(w)
	CT = Complex{T}
	f = convert.(T, w/(2π)) # note: NFFTPlan must have correct type
	p = NFFTPlan(f', N, nfft_m, nfft_sigma) # create plan

	# extra phase here because NFFT always starts from -N/2
	phasor = CT.(cis.(-w * (collect(N)/2. - n_shift)))
	phasor_conj = conj.(phasor)
	forw1 = x -> nfft(p, nufft_typer(CT, x)) .* phasor
	back1 = y -> nfft_adjoint(p, nufft_typer(CT, y .* phasor_conj))

	# no "many" for LinearMap:
	A = LinearMap{CT}(x -> forw1(reshape(x,N)), y -> back1(y)[:], M, prod(N))

	if do_many
		forw = x -> map_many(forw1, x, N)
		back = y -> map_many(back1, y, (M,))
	else
		forw = forw1
		back = back1
	end

	return (nufft=forw, adjoint=back, A=A)
end

"""
`nufft_test1(; M=30, N=20, n_shift=1.7, T=?, tol=?)`
simple 1D tests
"""
function nufft_test1( ;
		M::Int = 30, N::Int = 20, n_shift::Real = 1.7,
		T::DataType = Float64, tol::Real = 1e-6)
	seed!(0)
	w = (rand(M) .- 0.5) * 2 * pi
	w = T.(w)
	x = complex.(randn(T, N), randn(T, N))
	sd = dtft_init(w, N; n_shift=n_shift)
	sn = nufft_init(w, N, n_shift=n_shift)
	o0 = sd.dtft(x)
	y = Complex{T}.(o0 / norm(o0))
	a0 = sd.adjoint(y)
	o1 = sn.nufft(x)
	o2 = sn.A * x
	a1 = sn.adjoint(y)
	a2 = sn.A' * y
	@test norm(o1 - o0, Inf) / norm(o0, Inf) < tol
	@test isequal(o1, o2)
	@test norm(a1 - a0, Inf) / norm(a0, Inf) < tol
	@test isequal(a1, a2)
	@test isapprox(Matrix(sn.A)', Matrix(sn.A')) # 1D adjoint test

	sn = nufft_init(w, N, n_shift=n_shift, do_many=false)
	o3 = sn.nufft(x)
	@test norm(o3 - o0, Inf) / norm(o0, Inf) < tol

	sn.nufft(ones(Int,N)) # should produce a warning
	true
end

"""
`nufft_test2(; M=?, N=?, n_shift=?, T=?, tol=?)`
simple 2D test
"""
function nufft_test2(;
		M::Int = 31,
		N::Dims = (10,8),
		n_shift::AbstractVector{<:Real} = [4,3],
		T::DataType = Float64, tol::Real = 2e-6)

	seed!(0)
	w = []

#=
	# fft sampling
	M = prod(N)
	w1 = (2*pi) * (0:N[1]-1) / N[1]
    w2 = (2*pi) * (0:N[2]-1) / N[2]
    w1 = repeat(w1, 1, N[2])
    w2 = repeat(w2', N[1], 1)
    w = [w1[:] w2[:]]
	n_shift = [0,0]
=#

	w = (rand(M,2) .- 0.5) * 2 * pi

	w = T.(w)
	sd = dtft_init(w, N; n_shift=n_shift)
	sn = nufft_init(w, N, n_shift=n_shift, pi_error=false)

	x = complex.(randn(T, N), randn(T, N))
	o0 = sd.dtft(x)
	o1 = sn.nufft(x)
	@test norm(o1 - o0, Inf) / norm(o0, Inf) < tol

#=
	# fft test only
	o2 = fft(x)
	o0 = reshape(o0, N)
	o1 = reshape(o1, N)
	@show norm(o2 - o0, Inf) / norm(o0, Inf)
	@show norm(o1 - o0, Inf) / norm(o0, Inf)
	@show norm(o2 - o1, Inf) / norm(o0, Inf)
=#

	y = Complex{T}.(o0 / norm(o0))
	a0 = sd.adjoint(y)
	a1 = sn.adjoint(y)
	@test norm(a1 - a0, Inf) / norm(a0, Inf) < tol
	o2 = sn.A * x[:]
	@test isequal(o1, o2)
	a2 = sn.A' * y
	a2 = reshape(a2, N)
	@test isequal(a1, a2)
	@test isapprox(Matrix(sn.A)', Matrix(sn.A')) # 2D adjoint test

	@test isequal(sn.nufft(cat(dims=4, x, 2x)),
		cat(dims=3, sn.nufft(x), sn.nufft(2x)))
	@test isequal(sn.adjoint(cat(dims=3, y, 2y)),
			cat(dims=4, sn.adjoint(y), sn.adjoint(2y)))

	sn = nufft_init(w, N, n_shift=n_shift, pi_error=false, do_many=false)
	o3 = sn.nufft(x)
	@test norm(o3 - o0, Inf) / norm(o0, Inf) < tol
	true
end

"""
`w, errs = nufft_errors(; M=?, w=?, N=?, n_shift=?, ...)`

Compute worst-case errors for NUFFT (for signal of length N of unit norm)
"""
function nufft_errors(;
	M::Int = 401,
	N::Int = 512,
	w::AbstractArray{<:Real} = LinRange(0, 2π/N, M),
	n_shift::Real = 0,
	kwargs...)

	sd = dtft_init(w, N; n_shift=n_shift)
	sn = nufft_init(w, N; n_shift=n_shift, kwargs...)
	E = Matrix(sn.A - sd.A)
	return w, mapslices(norm, E, dims=2)[:] # [M]
end

"""
`nufft_plot1()`

Plot worst-case error over all frequencies w between 0 and 2pi/N for various N.
"""
function nufft_plot1()
	Nlist = 2 .^ (4:9)
	elist = zeros(length(Nlist))
	for ii=1:length(Nlist)
		w, errs = nufft_errors(; N = Nlist[ii])
	#	plot(w*N/2/π, tmp)
		elist[ii] = maximum(errs)
	end
	scatter(Nlist, elist, xtick=Nlist, label="", xlabel="N", ylabel="error")
end

"""
`nufft_plot_error_m(; mlist=?)`

plot error vs NFFT sigma
"""
function nufft_plot_error_m(;
	mlist::AbstractArray{<:Int} = 3:7)
	worst = zeros(length(mlist))
	for jm = 1:length(mlist)
		_, errs = nufft_errors(; nfft_m=mlist[jm])
		worst[jm] = maximum(errs)
	end
	scatter(mlist, worst, xlabel="m", ylabel="error", label="")
end

"""
`nufft_plot_error_s(; slist=?)`

plot error vs NFFT sigma
"""
function nufft_plot_error_s(;
	slist::AbstractArray{<:Real} = [1.5; 2:6])
	worst = zeros(length(slist))
	for is = 1:length(slist)
		_, errs = nufft_errors(; nfft_sigma=slist[is])
		worst[is] = maximum(errs)
	end
	scatter(slist, worst, xlabel="σ", ylabel="error", label="")
end

"""
`nufft_plots()`
various NUFFT error plots
"""
function nufft_plots()
	p1 = nufft_plot_error_s()
	p2 = nufft_plot_error_m()
	p3 = nufft_plot1()
	plot(p1, p2, p3)
end

"""
`nufft(:test)`
self tests
"""
function nufft(test::Symbol)
	test != :test && throw("bad symbol $test")
	@test_throws String nufft_eltype(ones(BigFloat,3))
	@test_throws String nufft_init([0], 2)
	@test_throws String nufft_init([0], 7)
	@test_throws ArgumentError nufft_init([2π], 8)
	@test nufft_plots() isa Plots.Plot
	@test nufft_test1()
	@test nufft_test1(; T=Float32)
	@test nufft_test2()
	@test nufft_test2(; T=Float32)
	true
end

#=
	# todo: 1d vs 2d
	M = 4
	N = (M,1)
	seed!(0)
	w = (0:(M-1))/M * 2 * pi
	w = [w zeros(M)]
	sd = dtft_init(w, N)
	Ad = Matrix(sd.A)
	sn1 = nufft_init(w[:,1], N[1]; pi_error=false)
	An1 = Matrix(sn1.A)
	@show maximum(abs.(An1 - Ad))
	sn2 = nufft_init(w, N; pi_error=false)
	An2 = Matrix(sn2.A)
#	@show o0 = sd.dtft(x)
#	@show o1 = sn.nufft(x)
#	@show o1-o0
#	@test maximum(abs.(o1 - o0)) < tol
=#

#=
	# todo MWE for 1D vs 2D
	M = 6
	x1 = collect((-Int(M/2)):(Int(M/2)-1))/M
	N1 = M
	p1 = NFFTPlan(x1, N1)
#	p1 = NFFTPlan(x1', (N1,)) # this works too
#	f1 = ones(N1)
	f1 = [1.,2,3,4,0,0]
	o1 = nfft(p1, complex(f1)) # slightly annoying to have to use complex() here

	N2 = (M,4) # works fine with "4" instead of "1" here and separable signal
	x2 = [x1 zeros(M)]
	p2 = NFFTPlan(x2', N2)
	f2 = f1 * [0,0,1,0]'
	o2 = nfft(p2, complex(f2)) # ditto

	display(round.(o1, digits=7))
	display(round.(o2, digits=7))

#	tmp = nufft_init(x1*2*pi, N1)
=#

#nufft(:test)
#nufft_test1(;N=41) # todo: inaccurate


text:
export @p_str, StaticPolynomial, differentiate, integrate, degree, coeffs, variable, isconstant

abstract type AbstractStaticPolynomial{T} <: Function end

struct StaticPolynomial{T<:Number, N} <: AbstractStaticPolynomial{T}
    coeffs::NTuple{N, T}
    var::Symbol
    function StaticPolynomial(coeffs::Tuple, var::Symbol = :x)
        T = supertype_all(coeffs...)
        isempty(coeffs) && return StaticPolynomial{T, 0}()
        N = length(coeffs)
        c = promote(coeffs...) |> NTuple{N, T}
        var = N > 1 ? var : :CONSTANT
        var = isconstant(coeffs) ? :CONSTANT : var
        return new{T, N}(c, var)
    end
    StaticPolynomial{T, N}(coeffs::NTuple{N, T}, var::Symbol) where {T<:Number, N}= new{T, N}(coeffs, var)
end

coeffs(p::AbstractStaticPolynomial) = p.coeffs
variable(p::AbstractStaticPolynomial) = p.var

(p::StaticPolynomial)(x) = evalpoly(x, coeffs(p))

StaticPolynomial(coeffs::AbstractVector, var::Symbol = :x) = StaticPolynomial(Tuple(coeffs), var)

# Constructors
function StaticPolynomial(spcoeffs::AbstractVector{T}, exponents::AbstractVector{<:Integer}, var::Symbol = :x) where T<:Number
    coeffs = zeros(T, maximum(exponents) + 1)
    for (c, e) in zip(spcoeffs, exponents)
        coeffs[e + 1] += c
    end
    return StaticPolynomial(coeffs, var)
end

function StaticPolynomial(spcoeffs::Tuple, exponents::Tuple, var::Symbol = :x)
    coeffs = zeros(supertype_all(spcoeffs...), maximum(exponents) + 1)
    for (c, e) in zip(spcoeffs, exponents)
        coeffs[e + 1] += c
    end
    return StaticPolynomial(coeffs, var)
end

StaticPolynomial(coeffs...) = StaticPolynomial(coeffs)

# Indexing and utility functions
degree(p::StaticPolynomial{T, N}) where {T<:Number, N} = N - 1
Base.length(p::StaticPolynomial{T, N}) where {T<:Number, N} = N

Base.getindex(p::StaticPolynomial{T, N}, inds...) where {T<:Number, N} = try
        coeffs(p)[inds...]
    catch BoundsError
        zero(T)
    end

# Conversion and promotion rules
Base.convert(::Type{T}, x::StaticPolynomial) where T<:StaticPolynomial = T(x)
Base.convert(::Type{T}, x::StaticPolynomial) where T<:Number = T(x)

StaticPolynomial(p::StaticPolynomial) = p

StaticPolynomial{T}(p::StaticPolynomial{S}) where {T<:Number, S<:Number} = StaticPolynomial(broadcast(T, p.coeffs))

function StaticPolynomial{T, N}(p::StaticPolynomial{S, M}) where {T<:Number, S<:Number, M, N}

    N < M && error("Cannot convert a polynomial of degree $(M-1) to degree $(N-1)")
    coeffs = zeros(T, N)
    for (i, c) in enumerate(p.coeffs)
        coeffs[i] = c
    end
    return StaticPolynomial(coeffs, p.var)
end

function StaticPolynomial{T, N}(p::StaticPolynomial{T, M}) where {T<:Number, N, M}
    N < M && error("Cannot convert a polynomial of degree $(M-1) to degree $(N-1)")
    p1 = StaticPolynomial(zeros(T, N))::StaticPolynomial{T, N}
    return p + p1
end

(::Type{T})(p::StaticPolynomial) where {T<:Number} =
    degree(p) |> iszero ? T(coeffs(p)[1])::T : throw(InexactError(nameof(T), T, p))

Base.promote_rule(::Type{StaticPolynomial{T, N}}, ::Type{S}) where {T<:Number, S<:Number, N} =
    StaticPolynomial{promote_type(T, S), N}

Base.promote_rule(::Type{StaticPolynomial{T, N}}, ::Type{StaticPolynomial{S, N}}) where {T<:Number, S<:Number, N} =
    StaticPolynomial{promote_type(T, S), N}

Base.promote_rule(::Type{StaticPolynomial{T, N}}, ::Type{StaticPolynomial{S, M}}) where {T<:Number, S<:Number, M, N} =
    StaticPolynomial{promote_type(T, S), max(N, M)}

# define basic arithmetic operations
Base.isapprox(p1::StaticPolynomial{T, N}, p2::StaticPolynomial{S, M}) where {T<:Number, S<:Number, M, N} =
    p1.var == p2.var && M == N && all([p1[i] == p2[i] for i in 1:M])

Base.:(==)(p1::StaticPolynomial, p2::StaticPolynomial) = (p1.coeffs == p2.coeffs && p1.var == p2.var) || (iszero(p1) && iszero(p2)) || (isone(p1) && isone(p2))

isconstant(p::StaticPolynomial) = isconstant(p.coeffs) && p.var == :CONSTANT
isconstant(coeffs) = all(iszero.(coeffs[2:end]))

function polynomial_arithmetic(operation::Function, p1::StaticPolynomial{T, N}, p2::StaticPolynomial{S, M}) where {T<:Number, S<:Number, M, N}
    R = promote_type(T, S)
    #@show p1
    #@show p2
    if p1.var == p2.var || isconstant(p1) || isconstant(p2)
        coeffs = (R(operation(p1[i], p2[i])) for i in 1:max(M, N)) |> collect
        if isconstant(p1)
            return StaticPolynomial(coeffs, p2.var)
        else
            return StaticPolynomial(coeffs, p1.var)
        end
    else
        error("Polynomials must be in terms of the same variable.")
    end
end

Base.:(+)(p1::StaticPolynomial, p2::StaticPolynomial) = polynomial_arithmetic(+, p1, p2)
Base.:(+)(p::StaticPolynomial, n::S) where {S<:Number} = p + StaticPolynomial(n)
Base.:(+)(n::S, p::StaticPolynomial) where {S<:Number} = p + n

Base.:(-)(p1::StaticPolynomial, p2::StaticPolynomial) = polynomial_arithmetic(-, p1, p2)
Base.:(-)(p::StaticPolynomial, n::S) where {S<:Number} = p - StaticPolynomial(n)
Base.:(-)(n::S, p::StaticPolynomial) where {S<:Number} = StaticPolynomial(n) - p

Base.:(/)(p::StaticPolynomial, n::S) where {S<:Number} = StaticPolynomial(coeffs(p)./n)

Base.:(//)(p::StaticPolynomial, n::S) where {S<:Number} = StaticPolynomial(coeffs(p).//n)

function Base.:(^)(p::StaticPolynomial, n::S) where {S<:Integer}
    if n == 1
        return p
    elseif n == 0
        return one(p)
    else
        return p * p^(n-1)
    end
end

Base.zero(p::StaticPolynomial{T, N}) where {T<:Number, N} = zero(StaticPolynomial{T, N}, p.var)
Base.zero(::Type{StaticPolynomial{T, N}}, var=:x) where {T, N} = StaticPolynomial(zeros(T, N), :CONSTANT)
Base.zero(::Type{StaticPolynomial{T}}) where T = zero(StaticPolynomial{T, 1})
Base.zero(::Type{StaticPolynomial}) = zero(StaticPolynomial{Int})

Base.one(p::StaticPolynomial{T, N}) where {T<:Number, N} = one(StaticPolynomial{T, N}, p.var)
Base.one(::Type{StaticPolynomial{T, N}}, var=:x) where {T, N} = zero(StaticPolynomial{T, N}, var) + 1
Base.one(::Type{StaticPolynomial{T}}) where T = one(StaticPolynomial{T, 1})
Base.one(::Type{StaticPolynomial}) = one(StaticPolynomial{Int})

function Base.:(*)(p1::StaticPolynomial{T, N}, p2::StaticPolynomial{S, M}) where {T<:Number, S<:Number, N, M}
    !(p1.var == p2.var || p1.var == :CONSTANT || p2.var == :CONSTANT) && error("Polynomials must be in terms of the same variable.")

    if M > 1 && N > 1
        degree = M + N - 2
        R = promote_type(T, S)
        cs = zeros(R, degree + 1)
        for i in 1:N
            for j in 1:M
                cs[i + j - 1] += p1[i] * p2[j]
            end
        end
        return StaticPolynomial(cs, p1.var)
    elseif N > 1
        return p1 * coeffs(p2)[1]
    else
        return p2 * coeffs(p1)[1]
    end
end

function Base.:(*)(p::StaticPolynomial{T}, n::S) where {T<:Number, S<:Number}
    if iszero(n)
        return StaticPolynomial(n)
    elseif isone(n)
        R = promote_type(T, S)
        return StaticPolynomial{R}(p)
    else
        return StaticPolynomial(n .* p.coeffs, p.var)
    end
end

Base.:(*)(n::S, p::StaticPolynomial) where {S<:Number} = p * n

Base.iszero(p::StaticPolynomial) = map(iszero, p.coeffs) |> all
Base.isone(p::StaticPolynomial) = isone(p.coeffs[1]) && map(iszero, p.coeffs[2:end]) |> all

# Differentiation and integration
function differentiate(p::StaticPolynomial, order = 1)
    if order == 0 || (iszero(p) && length(p) <= 1)
        return p
    elseif length(p) == 1 && !iszero(p)
        return zero(p)
    else
        coeffs = (i * p[i + 1] for i in 1:length(p)-1) |> collect
        var = length(coeffs) <= 1 ? :CONSTANT : p.var
        p = StaticPolynomial(coeffs, var)
        return differentiate(p, order - 1)
    end
end

function integrate(p::StaticPolynomial{T}; C = 0) where T<:AbstractFloat
    n = length(p)
    cs = zeros(T, n + 1)
    cs[1] = C
    cs[2:end] = [p[i] / i for i in 1:n]
    var = p.var == :CONSTANT ? :x : p.var
    return StaticPolynomial(cs, var)
end

function integrate(p::StaticPolynomial; C = 0)
    n = length(p)
    cs = zeros(Rational, n + 1)
    cs[1] = C
    cs[2:end] = [p[i] // i for i in 1:n]
    var = p.var == :CONSTANT ? :x : p.var
    return StaticPolynomial(cs, var)
end

#StaticPolynomial IO stuff
function exponent_string(exponent::Integer, mimetype=MIME"text/plain"())
    if exponent < 0
        return "⁻" * exponent_string(-exponent, mimetype)
    elseif exponent > 9
        return exponent_string(exponent ÷ 10, mimetype) * exponent_string(exponent % 10, mimetype)
    else
        if mimetype === MIME"text/latex"()
            return "^{" * string(exponent) * "}"
        else
            return exponent_dict[exponent]
        end
    end
end

const exponent_dict = Dict{Int, String}(
    0 => "⁰",
    1 => "¹",
    2 => "²",
    3 => "³",
    4 => "⁴",
    5 => "⁵",
    6 => "⁶",
    7 => "⁷",
    8 => "⁸",
    9 => "⁹",
)

poly_coeff_string(coeff, exponent) = isone(coeff) ? "" : (exponent == 0 ? string(coeff) : _poly_coeff_string(coeff))

_poly_coeff_string(coeff::(BigFloat | Float16 | Float32 | Float64 | Integer)) = string(coeff)
_poly_coeff_string(coeff) = string(coeff) |> wrap_string

wrap_string(str; left='(', right=')') = left * str * right

function poly_term_string(coeff, exponent, var, mimetype)
    if iszero(exponent)
        return isone(coeff) ? string(coeff) : poly_coeff_string(coeff, exponent)
    elseif isone(exponent)
        return poly_coeff_string(coeff, exponent) * string(var)
    else
        return poly_coeff_string(coeff, exponent) * string(var) * exponent_string(exponent, mimetype)
    end
end

signstring(n::Number) = signbit(n) ? "-" : "+"

function Base.show(io::IO, p::StaticPolynomial{T, N}) where {T<:Number, N}
    print(io, polyname(p), "(")
    print(io, printpoly(p))
    print(io, ")")
end

polyname(p) = string(typeof(p))

function printpoly(p::StaticPolynomial{T, N}, mimetype=MIME"text/plain"()) where {T<:Number, N}
    str = ""
    if iszero(p)
        str = string(zero(T))
    else
        if T <: Real
            for (i, c) in enumerate(coeffs(p))
                if c != 0
                    str *= isempty(str) ? poly_term_string(c, i - 1, p.var, mimetype) : " " * signstring(c) * " " * poly_term_string(abs(c), i-1, p.var, mimetype)
                end
            end
        else
            for (i, c) in enumerate(coeffs(p))
                if c != 0
                    str *= isempty(str) ? "" : " + "
                    str *= poly_term_string(c, i - 1, p.var, mimetype)
                end
            end
        end
    end
    return str
end

Base.string(p::StaticPolynomial) = printpoly(p)

Base.show(io::IO, mimetype::MIME"text/plain",  p::StaticPolynomial{T, N}) where {T<:Number, N} = begin
    Base.show(io, p)
end

function Base.show(io::IO, ::MIME"text/latex",  p::StaticPolynomial{T, N}) where {T, N}
    print(io, "\$")
    print(io, printpoly(p, MIME"text/latex"()))
    print(io, "\$")
end


text:
module CUSPARSE

using ..CuArrays
using ..CuArrays: unsafe_free!, @argout, @workspace

using CUDAapi

using CUDAdrv
using CUDAdrv: CUstream

using CUDAnative

using CEnum

const SparseChar = Char

const libcusparse = Ref("libcusparse")

# core library
include("libcusparse_common.jl")
include("error.jl")
include("libcusparse.jl")

# low-level wrappers
include("array.jl")
include("util.jl")
include("wrappers.jl")

# high-level integrations
include("interfaces.jl")

const created_handles = IdDict{CuContext,cusparseHandle_t}()
const active_handles = Vector{Union{Nothing,cusparseHandle_t}}()

function handle()
    tid = Threads.threadid()
    if @inbounds active_handles[tid] === nothing
        ctx = context()
        active_handles[tid] = get!(created_handles, ctx) do
            handle = cusparseCreate()
            atexit(()->CUDAdrv.isvalid(ctx) && cusparseDestroy(handle))
            handle
        end
    end
    @inbounds active_handles[tid]
end

function __init__()
    resize!(active_handles, Threads.nthreads())
    fill!(active_handles, nothing)

    CUDAnative.atcontextswitch() do tid, ctx
        # we don't eagerly initialize handles, but do so lazily when requested
        active_handles[tid] = nothing
    end
end

end


text:
module TestBlockFactorizations
using LinearAlgebra
using BlockFactorizations
using Test

@testset "block" begin

    element_types = (Float32, Float64, ComplexF32, ComplexF64)
    for elty_A in element_types
        @testset "eltype $elty_A" begin
            # strided block matrix
            di = 2
            dj = 3
            n = 5
            m = 7
            stride_matrices = [
                                [randn(elty_A, di, dj) for i in 1:n, j in 1:m], # generic strided matrix
                                Diagonal([randn(elty_A, di, dj) for i in 1:n]), # diagonal strided matrix
                                # IDEA: test with scalar elements as fallback
                              ]
            for A in stride_matrices
                # test general constructor
                di, dj = size(A[1])
                n, m = size(A)
                strided_factorizations = [BlockFactorization(A),
                                          BlockFactorization(A, isstrided = true),
                                          BlockFactorization(A, di, dj)]
                # stride_matrices = [A, A, A, D, D, D] # matrices corresponding to strided_factorizations
                for F in strided_factorizations
                    M = Matrix(F)
                    @test size(F) == (size(A, 1) * size(A[1], 1), size(A, 2) * size(A[1], 2)) # (n*di, m*dj)
                    @test F.nindices isa AbstractRange # testing that stridedness was correctly identified
                    @test F.mindices isa AbstractRange
                    @test eltype(F) == elty_A
                    @test eltype(M) == elty_A
                    for ni in 2:2, mi in 1:1
                        for i in 1:di, j in 1:dj
                            @test F[di*(ni-1)+i, dj*(mi-1)+j] == A[ni, mi][i, j]
                            @test M[di*(ni-1)+i, dj*(mi-1)+j] == A[ni, mi][i, j]
                        end
                    end
                    # matrix multiply
                    for elty_B in element_types
                        x = randn(elty_B, size(M, 2))
                        @test M*x ≈ F*x
                        k = 3
                        X = randn(elty_B, size(M, 2), k)
                        @test M*X ≈ F*X
                        # IDEA: test left multiply
                    end
                    # diagonal of factorization
                    diagF = diag(F)
                    @test length(diagF) == minimum(size(M))
                    @test diagF ≈ diag(M)
                end
                allo_1 = @allocated BlockFactorization(A, isstrided = true)
                allo_2 = @allocated BlockFactorization(A, isstrided = false)
                @test allo_1 < allo_2
            end

            # general, non-strided block matrix
            n, m = 4, 5 # this is the size of the non-block equivalent of the following matrices
            nindices = [1, 2, n+1]
            mindices = [1, 4, m+1]
            A = fill(randn(elty_A, 0, 0), 2, 2)
            A[1, 1] = randn(elty_A, 1, 3)
            A[1, 2] = randn(elty_A, 1, 2)
            A[2, 1] = randn(elty_A, 3, 3)
            A[2, 2] = randn(elty_A, 3, 2)
            # D = Diagonal([randn(elty_A, 1, 3), randn(elty_A, 3, 2)])
            matrices = [A, Diagonal(A)]

            for A in matrices
                factorizations = [BlockFactorization(A),
                                  BlockFactorization(A, isstrided = false),
                                  BlockFactorization(A, nindices, mindices)
                                  ]
                for F in factorizations
                    M = Matrix(F)
                    @test size(F) == size(M)
                    for i in 1:n, j in 1:m
                        @test F[i, j] == M[i, j]
                    end

                    # matrix multiply
                    for elty_B in element_types
                        x = randn(elty_B, size(M, 2))
                        @test M*x ≈ F*x
                        k = 3
                        X = randn(elty_B, size(M, 2), k)
                        @test M*X ≈ F*X
                        # IDEA: test left multiply
                    end
                end
            end

            # testing with vector elements
            n, m = 2, 5
            d = 3
            A = [randn(elty_A, d) for i in 1:n, j in 1:m]
            F = BlockFactorization(A)
            M = Matrix(F)
            for elty_B in element_types
                x = randn(elty_B, size(M, 2))
                @test M*x ≈ F*x
                k = 3
                X = randn(elty_B, size(M, 2), k)
                @test M*X ≈ F*X
                # IDEA: test left multiply
            end
            diagF = diag(F)
            @test length(diagF) == minimum(size(M))
            @test diagF ≈ diag(M)

            # testing with heterogeneous block types
            n = 3
            A11 = Diagonal(randn(elty_A, n))
            A12 = randn(elty_A, n, n)
            A21 = UpperTriangular(randn(elty_A, n, n))
            A22 = cholesky(Hermitian(A21'A21))
            A = reshape([A11, A12, A21, A22], 2, :)
            @test eltype(A) == Any
            B = BlockFactorization(A)
            @test eltype(B) == elty_A

            d, n = 2, 3
            A = Diagonal([randn(elty_A, d, d) for _ in 1:n])
            B = BlockFactorization(A)
            @test B isa BlockDiagonalFactorization
        end
    end
end

end # TestBlockFactorizations


text:
include("lib_for_practic_10.jl")

#r = Robot()
#r = BorderRobot{CoordRobot}(r)
#r = PutmarkersRobot{BorderRobot{Robot}}(r)
#show!(r.robot.robot)

function all_kvadro(r)

    x = movements!(r,West)
    y = movements!(r,Sud)
    snake!(r)

    movements!(r,West)
    movements!(r,Sud)
    movements!(r,Nord, y)
    movements!(r,Ost, x)

end

all_kvadro(r)
show(r.robot.robot)

max_stars_repo_name: AntonDudDm/mirea-progs

max_issues_repo_name: AntonDudDm/mirea-progs

max_forks_repo_name: AntonDudDm/mirea-progs


text:
module Sputnik

using Statistics
using Interact
using StatsPlots
using DataStructures
using JuliaDB
import Observables: AbstractObservable, observe
using GroupedErrors
using JuliaDB, IndexedTables
using Parameters
using Blink, Mux

import IndexedTables: AbstractIndexedTable
import TableWidgets, TableView

const tablefolder = joinpath(homedir(), ".sputnik", "tables")
const plotfolder = joinpath(homedir(), ".sputnik", "plots")

include("selectdata.jl")
include("defaultstyles.jl")
include("process.jl")
include(joinpath("techniques", "statplotsrecipe.jl"))
include(joinpath("techniques", "groupederrors.jl"))
include(joinpath("techniques", "table.jl"))

include("dropdownmenus.jl")
include("stylechooser.jl")
include("loadbutton.jl")
include("localsave.jl")
include("launch.jl")
include("plot.jl")

end # module


text:
function square_wave()

    g = kron([1, -1, 1, -1, 1], ones(5))

    M_bf = 5

    ζ_vec = [100, 30, 10, 3, 1, 0.3, 0.1, 0.03, 0.01]

    # For each m there are many non-unique solution
    I_sols =
    [Vector{Int}() for m=1:24]

    f_sols =
        [[24.0,
        18.10194908759502,
        15.757575757575758,
        6.602807597027249,
        4.84848484848485,
        4.285714285714288,
        2.4242424242424248,
        2.1356521739130434];
        zeros(Float64, 16)]

    return g, ζ_vec, I_sols, f_sols
end


text:
@testset "bigleaf_constants" begin
    cst = BigleafConstants()
    cst_moon = BigleafConstants(g = 9.81/6)
    f1(;cst = BigleafConstants()) = 2*cst.g
    @test f1() == 2*cst.g
    @test f1(;cst=cst_moon) == 2*cst_moon.g
    #@code_llvm f1()
    #@code_llvm f1(;cst=BigleafConstants(g = 9.81/6))
    #const cst_moon2 = BigleafConstants(g = 9.81/6)
    #@code_llvm f1(;cst=cst_moon2)
    #@code_llvm f1(;cst=BigleafConstants()) # not optimized
    f2(Tair;cst = BigleafConstants()) = Tair*oftype(Tair,cst.g)
    @test f2(2.0) == 2*cst.g
    #@code_llvm f2(Float32(2.0))
    Tair32 = Float32(2.0)
    @test f2(Tair32; cst=cst_moon) == Float32(2*cst_moon.g)
    #@code_llvm f2(Tair32;cst=cst_moon)
end


text:
function newbound_rounding_details(U,V,Xn)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)
  nU = size(U,1)
  nV = size(V,1)
  r = size(U,2)
  assert(r==size(V,2))

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

  P = spzeros(nU,nV)
  allrecoveries = zeros(size(U,2))
  allrecoveries_accum = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei1 = U_sortperm[1:lastidpos,i]
    ej1 = V_sortperm[1:lastidpos,i]
    ei2 = U_sortperm[nU-lneg:nU,i]
    ej2 = V_sortperm[nV-lneg:nV,i]
    ei = vcat(ei1,ei2)
    ej = vcat(ej1,ej2)
    allrecoveries[i] = evaluate_erdosreyni_experiment(A,B,ei,ej)
    P = P + sparse(ei,ej,1,n,n)

    Xsample = spones(P).*Xn
    maunion,mbunion = edge_list_local(bipartite_matching_local(Xsample))
    allrecoveries_accum[i] = evaluate_erdosreyni_experiment(A,B,maunion,mbunion)

  end
  return P,allrecoveries_accum,allrecoveries
end

function newbound_rounding_max_approx_eff(U,V)
  nU = size(U,1)
  nV = size(V,1)
  nU = size(U,1)
  nV = size(V,1)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)

  r = size(U,2)
  assert(r==size(V,2))

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

  weights = zeros(size(U,2))
  di = zeros(Float64,size(U,2))
  Umatching = Array{Int64,1}[]
  Vmatching = Array{Int64,1}[]

  for i = 1:size(U,2)
    U1 = Int64[]
    V1 = Int64[]
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    wi = U_weights[1:lastidpos,i]
    wj = V_weights[1:lastidpos,i]
    weights[i] += sum(wi.*wj)
    wi = U_weights[nU-lneg:nU,i]
    wj = V_weights[nV-lneg:nV,i]
    weights[i] += sum(wi.*wj)

    ei = U_sortperm[1:lastidpos,i]
    ej = V_sortperm[1:lastidpos,i]
    append!(U1,ei)
    append!(V1,ej)
    ei = U_sortperm[nU-lneg:nU,i]
    ej = V_sortperm[nV-lneg:nV,i]
    append!(U1,ei)
    append!(V1,ej)

    push!(Umatching,U1)
    push!(Vmatching,V1)
  end

  for i = 1:size(U,2)
    ei = Umatching[i]
    ej = Vmatching[i]

    for j = 1:size(U,2)
      num = weights[j]
      ui = U[:,j]
      vi = V[:,j]
      ui = ui[ei]
      vi = vi[ej]
      den = sum(ui.*vi)
      f = num/den
      di[i] = max(di[i],f)
    end
  end

  D = minimum(di)
  Dind = indmin(di)

  utouse = Umatching[Dind]
  vtouse = Vmatching[Dind]

  return utouse,vtouse,D
end

function newbound_rounding_max_approx(U,V)
  nU = size(U,1)
  nV = size(V,1)
  n = size(U,1)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)

  r = size(U,2)
  assert(r==size(V,2))

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

  weights = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    wi = U_weights[1:lastidpos,i]
    wj = V_weights[1:lastidpos,i]
    weights[i] += sum(wi.*wj)
    wi = U_weights[nU-lneg:nU,i]
    wj = V_weights[nV-lneg:nV,i]
    weights[i] += sum(wi.*wj)
  end

  di = zeros(Float64,size(U,2))
  for i = 1:size(U,2)

    P = spzeros(nU,nV)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei = U_sortperm[1:lastidpos,i]
    ej = V_sortperm[1:lastidpos,i]
    P = P + sparse(ei,ej,1,nU,nV)
    ei = U_sortperm[nU-lneg:nU,i]
    ej = V_sortperm[nV-lneg:nV,i]
    P = P + sparse(ei,ej,1,nU,nV)

    for j = 1:size(U,2)
      num = weights[j]
      ui = U[:,j]
      vi = V[:,j]
      M = ui*vi'
      den = sum(spones(P).*M)
      f = num/den
      di[i] = max(di[i],f)
    end
  end

  D = minimum(di)
  Dind = indmin(di)

  ui = U_weights[:,Dind]
  vi = V_weights[:,Dind]

  lastid_ui = findfirst(ui.<0)
  lastid_vi = findfirst(vi.<0)

  if lastid_ui == 0 && lastid_vi == 0
    lastidpos = d
    lneg = -1
  elseif lastid_vi == 0
    lastidpos = min(d,lastid_ui-1)
    lneg = -1
  elseif lastid_ui == 0
    lastidpos = min(d,lastid_vi-1)
    lneg = -1
  else
    lastidpos = min(lastid_ui,lastid_vi)-1
    lneg = min(nU-lastid_ui,nV-lastid_vi)
  end

  ei = U_sortperm[1:lastidpos,Dind]
  ej = V_sortperm[1:lastidpos,Dind]
  ei2 = U_sortperm[nU-lneg:nU,Dind]
  ej2 = V_sortperm[nV-lneg:nV,Dind]

  utouse = vcat(ei,ei2)
  vtouse = vcat(ej,ej2)

  return utouse,vtouse,D
end

function newbound_rounding(U,V)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)
  nU = size(U,1)
  nV = size(V,1)
  r = size(U,2)
  assert(r==size(V,2))

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

  P = spzeros(nU,nV)
  allrecoveries = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei1 = U_sortperm[1:lastidpos,i]
    ej1 = V_sortperm[1:lastidpos,i]
    ei2 = U_sortperm[nU-lneg:nU,i]
    ej2 = V_sortperm[nV-lneg:nV,i]
    ei = vcat(ei1,ei2)
    ej = vcat(ej1,ej2)
    # allrecoveries[i] = evaluate_erdosreyni_experiment(A,B,ei,ej)
    P = P + sparse(ei,ej,1,nU,nV) #+ generate_b_match_overlapping(ei,ej,2,nU,nV)

    # with bmatching: moved to newbound_rounding(U,V,bmatchval)
    # bmatchval = 20
    # for bm = 1:bmatchval
    # if !isempty(ej)
    #   shift!(ej)
    #   P = P + sparse(ei[1:end-bm],ej,1,nU,nV)
    # end
    # end

  end
  @show nnz(P)/prod(size(P))
  return P
end

function newbound_rounding_lowrank_evaluation(U,V)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)
  nU = size(U,1)
  nV = size(V,1)
  r = size(U,2)
  assert(r==size(V,2))

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

#   P = spzeros(nU,nV)
  U1 = []
  V1 = []
  allrecoveries = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei1 = U_sortperm[1:lastidpos,i]
    ej1 = V_sortperm[1:lastidpos,i]
    ei2 = U_sortperm[nU-lneg:nU,i]
    ej2 = V_sortperm[nV-lneg:nV,i]
    ei = vcat(ei1,ei2)
    ej = vcat(ej1,ej2)
    # allrecoveries[i] = evaluate_erdosreyni_experiment(A,B,ei,ej)
#     P = P + sparse(ei,ej,1,nU,nV)
    append!(U1,ei)
    append!(V1,ej)
  end

  all_matches = [U1 V1]
  unique_matches = unique(all_matches,1)
  U1unique = unique_matches[:,1]
  V1unique = unique_matches[:,2]
  uo = U[U1unique,:]
  vo = V[V1unique,:]
  weights = vec(sum(uo.*vo,2))
  X = sparse(U1unique,V1unique,weights,nU,nV)

  return X
end

function newbound_rounding_lowrank_evaluation_relaxed(U,V)
  U_sortperm = sortcolsperm(U,true)
  V_sortperm = sortcolsperm(V,true)
  nU = size(U,1)
  nV = size(V,1)
  r = size(U,2)
  assert(r==size(V,2))

  d = min(size(U_sortperm,1),size(V_sortperm,1))
  U_weights = sort(U,1,rev=true)
  V_weights = sort(V,1,rev=true)

#   U_weights = U_weights[1:d,:]
#   V_weights = V_weights[1:d,:]
#
#   U_sortperm = U_sortperm[1:d,:]
#   V_sortperm = V_sortperm[1:d,:]

#   P = spzeros(nU,nV)
  U1 = []
  V1 = []
  allrecoveries = zeros(size(U,2))
  for i = 1:size(U,2)
    ui = U_weights[:,i]
    vi = V_weights[:,i]
    lastid_ui = findfirst(ui.<0)
    lastid_vi = findfirst(vi.<0)

    if lastid_ui == 0 && lastid_vi == 0
      lastidpos = d
      lneg = -1
    elseif lastid_vi == 0
      lastidpos = min(d,lastid_ui-1)
      lneg = -1
    elseif lastid_ui == 0
      lastidpos = min(d,lastid_vi-1)
      lneg = -1
    else
      lastidpos = min(lastid_ui,lastid_vi)-1
      lneg = min(nU-lastid_ui,nV-lastid_vi)
    end

    ei1 = U_sortperm[1:lastidpos,i]
    ej1 = V_sortperm[1:lastidpos,i]
    ei2 = U_sortperm[nU-lneg:nU,i]
    ej2 = V_sortperm[nV-lneg:nV,i]
    ei = vcat(ei1,ei2)
    ej = vcat(ej1,ej2)
    # allrecoveries[i] = evaluate_erdosreyni_experiment(A,B,ei,ej)
#     P = P + sparse(ei,ej,1,nU,nV)

    # P = P + sparse(ei,ej,1,nU,nV) #+ generate_b_match_overlapping(ei,ej,2,nU,nV)
    # with bmatching:
    append!(U1,ei)
    append!(V1,ej)
    bmatchval = 6
    println("bmatchval is $bmatchval")
    for bm = 1:bmatchval
      if !isempty(ej)
        shift!(ej)
        # P = P + sparse(ei[1:end-bm],ej,1,nU,nV)
        append!(U1,ei[1:end-bm])
        append!(V1,ej)
      end
    end
  end

  all_matches = [U1 V1]
  unique_matches = unique(all_matches,1)
  U1unique = unique_matches[:,1]
  V1unique = unique_matches[:,2]
  uo = U[U1unique,:]
  vo = V[V1unique,:]
  weights = vec(sum(uo.*vo,2))
  X = sparse(U1unique,V1unique,weights,nU,nV)

  return X
end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME