text:
using NumericalIntegration

using LinearAlgebra
using BenchmarkTools
using Statistics
using PyPlot
pygui(true);

function mathieu(z, zτ, t, κ, δ, ϵ, b, τ, ω)
        dzdt=zeros(ComplexF64,size(z))
                dzdt[1]=z[2]
                dzdt[2]=-(δ+(ϵ * cos(ω * t)))*z[1] - κ * z[2] + b * zτ[1]
                return dzdt
end

function mathieu!(dzdt, z, zτ, t, κ, δ, ϵ, b, τ, ω)
        dzdt=zeros(ComplexF64,size(z))
                dzdt[1]=z[2]
                dzdt[2]=-(δ+(ϵ * cos(ω * t)))*z[1] - κ * z[2] + b * zτ[1]
                # return dzdt
end

#integration

function Sinterp(Sloc::Array{ComplexF64,2},indf::Float64)
        i_sub=Int64(floor(indf))::Int64
        remi=(indf-i_sub)::Float64
        # if abs(remi)<(eps(typeof(remi))*5.0)
        #         return S[:,i_sub] .* (1-remi)
        # else
        return Sloc[:,i_sub] .* (1.0-remi) .+ Sloc[:,i_sub+1]*remi
        # end
end
function Sinterp(Sloc,indf::Float64)
        i_sub=Int64(floor(indf))::Int64
        remi=(indf-i_sub)::Float64
        # if abs(remi)<(eps(typeof(remi))*5.0)
        #         return S[:,i_sub] .* (1-remi)
        # else
        return Sloc[:,i_sub] .* (1.0-remi) .+ Sloc[:,i_sub+1]*remi
        # end
end

function Sinterp(S::Array{ComplexF64,2},i::Integer)
        return S[:,i]
end
function Sinterp(S,i::Integer)
        return S[:,i]
end

# S=rand(ComplexF64,1,n+1+m)
# @time Sinterp(S,4.3)
# @time Sinterp(S,4)
# @time S[:,4]
# figure(1)
# plot(1:n+1+m,S[1,:])
# v=1:0.1:n+1+m
# scatter(v,[Sinterp(S[:,:,1],i) for i in v])

# n=2^5
# m=Int64(n*(1 / ω))
# RKtype="RK4"
function mathiestab( κ, δ, ϵ, b, τ, ω, n,m,mult,iter,RKtype::String)::Array{ComplexF64,1}

dt=τ/n
S=rand(ComplexF64,2,n+1+m,mult)
λj=Array{ComplexF64}(undef,mult)
#
# k=Array{ComplexF64}(undef,2)
# k1=Array{ComplexF64}(undef,2)
# k2=Array{ComplexF64}(undef,2)
# k3=Array{ComplexF64}(undef,2)
# k4=Array{ComplexF64}(undef,2)

for kinter in 1:iter
        #Threads.@threads
        for mi in 1:mult
                for kt in n+1:n+m
                if RKtype=="EE"
                        @views k1=mathieu(S[:,kt,mi], Sinterp(S[:,:,mi],kt-n),(kt-(n+1))*dt, κ, δ, ϵ, b, τ, ω)
                        k=k1 ./ 1.0
                # elseif RKtype=="RK2-fast"
                #         @views k1=mathieu(S[:,kt,mi]      ,S[:,kt-n  ,mi],(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                #         @views k2=mathieu(S[:,kt,mi]+k1*dt,S[:,kt-n+1,mi],(kt-(n+1)+1)*dt, κ, δ, ϵ, b, τ, ω)
                #          k=(k1+k2) ./ 2.0
                elseif RKtype=="RK2-Heun"
                        @views k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)  ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                        @views k2=mathieu(S[:,kt,mi]+k1*dt,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                        k=(k1+k2) ./ 2.0
                elseif RKtype=="RK2-Ralston"
                        @views k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)  ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                        @views k2=mathieu(S[:,kt,mi]+k1*dt*2.0/3.0,Sinterp(S[:,:,mi],kt-n+2.0/3.0),(kt-(n+1)+2.0/3.0)*dt, κ, δ, ϵ, b, τ, ω)
                        k=(k1+ 3* k2) ./ 4.0
                elseif RKtype=="RK3"
                        @views k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)  ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                        @views k2=mathieu(S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                        @views k3=mathieu(S[:,kt,mi]-k1*dt*1.0+k2*dt*2.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                        k=(k1+4*k2+k3) ./ 6.0
                elseif RKtype=="RK3-strong"
                        @views k1=mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)  ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                        @views k2=mathieu(S[:,kt,mi]+k1*dt*1.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                        @views k3=mathieu(S[:,kt,mi]+k1*dt*0.25+k2*dt*0.25,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                        k=(k1+k2+4*k3) ./ 6.0
                elseif RKtype=="RK4"
                         @views k1 = mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                         @views k2 = mathieu(S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                         @views k3 = mathieu(S[:,kt,mi]+k2*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                         @views k4 = mathieu(S[:,kt,mi]+k3*dt*1.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                          # mathieu!(k1,S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                          # mathieu!(k2,S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                          # mathieu!(k3,S[:,kt,mi]+k2*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                          # mathieu!(k4,S[:,kt,mi]+k3*dt*1.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                         k = (k1+2*k2+2*k3+k4) ./ 6.0
                 elseif RKtype=="RK5"
                          @views k1 = mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                          @views k2 = mathieu(S[:,kt,mi]+k1*dt*0.25,Sinterp(S[:,:,mi],kt-n+0.25),(kt-(n+1)+0.25)*dt, κ, δ, ϵ, b, τ, ω)
                          @views k3 = mathieu(S[:,kt,mi]+k1*dt*0.125+k2*dt*0.125,Sinterp(S[:,:,mi],kt-n+0.25),(kt-(n+1)+0.25)*dt, κ, δ, ϵ, b, τ, ω)
                          @views k4 = mathieu(S[:,kt,mi]-k2*dt*0.5+k3*dt,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                          @views k5 = mathieu(S[:,kt,mi]+k1*dt*3.0/16.0+k4*dt*9.0/16.0,Sinterp(S[:,:,mi],kt-n+0.75),(kt-(n+1)+0.75)*dt, κ, δ, ϵ, b, τ, ω)
                          @views k6 = mathieu(S[:,kt,mi]-k1*dt*3.0/7.0+k2*dt*2.0/7.0+k3*dt*12.0/7.0-k4*dt*12.0/7.0+k5*dt*8.0/7.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                           # mathieu!(k1,S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                           # mathieu!(k2,S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                           # mathieu!(k3,S[:,kt,mi]+k2*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                           # mathieu!(k4,S[:,kt,mi]+k3*dt*1.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                          k = (7*k1    +32*k3+12*k4+32*k5+7*k6) ./ 90.0
                  elseif RKtype=="RK5-2"
                           @views k1 = mathieu(S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                           @views k2 = mathieu(S[:,kt,mi]+k1*dt*1/5,Sinterp(S[:,:,mi],kt-n+1/5),(kt-(n+1)+1/5)*dt, κ, δ, ϵ, b, τ, ω)
                           @views k3 = mathieu(S[:,kt,mi]+k2*dt*2/5,Sinterp(S[:,:,mi],kt-n+2/5),(kt-(n+1)+2/5)*dt, κ, δ, ϵ, b, τ, ω)
                           @views k4 = mathieu(S[:,kt,mi]+k1*dt*9/4-k2*dt*5+k3*dt*15/4,Sinterp(S[:,:,mi],kt-n+1),(kt-(n+1)+1)*dt, κ, δ, ϵ, b, τ, ω)
                           @views k5 = mathieu(S[:,kt,mi]-k1*dt*63/100+k2*dt*9/5-k3*dt*13/20+k4*dt*2/25,Sinterp(S[:,:,mi],kt-n+3/5),(kt-(n+1)+3/5)*dt, κ, δ, ϵ, b, τ, ω)
                           @views k6 = mathieu(S[:,kt,mi]-k1*dt*6/25+k2*dt*4/5+k3*dt*2/15+k4*dt*8/75,Sinterp(S[:,:,mi],kt-n+4/5),(kt-(n+1)+4/5)*dt, κ, δ, ϵ, b, τ, ω)
                            # mathieu!(k1,S[:,kt,mi]      ,Sinterp(S[:,:,mi],kt-n)        ,(kt-(n+1)  )*dt, κ, δ, ϵ, b, τ, ω)
                            # mathieu!(k2,S[:,kt,mi]+k1*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                            # mathieu!(k3,S[:,kt,mi]+k2*dt*0.5,Sinterp(S[:,:,mi],kt-n+0.5),(kt-(n+1)+0.5)*dt, κ, δ, ϵ, b, τ, ω)
                            # mathieu!(k4,S[:,kt,mi]+k3*dt*1.0,Sinterp(S[:,:,mi],kt-n+1.0),(kt-(n+1)+1.0)*dt, κ, δ, ϵ, b, τ, ω)
                           k = (17/144*k1+25/36*k3+1/72*k4-25/72*k5+25/48*k6)
                else
                        println("RKtype")
                        error("nincs ilyen integrator")
                end
                 S[:,kt+1,mi] = S[:,kt,mi]+dt*k;
                end
        end

                Sj=reshape(S[:,1:n+1,:],(2*(n+1),mult))
                Vj=reshape(S[:,end-n:end,:],(2*(n+1),mult))
                        # Sj=S[:,1:n+1,1]
                        # Vj=S[:,end-n:end,1]

        # SS=[integrate(1:size(Sj,1), Sj[:,m1] .* Sj[:,m2], SimpsonEvenFast()) for m1 in 1:mult, m2 in 1:mult]
        # SV=[integrate(1:size(Sj,1), Sj[:,m1] .* Vj[:,m2], SimpsonEvenFast()) for m1 in 1:mult, m2 in 1:mult]

        # SS=[integrate(1:size(Sj,1), Sj[:,m1] .* Sj[:,m2], TrapezoidalEvenFast()) for m1 in 1:mult, m2 in 1:mult]
        # SV=[integrate(1:size(Sj,1), Sj[:,m1] .* Vj[:,m2], TrapezoidalEvenFast()) for m1 in 1:mult, m2 in 1:mult]
                # println(n)
                # println(size(Sj,1))
                # SS=[integrate(1:size(Sj,1), Sj[:,m1] .* Sj[:,m2], RombergEven()) for m1 in 1:mult, m2 in 1:mult]
                # SV=[integrate(1:size(Sj,1), Sj[:,m1] .* Vj[:,m2], RombergEven()) for m1 in 1:mult, m2 in 1:mult]

                        SS=(transpose(Sj)*Sj)
                        SV=(transpose(Sj)*Vj)

        #Hj=(Sj\Vj)
        Hj=(SS\SV)
        # Hj=((transpose(Sj)*Sj)\(transpose(Sj)*Vj))
        (λj,Gj)=eigen(Hj)
                # end

        S[:,1:n+1,:] = (Vj*(Gj/diagm(λj)))[:]
        # S[:,1:n+1,1] = (Vj*(Gj/diagm(λj)))[:]
        # S[:,1:n+1,:] = (Vj*(Gj))[:]

        # figure(3)
        # clf()
        # plot((1:n+1+m).*dt,S[2,1:n+1+m,:])
        # plot((1:n+1).*dt,S[2,1:n+1,:])
        # pause(0.01)
        #
        # figure(4)
        # scatter(kinter .* ones(mult) .+ 0.1,abs.(λj))
        # pause(0.01)
        #
        # figure(8)
        # clf()
        # surf(abs.(Gj))
end

return λj
end

println("----------=======<<<<<<<<<<<<<#######   START   #######>>>>>>>>>>>>>>========----------")
 const κ=0.05
 const δ=5
 const ϵ=1.0
 const b=0.1
 const τ=2pi
 const ω=0.25

 mult=8 # number of computed eigenvalues
 iter=8 # number of iteration in issi

@time mathiestab( κ, δ, ϵ, b, τ, ω, 2^5, Int64(2^5*(1 / ω)),mult,iter,"RK4")
 #
 # for NN in 5:12
 # @time mathiestab( κ, δ, ϵ, b, τ, ω, 2^NN, Int64(2^NN*(1 / ω)),mult,iter,"RK4")
 # end
 #
 # @btime mathiestab( κ, δ, ϵ, b, τ, ω, 2^6, Int64(2^6*(1 / ω)),mult,iter,"RK4") #First run "time" is too long

# figure(4)
# clf()
@elapsed mathiestab( κ, δ, ϵ, b, τ, ω, 2^5, Int64(2^5*(1 / ω)),mult,iter,"RK2-Ralston")
@debug mathiestab( κ, δ, ϵ, b, τ, ω, 2^5, Int64(2^5*(1 / ω)),mult,iter,"RK2-Ralston") #First run "time" is too long
# @elapsed mu_ref=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω,10^6,m,mult,iter,"RK4")))
# @code_warntype mathiestab( κ, δ, ϵ, b, τ, ω, 2^4, Int64(2^4*(1 / ω)),mult,iter,"RK4") #First run "time" is too long

println("Computing the reference")
# n=2^19 # 16 sec
# n=2^19#16
n=2^15#16
m=Int64(n*(1 / ω))
@time mu_ref=mathiestab( κ, δ, ϵ, b, τ, ω, n,m,mult,iter,"RK4")
if imag(mu_ref[end])<0
        conj!(mu_ref)
end

println("--------------------===================------------------------------")
integrators1=["EE"]
integrators2=["RK2-Ralston","RK2-Heun"]
integrators3=["RK3","RK3-strong"]
integrators4=["RK4"]
integrators5=["RK5","RK5-2"]
integrators=[integrators1...,integrators2...,integrators3...,integrators4...,integrators5...]

# ni=2 .^ (4:19)
ni=2 .^ (2:13)
mumax=zeros(ComplexF64,length(ni),length(integrators),mult)
timelength=zeros(length(ni),length(integrators))
timestd=zeros(length(ni),length(integrators))

# figure(1)
# clf()
# figure(2)
# clf()
# figure(3)
# clf()

figure(1)
# clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("n / error")
xlabel("n")
ylabel("error")
#
figure(2)
# clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("n / time")
xlabel("n")
ylabel("time")
#
figure(3)
# clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("time / error")
xlabel("time")
ylabel("error")

Nsample=1

for kn in 1:length(ni)
        n=ni[kn]
        m=Int64(n*(1 / ω))

        for ksample in 1:Nsample#maximum([Nsample-Int(log(ni[kn])/log(2)),1])
                println("<<<<<<<<<<<<<<>>>>>>>>>>>>>>>>>>>>")
                println([ksample,kn,ni[kn]])
                for kRKi in 1:length(integrators)
                        integratortype=integrators[kRKi]
                        print(integratortype)
                        print("  -  ")

        # timeing = @benchmark (mumax[kn]=maximum(abs.(mathiestab( κ, δ, ϵ, b, τ, ω,n,m,mult,iter,"RK4"))))
        # timelength[kn]=mean(timeing.times)
        # timestd[kn]=std(timeing.times)
        timeing = @elapsed (mui_set=mathiestab( κ, δ, ϵ, b, τ, ω,n,m,mult,iter,integratortype))
        if imag(mui_set[end])<0
                conj!(mui_set)
        end
        mumax[kn,kRKi,:] .= mui_set
        timelength[kn,kRKi]+=timeing
        println(timeing)

        figure(1)
        scatter(n,abs(mumax[kn,kRKi,end] - mu_ref[end]))
        # scatter(n,abs(mumax[kn,kRKi,end]) - abs(mu_ref[end]))
        figure(2)
        scatter(n,timeing)
        figure(3)
        scatter(timeing,abs(mumax[kn,kRKi,end] - mu_ref[end]))
        # scatter(timeing,abs(mumax[kn,kRKi,end]) - abs(mu_ref[end]))

        # figure(5)
        # scatter(abs(mumax[kn,kRKi,end] - mu_ref[end]),abs(mumax[kn,kRKi,end]) - abs(mu_ref[end]))
                end
        end

        figure(1)
        plot(ni[1:kn,:],abs.(mumax[1:kn,:,end] .- mu_ref[end]))
        savefig("1_v_nt")
        # savefig("1_t")

        figure(2)
        plot(ni[1:kn,:],timelength[1:kn,:]./Nsample)
        savefig("2_v_nt")
        # savefig("2_t")

        figure(3)
        plot(timelength[1:kn,:]./Nsample,abs.(mumax[1:kn,:,end] .- mu_ref[end]))
        savefig("3_v_nt")
        # savefig("3_t")

        # figure(5)
        # # clf()
        # ax=gca()
        # ax.set_xscale("log")
        # ax.set_yscale("log")
        # grid("on")
        # savefig("5_nt")
end

figure(1)
plot(ni,abs.(mumax[:,:,end] .- mu_ref[end]))
savefig("1_v_nt")
# savefig("1_t")

figure(2)
plot(ni,timelength./Nsample)
savefig("2_v_nt")
# savefig("2_t")

figure(3)
plot(timelength./Nsample,abs.(mumax[:,:,end] .- mu_ref[end]))
savefig("3_v_nt")
# savefig("3_t")

# figure(5)
# savefig("5_nt")
println("---------   <<<<<<<<  END >>>>>>>>>   ----------")
                # figure(1)
                # plot(log.(ni)./log(2),log.(abs.(mumax .- mu_ref))./log(2))
                # figure(2)
                # plot(log.(ni)./log(2),log.(timelength)./log(2))
                # # plot(log.(timelength .+ timestd))
                # # plot(log.(timelength .- timestd))

# #----------------------------
# using MDBM
#
# foo=(x,y)->(maximum(abs.(mathiestab( κ, x, ϵ, y, τ, ω,n,m,mult,iter,"RK4")))-1)
#
# δv=-0.1:0.2:10
# # ϵv=-3:0.1:3
# bv=-1:0.2:1
#
#
# mymdbm=MDBM_Problem(foo,Axis.([δv,bv]))
# MDBM.solve!(mymdbm,3)
#
# x_eval,y_eval=getevaluatedpoints(mymdbm)
# x_sol,y_sol=getinterpolatedsolution(mymdbm)
#
# fig = figure(11);clf()
# PyPlot.scatter(x_eval,y_eval,s=5)
# PyPlot.scatter(x_sol,y_sol)
#
# println("cuccma")

figure(11)
 clf()
ax=gca()
ax.set_xscale("log")
ax.set_yscale("log")
grid("on")
title("n / error")
xlabel("n")
ylabel("error")

A2=hcat([-1,1,-1,1],[1,1,1,1],[3,-2,1,0],[3,2,1,0])'
A1=hcat([-1,1],[1,1])'

w2=[0,2/3,0,2]'*inv(A2)
w1=[0,2]'*inv(A1)

# t=0:0.1:pi/2
# y=sin.(try)
# dydy=cos.(t)
for kk in 3:20
        dt=pi/(2.0 ^ (kk))

        t=0:dt:pi/2
        # y=t
        # dydt=1 .+ 0 .* (t)
        y=sin.(t)
        dydt=cos.(t)

        # # y=t
        # # dydt=1 .+ 0 .* (t)
        # y=t.^7.0
        # dydt=0.0 .+ 7.0 .* t .^6.0
        I=0.0
        for l in 1:(length(t)-1)
                # global I
                # I += dt*y[l]
                # I += dt*y[l+1]
                 # I += dt/2*(y[l]+y[l+1])
                  # I += (dt/2)*w1*[y[l],y[l+1]]
                I += (dt/2)*w2*[y[l],y[l+1],dydt[l]*(dt/2),dydt[l+1]*(dt/2)]
        end

        figure(11)
        scatter(length(t),abs(I-1.0))
        println(I)

        IS=integrate(t, y, SimpsonEvenFast())
        IT=integrate(t, y, TrapezoidalEvenFast())
        IR=integrate(t, y, RombergEven())
        scatter(length(t),abs(sum(y)*dt-1.0), color=:green)
        scatter(length(t),abs(IS-1.0), color=:red)
        scatter(length(t),abs(IT-1.0), color=:blue)
        scatter(length(t),abs(IR-1.0), color=:black)
end

length(dt)


text:
# script that tests all packages and writes the results to a database file

using DataFrames
using NewPkgEval
using SQLite
using Dates

function main(;julia_releases=["1.2"], pkg_names=["Example"], registry="General",
               dbfile=nothing)
    db = if dbfile === nothing
        # in memory
        SQLite.DB()
    else
        SQLite.DB(dbfile)
    end

    # get the Julia versions
    julia_versions = Dict(NewPkgEval.download_julia(julia_release) => julia_release
                          for julia_release in julia_releases)
    NewPkgEval.prepare_julia.(keys(julia_versions))

    # get the packages
    NewPkgEval.prepare_registry(update=true)
    pkgs = NewPkgEval.read_pkgs(pkg_names)

    NewPkgEval.prepare_runner()

    # prepare the database
    SQLite.execute!(db, """
        CREATE TABLE IF NOT EXISTS builds (package_name TEXT,
                                           package_version TEXT,
                                           julia_release TEXT,
                                           julia_version TEXT,
                                           run INT,
                                           status TEXT,
                                           reason TEXT,
                                           log TEXT,
                                           datetime TEXT,
                                           duration REAL)
        """)

    # find a unique run identifier
    run = (SQLite.Query(db, "SELECT COALESCE(MAX(run), 0) FROM builds") |> DataFrame)[1,1] + 1

    function store_result(julia_version, package_name, package_version, start, status, reason, log)
        stop = now()
        elapsed = (stop-start) / Millisecond(1000)
        julia_release = julia_versions[julia_version]

        # stringify all values to prevent serialization, but keep `missing`s
        string_or_missing(obj) = ismissing(obj) ? missing : string(obj)

        SQLite.Query(db, "INSERT INTO builds VALUES (?, ?, ?, ?, ?, ?, ?, ?, ?, ?)";
                        values=string_or_missing.([package_name, package_version,
                                                julia_release, julia_version,
                                                run, status, reason,
                                                log, now(), elapsed]))
    end

    # test!
    NewPkgEval.run(collect(keys(julia_versions)), pkgs; callback=store_result)

    return
end

isinteractive() || main(julia_releases=["1.0", "1.2", "latest"], pkg_names=String[],
                        dbfile=joinpath(@__DIR__, "test.db"))

ext: jl


text:
using Discord
include("base.jl")

function la(client::Client,m::Message;over_time="")
    global  management_dict

    channel_id_str = string(m.channel_id)

    if !haskey(management_dict,channel_id_str)
        return reply(client,m,"初期設定が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if !haskey(management_dict[channel_id_str],"remain_attack_channel")
        return reply(client,m,"残凸把握板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end

    if  !haskey(management_dict[channel_id_str],"output_channel")
        return reply(client,m,"予約確認板の登録が完了しておりません．`.set`コマンドにて初期設定を完了させてください")
    end
    if !management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["over"] && length(over_time) == 0
        return reply(client,m,"持ち越しした場合には`.la [持ち越し時間]`を送信して下さい")
    end
    boss_number = string(management_dict[channel_id_str]["boss_number"])
    reserve_list =  management_dict[channel_id_str]["reserve"][boss_number]["reserve_list"]

    dict_number = 0

    #予約ずみかどうか事前に凸宣言をしていたかどうかをチェックする
    for i = 1:length(reserve_list)
        if reserve_list[i]["id"] == m.author.id
            if !reserve_list[i]["attack"]
                return reply(client,m,"凸宣言がされておりません．`.attack`コマンドにて宣言を行って下さい")
            end
            dict_number = i
        end
    end

    #予約していなかっtら怒る
    if dict_number == 0
        return reply(client,m,"凸予約がされていません.")
    end

    next_boss_number = management_dict[channel_id_str]["boss_number"]+1
    management_dict[channel_id_str]["reserve"][boss_number] =  Dict(
        "reserve_list" => [],
        "plan_remain_hp" => get_default_hp(management_dict[channel_id_str]["lap"]+1,management_dict[channel_id_str]["boss_number"]),
        "remain_hp" => get_default_hp(management_dict[channel_id_str]["lap"]+1,management_dict[channel_id_str]["boss_number"])
    )

    if next_boss_number == 6
        next_boss_number =1
        management_dict[channel_id_str]["lap"] += 1
    end

    management_dict[channel_id_str]["boss_number"] =next_boss_number

    output_message_content = create_output_message(management_dict[channel_id_str])
    edit_message(client,management_dict[channel_id_str]["output_channel"],management_dict[channel_id_str]["output_message"],content=output_message_content)

    #次の人に通知を行う
    notify_message = ""

    for next_reserve_dict in management_dict[channel_id_str]["reserve"][string(next_boss_number)]["reserve_list"]
        notify_message *= "<@!$(next_reserve_dict["id"])> "
    end

    if notify_message != ""
        notify_message *= "\n$(Boss_dict["normal"][next_boss_number].name)になりました．"
        reply(client,m,notify_message)
    end

    #残凸状況を修正，持ち越ししていることを把握する
    if management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["remain"] > 0 &&  length(over_time) >0
        management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["over"] = true
        management_dict[channel_id_str]["remain_attack"][string(m.author.id)]["over_content"]  = "$(Boss_dict["normal"][int(boss_number)].name) $(over_time)"
        remain_attack_message_content = make_remain_attack_message_content(management_dict[channel_id_str]["remain_attack"],management_dict[channel_id_str]["today"])
        edit_message(client,management_dict[channel_id_str]["remain_attack_channel"],management_dict[channel_id_str]["remain_attack_message"],content=remain_attack_message_content)
    end

    save_management_dict(management_dict)
    create(client, Reaction, m, '👍')

end

max_stars_repo_name: hibibol/clan-battle

max_issues_repo_name: hibibol/clan-battle

max_forks_repo_name: hibibol/clan-battle


text:
module JOLI4Flux
    using JOLI: joAbstractLinearOperator
    using Zygote
    import Base.*
    using Zygote: @adjoint

    @adjoint a::joAbstractLinearOperator * b::AbstractVecOrMat =
        *(a, b), Δ -> (nothing, transpose(a) * Δ)
end


text:
float1zero = [
:ceil,
:floor,
:isfinite,
:isinf,
:isnan,
:round,
:trunc,
]
for f in float1zero; @eval @zerograd $f(x); end

float2zero = [
:div, # The quotient from Euclidean division. Computes x/y, truncated to an integer; operators
#:÷,  # Same as div.
]
for f in float2zero; @eval @zerograd $f(x1,x2); end

float1arg = [
(:+, :dy),  # supports (N,) (A,) (N,N) (N,A) (A,N) (A,A); float arraymath abstractarraymath operators
(:-, :(-dy)),  # supports (N,) (A,) (N,N) (N,A) (A,N) (A,A); float arraymath
(:float, :dy),          # float
]

for (f,g) in float1arg
    @eval @primitive $f(x),dy,y $g
    addtest1(f,(-Inf,Inf))
end

float2arg = [
(:+, :dy, :dy), # extra (N,) (A,); float,arraymath,abstractarraymath,operators
(:-, :dy, :(-dy)),                    # extra (N,) (A,); float,arraymath
#(:rem, :dy, :(-dy.*trunc(x1./x2))),   # Remainder from Euclidean division, return same sign as x, BUG: missing (A,A); BUG: WARNING: (:check_grads,(:sum,:rem),:args,(-0.13412338383912367,[-0.00025363687477246275,0.4389355563026644]),:exact,(2.0,[-528.0,0.0]),:numeric,(1.9999999999997797,[0.8920182562441314,0.0])); float,arraymath
#(:%, :dy, :(-dy.*trunc(x1./x2))),     # same as rem
#(:mod, :dy, :(-dy.*floor(x1./x2))),   # BUG: WARNING: (:check_grads,(:sum,:mod),:args,([-1.850960311615227,-1.282613199024709],[0.8035558268314972,-0.32067619631949534]),:exact,([1.0,1.0],[3.0,-3.0]),:numeric,([1.0000000000021103,1.0000000000021103],[2.9999999999996696,3203.2619631949538])); Modulus after division, return same sign as y; float,arraymath
#(:(==), 0,                      # BUG: StackOverflowError(). supports any pair of values; float,operators,abstractarray
]

for (f,g1, g2) in float2arg
    @eval @primitive $f(x1,x2),dy,y unbroadcast(x1,$g1) unbroadcast(x2,$g2)
    addtest2(f,(-Inf,Inf))
end

# Methods for multiplication:
# *(x::Float64, y::Float64) at float.jl:212  (same as .*)
# *(A::Number, B::AbstractArray{T,N}) at abstractarraymath.jl:54  (calls .*)
# *(A::AbstractArray{T,N}, B::Number) at abstractarraymath.jl:55  (calls .*)
# The Array-Array case is handled by linalg/matmul.

# To distinguish these from the Array-Array case, we type the scalars:
@primitive (*)(x1::Number,x2::Number),dy,y  dy*x2  dy*x1
@primitive (*)(x1::Number,x2),dy,y  unbroadcast(x1,dy.*x2)  dy*x1
@primitive (*)(x1,x2::Number),dy,y  dy*x2  unbroadcast(x2,dy.*x1)
addtest(*,randn(),randn())
addtest(*,randn(),randn(2))
addtest(*,randn(2),randn())

# Methods for division:
# /(x::Float64, y::Float64) at float.jl:214
# /(A::AbstractArray{T,N}, B::Number) at abstractarraymath.jl:57
# The Array-Array case is handled by linalg/generic.
# There is no Number-Array support.

@primitive (/)(x1,x2::Number),dy,y  (dy/x2)  unbroadcast(x2,-dy.*x1./abs2(x2))
x = randn(); a = randn(2)
addtest(/,randn(),randn())
addtest(/,randn(2),randn())

# These are defined in terms of isless which is handled in interfaces.jl
# float2arg1 = Dict{Symbol,Any}(
# :<, 0,                         # only supports (N,N), arrays not supported; float,operators
# :<=, 0,                        # only supports (N,N), arrays not supported; float,operators
# :>, 0,                         # only supports (N,N), arrays not supported; operators
# :>=, 0,                        # only supports (N,N), arrays not supported; operators
# )

# These operators take bool or int as input and are excluded:

# bool1arg = Dict{Symbol,Any}(
# #:(~), :todo, # bitwise not, domain=Integer; bool,int,gmp,bitarray,arraymath
# #:(!), :todo, # boolean not, domain=Bool; bool,bitarray,arraymath
# )

# int1arg = Dict{Symbol,Any}(
# #:(~), :todo, # bitwise not, domain=Integer; bool,int,gmp,bitarray,arraymath
# )
# int2arg = Dict{Symbol,Any}(
# #:$, :todo,                     # domain: Integers, bitwise xor; int,bool,bitarray,gmp,operators,promotion,arraymath,sparsematrix
# #:&, :todo,                     # domain: Integers, bitwise and
# #:|, :todo,                     # domain: Integers, bitwise or
# )

# TODO:

# eval
# convert
# promote_rule
# float
# unsafe_trunc
# trunc
# floor
# ceil
# round
# widen
# -
# +
# *
# /
# muladd
# rem
# cld
# mod
# ==
# !=
# <
# <=
# isequal
# isless
# cmp
# abs
# isnan
# isfinite
# isinf
# hx: Not exported
# hash
# precision
# float_lex_order: Not exported
# nextfloat
# prevfloat
# issubnormal
# typemin
# typemax
# realmin
# realmax
# eps
# bswap
# reinterpret
# sign_mask: Not exported
# exponent_mask: Not exported
# exponent_one: Not exported
# exponent_half: Not exported
# significand_mask: Not exported
# significand_bits: Not exported
# exponent_bits: Not exported
# exponent_bias: Not exported
# $(Expr(:$, :fn)): Not a symbol
# big
# (abs, :(dy.*sign(x))) moved to number.jl


text:
function analyzer(mat::Dict, outdir::AbstractString; dx=1, makeplots=true, isExt=false)

  R1map = mat["R10"]
  S0map = mat["S0"]
  modelmap = mat["modelmap"]
  Ct = mat["Ct"]
  Kt = mat["Kt"]
  ve = mat["ve"]
  vp = mat["vp"]
  resid = mat["resid"]
  q = quantile(S0map[:], 0.99)
  S0map[S0map .> q] .= q
  back = (S0map .- minimum(S0map)) / (maximum(S0map) - minimum(S0map))
  mask = convert(Array{Bool,2}, mat["mask"])

  # compare to known truths
  if ~isExt
    u = ones(div(10,dx),div(10,dx))
    Kt_truth = repeat([0.01*u; 0.02*u; 0.05*u; 0.1*u; 0.2*u; 0.35*u], 1, 5)
    ve_truth = repeat([0.01*u; 0.05*u; 0.1*u; 0.2*u; 0.5*u]', 6, 1)
  else
    u = ones(div(180,dx),div(10,dx))
    Kt_truth = hcat(0.01u, 0.02u, 0.05u, 0.1u, 0.2u)

    u = ones(div(10,dx),div(50,dx))
    u = vcat(0.1u, 0.2u, 0.5u)
    ve_truth = vcat(u, u, u, u, u, u)

    u = ones(div(30,dx),div(50,dx))
    vp_truth = vcat(0.001u, 0.005u, 0.01u, 0.02u, 0.05u, 0.1u)
  end

  Kt_error = clamp.(100.0*(Kt - Kt_truth) ./ (Kt_truth .+ eps()), -100.0, 100.0)
  ve_error = clamp.(100.0*(ve - ve_truth) ./ (ve_truth .+ eps()), -100.0, 100.0)
  cccKt = ccc(Kt_truth, Kt)
  cccve = ccc(ve_truth, ve)
  cccvp = 0.0 # ccc for ExtTofts is calculated in next block
  printstyled("Kt\n\tRMSE:\t$(sqrt(norm(Kt_error)^2 / length(Kt_error))) %\n", color=:green)
  printstyled("\terrmax:\t$(maximum(abs.(Kt_error)))\n", color=:green)
  printstyled("\tCCC:\t$cccKt\n", color=:green)
  printstyled("ve\n\tRMSE:\t$(sqrt(norm(ve_error)^2 / length(ve_error))) %\n", color=:green)
  printstyled("\terrmax:\t$(maximum(abs.(ve_error)))\n", color=:green)
  printstyled("\tCCC:\t$cccve\n", color=:green)

  if isExt
    vp_error = clamp.(100.0*(vp - vp_truth) ./ (vp_truth .+ eps()), -100.0, 100.0)
    cccvp = ccc(vp_truth, vp)
    printstyled("vp\n\tRMSE:\t$(sqrt(norm(vp_error)^2 / length(vp_error))) %\n", color=:green)
    printstyled("\terrmax:\t$(maximum(abs.(vp_error))) %\n", color=:green)
    printstyled("\tCCC:\t$cccvp\n", color=:green)
  end

  if !makeplots
    return (cccKt, cccve, cccvp)
  end

  if ~isExt
    ytpos = collect((0+floor(Integer, 5/dx)):div(10,dx):(div(60,dx)-1))
    xtpos = collect((0+floor(Integer, 5/dx)):div(10,dx):(div(50,dx)-1))
    ytlabels = [string(x) for x in [0.01,0.02,0.05,0.1,0.2,0.35]]
    xtlabels = [string(x) for x in [0.01,0.05,0.1,0.2,0.5]]
    # Size of figures that contain 2D maps and their x/y labels
    mapWidth = 4.5
    mapHeight = 4.5
    mapLabelX = "\$v_\\mathrm{e}\$"
    mapLabelY = "\$K^\\mathrm{trans}\$"
  else
    ytpos = collect((div(10,dx)+floor(Integer, 5/dx)):div(30,dx):(div(180,dx)-1))
    xtpos = collect((0+floor(Integer, 5/dx)):div(10,dx):(div(50,dx)-1))
    ytlabels = [string(x) for x in [0.001, 0.005, 0.01, 0.02, 0.05, 0.1]]
    xtlabels = [string(x) for x in [0.01,0.02,0.05,0.1,0.2]]
    mapWidth = 3.5
    mapHeight = 6
    mapLabelX = "\$K^\\mathrm{trans}\$"
    mapLabelY = "\$v_\\mathrm{p}\$"
  end

  println("Plotting results ...")

  # AIF
  figure(figsize=(4.5,4.5))
  clf()
  plot(mat["t"], mat["Cp"], "ko-")
  xlabel("time (min)")
  yticks(collect(0:2:10)) # This produces an error
  ylim(0,10)
  ylabel("[Gd-DTPA] (mM)")
  title("arterial input function, \$C_p\$")
  savefig("$outdir/aif.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(modelmap, interpolation="nearest", cmap="cubehelix")
  title("model used")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/modelmap.pdf",bbox_inches="tight")

  # PARAMETER MAPS
  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(Kt, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=maximum(Kt_truth))
  title("\$K^\\mathrm{trans}\$ (min\$^{-1}\$)")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/Kt.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(ve, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=maximum(ve_truth))
  title("\$v_\\mathrm{e}\$")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/ve.pdf",bbox_inches="tight")

  if isExt
    figure(figsize=(mapWidth, mapHeight))
    clf()
    imshow(vp, interpolation="nearest", cmap="cubehelix", vmin=0, vmax=maximum(vp_truth))
    title("\$v_p\$")
    xticks(xtpos, xtlabels, fontsize=8)
    yticks(ytpos, ytlabels)
    xlabel(mapLabelX)
    ylabel(mapLabelY)
    colorbar()
    savefig("$outdir/vp.pdf",bbox_inches="tight")
  end

  figure(figsize=(mapWidth, mapHeight))
  clf()
  imshow(resid, interpolation="nearest", cmap="cubehelix", vmin=0)
  title("residual")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/resid.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  m = maximum(abs.(Kt_error))
  imshow(Kt_error, interpolation="nearest", cmap="PiYG", vmin=-m, vmax=m)
  title("% error in \$K^\\mathrm{trans}\$")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/Kt_error.pdf",bbox_inches="tight")

  figure(figsize=(mapWidth, mapHeight))
  clf()
  m = maximum(abs.(ve_error))
  imshow(ve_error, interpolation="nearest", cmap="PiYG", vmin=-m, vmax=m)
  title("% error in \$v_\\mathrm{e}\$")
  xticks(xtpos, xtlabels, fontsize=8)
  yticks(ytpos, ytlabels)
  xlabel(mapLabelX)
  ylabel(mapLabelY)
  colorbar()
  savefig("$outdir/ve_error.pdf",bbox_inches="tight")

  if isExt
    figure(figsize=(mapWidth, mapHeight))
    clf()
    m = maximum(abs.(vp_error))
    imshow(vp_error, interpolation="nearest", cmap="PiYG", vmin=-m, vmax=m)
    title("% error in \$v_p\$")
    xticks(xtpos, xtlabels, fontsize=8)
    yticks(ytpos, ytlabels)
    xlabel(mapLabelX)
    ylabel(mapLabelY)
    colorbar()
    savefig("$outdir/vp_error.pdf",bbox_inches="tight")
  end

  return (cccKt, cccve, cccvp)
end

function analyze(n, mat::Dict, outdir::AbstractString; kwargs...)
  if n == 4
    return analyzer(mat, outdir; isExt=true, kwargs...)
  elseif n == 6
    return analyzer(mat, outdir; isExt=false, kwargs...)
  end
end

function validate(n, outdir::AbstractString; kwargs...)
  @assert n == 4 || n == 6 "n must be 4 or 6"
  cd("$(@__DIR__)/../test/q$n")

  println("Running analysis of noise-free QIBA v$n data ...")
  isdir("$outdir/results") || mkdir("$outdir/results")
  results = fitdata(datafile="qiba$n.mat",outfile="$outdir/results/results.mat",save=false)
  ccc = analyze(n, results, "$outdir/results", dx=10; kwargs...)

  println("Running analysis of noisy QIBA v$n data ...")
  isdir("$outdir/results_noisy") || mkdir("$outdir/results_noisy")
  results = fitdata(datafile="qiba$(n)noisy.mat", outfile="$outdir/results_noisy/results.mat",save=false)

  cccnoisy = analyze(n, results, "$outdir/results_noisy"; kwargs...)
  println("Validation complete. If saved, results can be found in $outdir.")
  return (ccc, cccnoisy)
end

validate(n; kwargs...) = validate(n, "$(@__DIR__)/../test/q$n"; kwargs...)
function validate(; kwargs...)
  ccc6, cccnoisy6 = validate(6; kwargs...)
  ccc4, cccnoisy4 = validate(4; kwargs...)
  return (ccc6, cccnoisy6, ccc4, cccnoisy4)
end


text:
"""
$(SIGNATURES)

Solve the Lindblad equation defined by `A` for a total evolution time `tf`.

...
# Arguments
- `A::Annealing`: the Annealing object.
- `tf::Real`: the total annealing time.
- `tspan` = (0, tf): time interval to solve the dynamics.
- `vectorize::Bool = false`: whether to vectorize the density matrix.
- `kwargs` : other keyword arguments supported by `DifferentialEquations`.
...
"""
function solve_lindblad(
    A::Annealing,
    tf::Real;
    tspan=(0.0, tf),
    vectorize::Bool=false,
    kwargs...,
)
    u0 = build_u0(A.u0, :m, vectorize=vectorize)
    L = OpenQuantumBase.lindblad_from_interactions(A.interactions)
    if vectorize
        error("Vectorization is not yet supported for Lindblad equation.")
    end
    f = DiffEqLiouvillian(A.H, [], L, size(A.H, 1))
    p = ODEParams(f, float(tf), A.annealing_parameter)
    prob = ODEProblem(f, u0, tspan, p)
    solve(prob; alg_hints=[:nonstiff], kwargs...)
end

function build_ensemble_lindblad(
    A::Annealing,
    tf::Real,
    output_func,
    reduction;
    tspan=(0.0, tf),
    kwargs...,
)
    u0 = build_u0(A.u0, :v)
    L = DiffEqLiouvillian(A.H, [], OpenQuantumBase.lindblad_from_interactions(A.interactions), size(A.H, 1))
    cb = LindbladJumpCallback()
    p = ODEParams(L, tf, A.annealing_parameter)
    update_func = function (cache, u, p, t)
        update_cache!(cache, p.L, p, t)
    end
    cache = get_cache(A.H)
    diff_op = DiffEqArrayOperator(cache, update_func=update_func)
    jac_cache = similar(cache)
    jac_op = DiffEqArrayOperator(jac_cache, update_func=update_func)
    ff = ODEFunction(diff_op; jac_prototype=jac_op)

    prob = ODEProblem{true}(ff, u0, tspan, p, callback=cb)
    ensemble_prob = EnsembleProblem(
        prob;
        output_func=output_func,
        reduction=reduction,
        kwargs...
    )
    ensemble_prob
end

max_stars_repo_name: huangjingzheng/OpenQuantumTools

max_issues_repo_name: huangjingzheng/OpenQuantumTools

max_forks_repo_name: huangjingzheng/OpenQuantumTools


text:
using Images

struct star
  maxval
  x::Int32
  y::Int32
end

# pseudo-constants; they would be real constants except that they can be altered from the command line

EDGE::Int32 = 25        # size of one cell and of one centroid for isolating stars
                        # This is also the size of the unused margin around the photo (so that centroid seeking doesn't wander off the picture)
CUTMAX::Float16 = 1.00  # don't use maxima above this because they're probaby wrong
CUTMIN::Float16 = 0.01  # don't use maxuma below this because they're probably just noise
STARCOUNT::Int32 = 40   # How many stars are we looking for?
EPSILON::Int32 = 2      # How close do two potential stars have to be to count as the same

function center_of_mass(x::Int, y::Int)
  offset::Int32 = floor(EDGE / 2)
  cells = img[y - offset:y + offset, x - offset:x + offset]
  weights::Float32 = 0.0
  xsum::Float32 = 0.0
  ysum::Float32 = 0.0
  for ix = -offset:offset
    for iy = -offset:offset
      pixel = cells[ix + offset + 1, iy + offset + 1]
      pixval = min(pixel.r, pixel.g, pixel.b)
      weights += pixval
      xsum += pixval * ix
      ysum += pixval * iy
    end
  end
  centerx = xsum / weights
  centery = ysum / weights
  rx = x + centerx
  ry = y + centery
  return (rx, ty)
end

filen = "test.png"

print("Loading $filen... ")
img = load(filen)
if img == nothing
  println("Image failed to load.")
  exit(-1)
end

println("Done!")
sz = size(img)
px::Int32 = sz[1]
py::Int32 = sz[2]
println("Image size: ", px, "x", py, " pixels")
imagetype::Type = typeof(img[1, 1].r)
println("Image kind ", imagetype)
colcells::Int32 = floor((px - 2 * EDGE) / EDGE) # leave a margin of at least EDGE on all sides
rowcells::Int32 = floor((py - 2 * EDGE) / EDGE) # this allows the center_of_mass routing room to move

maxval = zeros(imagetype, rowcells, colcells)
maxx = zeros(UInt32, rowcells, colcells)
maxy = zeros(UInt32, rowcells, colcells)

# find max in each cell
yc::Int32 = EDGE
yptr::Int32 = 1
ptr::Int32 = 1
for bigy = 1:rowcells
  yp = EDGE * bigy
  for littley = 1:EDGE
    for bigx = 1:colcells
      xp = EDGE * bigx
      for littlex = 1:EDGE
        pixel = img[xp, yp]
        pixval = min(pixel.r, pixel.g, pixel.b)
        if pixval < CUTMAX & pixval > CUTMIN & pixval > maxval[bigx, bigy]
            maxval[bigx, bigy] = pixval
            maxx[bigx, bigy] = x
            maxy[bigx, bigy] = y
        end
        xp += 1
      end
    end
    yp += 1
  end
end

# Start processing with the brightest objects

order = sortperm(maxval, rev = true)
stars = Array{star}[]

while size(stars) < STARCOUNT && size(order) > 0
  starval = maxval[order[1]]
  starx = maxx[order[1]]
  stary = maxy[order[1]]
  onestar::star = refine(starx, stary)
  add = true
  if onestar.maxval > 0 # refine() returns a zero in maxval if convergence failed or shape test failed
    for i = 1:size(stars)
      if abs(starx - stars[i].x) <= EPSILON ||
         abs(stary - stars[i].y) <= EPSILON
         add = false
      end
    end
    if add
      push!(stars, onestar)
    end
  end
  shift!(order)
end

if size(stars) < 5
  println("Warning: Fewer than 5 stars found")
end

max_stars_repo_name: FatherGeoffHorton/AstroTools

max_issues_repo_name: FatherGeoffHorton/AstroTools

max_forks_repo_name: FatherGeoffHorton/AstroTools


text:
@info "Test decorator pattern..."

widget = VerticalScrollBar(Border(Window(2.5), 5.0))
draw(widget)


text:
using FeatureSelect
using Test

@testset "FeatureSelect.jl" begin
    # Write your tests here.
end


text:
@testset "Basic Utilities Tests" begin
  @testset "Utils" begin
    ϵ = log(3); δ = 0.05

    @test gaussianMechConstant(ϵ, δ) ≈ 1.7563398731147597 atol=1e-10
    tmp = gaussianMechConstant2(ϵ, δ)
    @test tmp[1] ≈ 1.7563398731147597 atol=1e-10
    @test tmp[2] ≈ 2.3095249420840473 atol=1e-10
  end

  @testset "Basic mechanisms" begin
    ϵ = log(2); δ = 0.04

    Random.seed!(12345)  # reseed the RNG
    D = Random.rand(100)
    @test gaussianMech(D, Distributions.mean, 2.0/100, ϵ, δ)[1] == 0.4764940748202759
    @test laplaceMech(D, Distributions.mean, 2.0/100, ϵ)[1] == 0.5366973705020675
    @test mean_dp(D, ϵ, δ) == 0.47514367413387
    @test mean_dp(D, ϵ, 0.0) == 0.4920613726730682
    @test mean_dp(D, ϵ) == 0.4886977015049402

    Random.seed!(12345)  # reseed the RNG
    f(x) = [x[1]^2*x[2], sqrt(x[1]), x[2]^3]
    x0 = [3, 4]
    l1sens = 0.5
    r,a = truncatedLaplaceMech(x0, f, l1sens, ϵ, δ)
    @test r[1] == 36.090422106340114
    @test r[2] == 2.50473361660456
    @test r[3] == 63.800895516558114
    @test a == 2.020992083320228
  end
end


text:
@testset "CUDA: Integration time" begin
    gpsl1 = GPSL1(use_gpu = Val(true))
    galileo_e1b = GalileoE1B(use_gpu = Val(true))
    max_integration_time = 2ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found)
    @test time == 1ms

    max_integration_time = 2ms
    bit_found = true
    time = @inferred Tracking.get_integration_time(gpsl1, max_integration_time, bit_found)
    @test time == 2ms

    max_integration_time = 2ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(
        galileo_e1b,
        max_integration_time,
        bit_found
    )
    @test time == 2ms

    max_integration_time = 10ms
    bit_found = false
    time = @inferred Tracking.get_integration_time(
        galileo_e1b,
        max_integration_time,
        bit_found
    )
    @test time == 4ms
end

@testset "CUDA: Number of chips to integrate GPU" begin
    gpsl1 = GPSL1(use_gpu = Val(true))
    max_integration_time = 1ms
    code_phase = 200
    bit_found = true
    chips = @inferred Tracking.get_num_chips_to_integrate(
        gpsl1,
        max_integration_time,
        code_phase,
        bit_found
    )
    @test chips == 1023 - 200

    code_phase = 1200
    chips = @inferred Tracking.get_num_chips_to_integrate(
        gpsl1,
        max_integration_time,
        code_phase,
        bit_found
    )
    @test chips == 2046 - 1200
end

@testset "CUDA: Number of samples to integrate" begin
    gpsl1 = GPSL1(use_gpu = Val(true))
    max_integration_time = 1ms
    sampling_frequency = 4e6Hz
    current_code_doppler = 0Hz
    current_code_phase = 0
    bit_found = true
    samples = @inferred Tracking.get_num_samples_left_to_integrate(
        gpsl1,
        max_integration_time,
        sampling_frequency,
        current_code_doppler,
        current_code_phase,
        bit_found
    )
    @test samples == 4000
end

@testset "CUDA: Code phase delta" begin
    gpsl1 = GPSL1(use_gpu = Val(true))
    code_doppler = 1Hz
    frequency = @inferred Tracking.get_current_code_frequency(gpsl1, code_doppler)
    @test frequency == 1023e3Hz + 1Hz
end

@testset "CUDA: Doppler aiding" begin
    gpsl1 = GPSL1(use_gpu = Val(true))
    init_carrier_doppler = 10Hz
    init_code_doppler = 1Hz
    carrier_freq_update = 2Hz
    code_freq_update = -0.5Hz
    velocity_aiding = 3Hz

    carrier_freq, code_freq = @inferred Tracking.aid_dopplers(
        gpsl1,
        init_carrier_doppler,
        init_code_doppler,
        carrier_freq_update,
        code_freq_update,
        velocity_aiding
    )

    @test carrier_freq == 10Hz + 2Hz + 3Hz
    @test code_freq == 1Hz + (2Hz + 3Hz) / 1540 - 0.5Hz
end

@testset "CUDA: Tracking with signal of type StructArray of CuArrays ComplexF32" begin
    CUDA.allowscalar() do
        gpsl1 = GPSL1(use_gpu = Val(true))
        carrier_doppler = 200Hz
        start_code_phase = 100
        code_frequency = carrier_doppler / 1540 + 1023e3Hz
        sampling_frequency = 4e6Hz
        prn = 1
        range = 0:3999
        start_carrier_phase = π / 2
        state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase, num_samples = 4000)

        carrier_phases = cu(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase)
        code_phases = get_code_frequency(gpsl1) / sampling_frequency .* range .+ start_code_phase
        upsampled_code = gpsl1.codes[1 .+ mod.(floor.(Int, code_phases), 1023), prn]
        signal = StructArray{ComplexF32}((cos.(carrier_phases), sin.(carrier_phases)))
        @. signal *= upsampled_code

        track_result = @inferred track(signal, state, sampling_frequency)

        iterations = 2000
        code_phases = zeros(iterations)
        carrier_phases = zeros(iterations)
        tracked_code_phases = zeros(iterations)
        tracked_carrier_phases = zeros(iterations)
        tracked_code_dopplers = zeros(iterations)
        tracked_carrier_dopplers = zeros(iterations)
        tracked_prompts = zeros(ComplexF64, iterations)
        for i = 1:iterations
            carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency +
                start_carrier_phase + π) - π
            code_phase = mod(
                code_frequency * 4000 * i / sampling_frequency + start_code_phase,
                1023
            )
            signal_temp = CuArray{ComplexF32}(cis.(
                    2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase
                ) .*
                get_code.(
                    gpsl1,
                    code_frequency .* range ./ sampling_frequency .+ code_phase,
                    prn
                )
            )
            signal = StructArray{ComplexF32}((real.(signal_temp), imag.(signal_temp)))
            track_result = @inferred track(
                signal,
                get_state(track_result),
                sampling_frequency
            )
            comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) /
                sampling_frequency + start_carrier_phase + π) - π
            comp_code_phase = mod(
                code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase,
                1023
            )
            tracked_code_phases[i] = get_code_phase(track_result)
            tracked_carrier_phases[i] = get_carrier_phase(track_result)
            tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz
            tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz
            tracked_prompts[i] = get_prompt(track_result)
            code_phases[i] = comp_code_phase
            carrier_phases[i] = comp_carrier_phase
        end

    @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4
    @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2
    end
#    using PyPlot
#    pygui(true)
#    figure("carrier_phases")
#    plot(tracked_carrier_phases)
#    plot(carrier_phases)
#    grid(true)
#    figure("Code Phases")
#    plot(300 * (tracked_code_phases .- code_phases))
#    figure("Carrier Doppler")
#    plot(tracked_carrier_dopplers)
#    figure("Code Doppler")
#    plot(tracked_code_dopplers)
#    figure("Prompt")
#    plot(real.(tracked_prompts))
#    plot(imag.(tracked_prompts))
end

@testset "CUDA: Track multiple signals with signal of type StructArray of CuArrays ComplexF32" begin
    CUDA.allowscalar() do
        gpsl1 = GPSL1(use_gpu = Val(true))
        carrier_doppler = 200Hz
        start_code_phase = 100
        code_frequency = carrier_doppler / 1540 + 1023e3Hz
        sampling_frequency = 4e6Hz
        prn = 1
        range = 0:3999
        start_carrier_phase = π / 2
        state = TrackingState(prn, gpsl1, carrier_doppler - 20Hz, start_code_phase, num_ants = NumAnts(3), num_samples = 4000)

        carrier_phases = cu(2π .* carrier_doppler .* range ./ sampling_frequency .+ start_carrier_phase)
        code_phases = get_code_frequency(gpsl1) / sampling_frequency .* range .+ start_code_phase
        upsampled_code = gpsl1.codes[1 .+ mod.(floor.(Int, code_phases), 1023), prn]
        carrier_phases_mat = carrier_phases * CUDA.ones(3)'
        signal_mat = StructArray{ComplexF32}((cos.(carrier_phases_mat), sin.(carrier_phases_mat)))
        @. signal_mat *= upsampled_code

        track_result = @inferred track(signal_mat, state, sampling_frequency)

        iterations = 2000
        code_phases = zeros(iterations)
        carrier_phases = zeros(iterations)
        tracked_code_phases = zeros(iterations)
        tracked_carrier_phases = zeros(iterations)
        tracked_code_dopplers = zeros(iterations)
        tracked_carrier_dopplers = zeros(iterations)
        for i = 1:iterations
            carrier_phase = mod2pi(2π * carrier_doppler * 4000 * i / sampling_frequency +
                start_carrier_phase + π) - π
            code_phase = mod(
                code_frequency * 4000 * i / sampling_frequency + start_code_phase,
                1023
            )
            signal_temp = CuArray{ComplexF32}(cis.(
                    2π .* carrier_doppler .* range ./ sampling_frequency .+ carrier_phase
                ) .*
                get_code.(
                    gpsl1,
                    code_frequency .* range ./ sampling_frequency .+ code_phase,
                    prn
                )
            )
            signal_temp_mat = repeat(signal_temp, outer = (1,3))
            signal_mat = StructArray(signal_temp_mat)
            track_result = @inferred track(
                signal_mat,
                get_state(track_result),
                sampling_frequency
            )
            comp_carrier_phase = mod2pi(2π * carrier_doppler * 4000 * (i + 1) /
                sampling_frequency + start_carrier_phase + π) - π
            comp_code_phase = mod(
                code_frequency * 4000 * (i + 1) / sampling_frequency + start_code_phase,
                1023
            )
            tracked_code_phases[i] = get_code_phase(track_result)
            tracked_carrier_phases[i] = get_carrier_phase(track_result)
            tracked_carrier_dopplers[i] = get_carrier_doppler(track_result)/Hz
            tracked_code_dopplers[i] = get_code_doppler(track_result)/Hz
            code_phases[i] = comp_code_phase
            carrier_phases[i] = comp_carrier_phase
        end

    @test tracked_code_phases[end] ≈ code_phases[end] atol = 2e-4
    @test tracked_carrier_phases[end] + π ≈ carrier_phases[end] atol = 5e-2
    end
#    using PyPlot
#    pygui(true)
#    figure("carrier_phases")
#    plot(tracked_carrier_phases)
#    plot(carrier_phases)
#    grid(true)
#    figure("Code Phases")
#    plot(300 * (tracked_code_phases .- code_phases))
#    figure("Carrier Doppler")
#    plot(tracked_carrier_dopplers)
#    figure("Code Doppler")
#    plot(tracked_code_dopplers)
end


text:
module PlutoWebcam

include("./notebook.jl")

end


text:
module CUFFT

using CUDAapi

using ..CuArrays
import ..CuArrays: unsafe_free!

using CUDAdrv
using CUDAdrv: CUstream

using CEnum

const libcufft = if Sys.iswindows()
    # no ccall by soname, we need the filename
    # NOTE: we discover the full path here, while only the wordsize and toolkit versions
    #       would have been enough to construct "cufft64_10.dll"
    toolkit = find_toolkit()
    path = find_cuda_library("cufft", toolkit)
    if path === nothing
        error("Could not find libcufft")
    end
    basename(path)
else
    # ccall by soname; CuArrays.__init__ will have populated Libdl.DL_LOAD_PATH
    "libcufft"
end

# core library
include("libcufft_common.jl")
include("error.jl")
include("libcufft.jl")

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

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

end


text:
"""
    capacitor(freq, C)
Create a `CircuitData` for a capacitor of capacitance `C` at frequencies
specified by `freq`. It uses the formula Z(ω) = 1/(jωC) where ω=2πf.
"""
capacitor(freq::Real, C::Number) = Impedance(1/(im * (2π * freq) * C))
capacitor(freq::AbstractVector{T}, C::Number) where {T<:Real} =
    CircuitData(freq, [capacitor(f, C) for f in freq])

"""
    inductor(freq, C)
Create a `CircuitData` for an inductor of inductance `L` at frequencies
specified by `freq`. It uses the formula Z(ω) = jωL where ω=2πf.
"""
inductor(freq::Real, L::Number) = Impedance(im * (2π * freq) * L)
inductor(freq::AbstractVector{T}, L::Number) where {T<:Real} =
    CircuitData(freq, [inductor(f, L) for f in freq])

"""
    resistor(freq, R)
Create a `CircuitData` for a resistor of resistance `R` at frequencies
specified by `freq`. It uses the formula Z(ω) = R where ω=2πf.
"""
resistor(freq::Real, R::Number) = Impedance(R)
resistor(freq::AbstractVector{T}, R::Number) where {T<:Real} =
    CircuitData(freq, [Impedance(R) for f in freq])


text:
## import Helmholtz-solver library
using WavesAndEigenvalues.Helmholtz
## load mesh from file
mesh=Mesh("../NTNU/NTNU.msh",scale=1.0)
#create unit and full mesh from half-cell
# Important: full merges all the domains that are copied from the unit cell
# unit: keeps the domains separate and appends numbers!
doms=[("Interior",:full),("Inlet",:full), ("Outlet_high",:full), ("Outlet_low",:full), ("Flame",:unit),]#
collect_lines!(mesh)
unit_mesh=extend_mesh(mesh,doms,unit=true)
full_mesh=extend_mesh(mesh,doms,unit=false)

## Helmholtz solver
#define speed of sound field
function speedofsound(x,y,z)
    if z<0.415
        return 347.0#m/s
    else
        return 850.0#m/s
    end
end

c=generate_field(unit_mesh,speedofsound,0)
C=generate_field(full_mesh,speedofsound,0)
#D=Dict()

#vtk_write("NTNU_c",full_mesh,D)
##describe model
D=Dict()
D["Interior"]=(:interior,())
D["Outlet_high"]=(:admittance, (:Y_in,1E15))
D["Outlet_low"]=(:admittance, (:Y_out,1E15))
##discretize models
L=discretize(full_mesh, D, C, order=:1,)
l=discretize(unit_mesh, D, c, order=:1, b=:b)
## solve model using beyn
#(type Γ by typing \Gamma and Enter)
Γ=[150.0+5.0im, 150.0-5.0im, 1000.0-5.0im, 1000.0+5.0im].*2*pi #corner points for the contour (in this case a rectangle)
Ω, P = beyn(L,Γ,l=10,N=64, output=true)
ω, p = beyn(l,Γ,l=10,N=128, output=true)
##verify beyn's solution by local householder iteration
tol=1E-10
D=Dict()
for idx=1:length(Ω)
    sol,nn,flag=householder(L,Ω[idx],maxiter=6,order=5,tol=tol,n_eig_val=3,output=false);
    mode="[$(string(round((Ω[idx]/2/pi),digits=2)))]Hz"
    if abs(sol.params[:ω]-Ω[idx])<5*tol
        println("Kept:$mode")
        D["abs:"*mode]=abs.(sol.v)./maximum(abs.(sol.v)) #note that in Julia any function can be performed elementwise by putting a . before the paranthesis
        D["phase:"*mode]=angle.(sol.v)
    else
        convmode="[$(string(round((sol.params[:ω]/2/pi),digits=2)))]Hz"
        println("Excluded:$mode because converged to $convmode")
    end
end
##
D=Dict()
D["speedofsound"]=C
vtk_write("NTNU",full_mesh,D)
##

sol,nn,flag=householder(l,1000*2*pi,maxiter=6,order=1,tol=1E-9,n_eig_val=1,output=true);
mode="[$(string(round((sol.params[:ω]/2/pi),digits=2)))]Hz"
v=bloch_expand(unit_mesh,sol)
D["abs Bloch:"*mode]=abs.(v)./maximum(abs.(v)) #note that in Julia any function can be performed elementwise by putting a . before the paranthesis
D["phase Bloch:"*mode]=angle.(v)
#ol,nn,flag=householder(L,sol.params[:ω],maxiter=6,order=5,tol=1E-9,n_eig_val=3,output=true);
Sol,nn,flag=householder(L,sol.params[:ω],maxiter=6,order=5,tol=1E-5,n_eig_val=3,output=true);

mode="[$(string(round((Sol.params[:ω]/2/pi),digits=2)))]Hz"
D["abs:"*mode]=abs.(Sol.v)./maximum(abs.(Sol.v)) #note that in Julia any function can be performed elementwise by putting a . before the paranthesis
D["phase:"*mode]=angle.(Sol.v)
vtk_write("NTNU",full_mesh,D)

##

##

D=Dict()
D["speedofsound"]=C
Sol,nn,flag=householder(L,500*2*pi,maxiter=6,order=5,tol=1E-5,n_eig_val=3,output=true);
mode="[$(string(round((Sol.params[:ω]/2/pi),digits=2)))]Hz"
D["abs:"*mode]=abs.(Sol.v)./maximum(abs.(Sol.v)) #note that in Julia any function can be performed elementwise by putting a . before the paranthesis
D["phase:"*mode]=angle.(Sol.v)
vtk_write("NTNU",full_mesh,D)

##


text:
module TestOptional
using StaticArrays
using Rotations
using CoordinateTransformations
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using Accessors
using PhaseSpaceIO
using PhaseSpaceIO: propagate_z, position
using PhaseSpaceIO.Testing
using Accessors
using Test
using LinearAlgebra

@testset "SetZ" begin
    p = arbitrary(EGSParticle{Float32})
    z = randn(Float32)
    p2 = SetZ(z)(p)
    @test p2 isa EGSParticleZ
    @test p2.z === z
    @test p2 === @set p.z = z
end

@testset "propagate_z" begin
    for p in [
        arbitrary(EGSParticleZ{Float32}),
        arbitrary(EGSParticleZ{Nothing}),
        arbitrary(IAEAParticle{1,1}),
        arbitrary(IAEAParticle{0,1}),
        arbitrary(IAEAParticle{0,0}),
        (
         position=@SVector(randn(3)),
         direction = normalize(@SVector(randn(3))),
        ),
        (
         position=@SVector(randn(3)),
         direction = normalize(@SVector(randn(3))),
         E = rand(),
         weight = rand(),
        )
        ]

        p = @set direction(p) = @SVector[0,0,1]
        z_to = randn()
        p2 = @inferred propagate_z(p, z_to)
        p2_expected =@set position(p)[3] = position(p2)[3]
        @test position(p2)[3] ≈ z_to
        @test p2 === p2_expected
    end

    dir = normalize([0,1,1])
    pos = [0,0,-100]
    pos_new = [0,100,0]
    p = arbitrary(IAEAParticle{0,0})
    p = @set position(p) = pos
    p = @set direction(p) = dir
    p2 = propagate_z(p, 0.0)
    @test p2 == (@set position(p) = pos_new)

    dir = normalize([0,1,1])
    pos = randn(Float32, 3)
    z_to = pos[3] + 100
    p = arbitrary(EGSParticleZ{Nothing})
    p = @set position(p) = pos
    p = @set direction(p) = dir
    p2 = propagate_z(p, z_to)
    x,y,z=pos
    pos_new = [x, y+100, z+100]
    @test position(p2) == pos_new

    p = arbitrary(EGSParticle{Nothing})
    @test_throws ArgumentError propagate_z(p, 1f0)
end

@testset "position, direction" begin
    dir = StaticArrays.normalize(@SVector(randn(3)))
    pos_iaea = @SVector(randn(3))
    pos_egs = @SVector(randn(2))
    for (p, pos, dir) in [
                          (arbitrary(EGSParticle{Float32}), pos_egs, dir),
                          (arbitrary(IAEAParticle{2,3}), pos_iaea, dir),
           ]
        @inferred position(p)
        @inferred direction(p)
        q = @set position(p) = pos
        @test typeof(q) == typeof(p)
        @test q.E === p.E
        @test position(q) ≈ pos
        q = @set direction(p) = dir
        @test typeof(q) == typeof(p)
        @test q.E === p.E

        @test direction(q) ≈ dir
        @test dir ≈ [q.u, q.v, q.w]
    end

    p = arbitrary(EGSParticle{Nothing})
    @test position(p)   === @SVector[p.x, p.y]
    @test position(p,3) === @SVector[p.x, p.y, 3]
    p = arbitrary(EGSParticleZ{Float32})
    @test position(p)   === @SVector[p.x, p.y, p.z]
    p = arbitrary(IAEAParticle{1,2})
    @test position(p)   === @SVector[p.x, p.y, p.z]
end

@testset "coordinate transformations" begin
    p = arbitrary(IAEAParticle{1,2})

    t = Translation(randn(3))
    @test position(t(p)) ≈ t(position(p))
    @test direction(p) ≈ direction(t(p))

    l = LinearMap(rand(AngleAxis))
    @test position(l(p)) ≈ l(position(p))
    @test l(direction(p)) ≈ direction(l(p))

    tl = t ∘ l
    @test position(tl(p)) ≈ position(t(l(p)))
    @test direction(tl(p)) ≈ direction(l(p))
end

end#module


text:
using BenchmarkTools
using JultraDark
using JultraDark: psi_half_step!, psi_whole_step!, phi_whole_step!, take_steps!
using MPI

if ~MPI.Initialized()
    MPI.Init()
end

n = 64
Δt = 0.1

suite = BenchmarkGroup()

for (label, grids) in [("Grids", Grids(1.0, n)), ("PencilGrids", PencilGrids(1.0, n))]

    suite[label] = BenchmarkGroup()

    suite[label]["psi_half_step"] = @benchmarkable psi_half_step!($Δt, $grids)
    suite[label]["psi_whole_step"] = @benchmarkable psi_whole_step!($Δt, $grids)
    suite[label]["phi_whole_step"] = @benchmarkable phi_whole_step!($Δt, $grids)

    suite[label]["take_steps"] = @benchmarkable take_steps!($grids, $1.0, $Δt, 1, $(OutputConfig("output", 1:2)), $(t -> 1))
    suite[label]["take_steps 10"] = @benchmarkable take_steps!($grids, $1.0, $Δt, 10, $(OutputConfig("output", 1:2)), $(t -> 1))

end

@show tune!(suite)

results = run(suite)

@show res_min = minimum(results)
@show res_med = median(results)

judge(res_med["Grids"], res_med["PencilGrids"])


text:
using MacroTools
using Poly

export compile

"""
SYMBOLS FUNCTIONS
"""

"""
get all symbols in a symbol
"""
get_all_symbols(s::Symbol)::Set{Symbol} = Set([s])

"""
get all symbols in a number
"""
get_all_symbols(n::Number)::Set{Symbol} = Set{Symbol}()

"""
get all symbols in a LineNumberNode
"""
get_all_symbols(n::LineNumberNode)::Set{Symbol} = Set{Symbol}()

"""
get all symbols in Any
"""
get_all_symbols(::Any)::Set{Symbol} = Set{Symbol}()

"""
get all symbols in an expression
"""
function get_all_symbols(expr::Expr)::Set{Symbol}
    symbols = Set{Symbol}()
    if expr.head == :call
        # if function call, include all function args
        for arg in expr.args[2:end]
            union!(symbols, get_all_symbols(arg))
        end
    elseif expr.head == :macrocall
        for arg in expr.args[2:end]
            union!(symbols, get_all_symbols(arg))
        end
    else
        for arg in expr.args
            union!(symbols, get_all_symbols(arg))
        end
    end
    return symbols
end

"""
get all symbols in a kernel
"""
function get_all_symbols(kernel::LoopKernel)::Set{Symbol}
    symbols = Set{Symbol}()

    for instruction in kernel.instructions
        union!(symbols, get_all_symbols(instruction.body))
    end

    for domain in kernel.domains
        union!(symbols, get_all_symbols(domain.step))
        union!(symbols, get_all_symbols(domain.lowerbound))
        union!(symbols, get_all_symbols(domain.upperbound))
        push!(symbols, domain.iname)
    end

    return symbols
end

"""
KERNEL ARG FUNCTIONS
"""

"""
get all function arguments to kernel
"""
function get_kernel_args(kernel::LoopKernel; use_defined=false)::Set{Symbol}
    all_symbols = get_all_symbols(kernel)
    defined_symbols = Set{Symbol}()

    # get all symbols defined in LHS of instructions
    for instruction in kernel.instructions
        lhs = instruction.body.args[1]
        if typeof(lhs) == Symbol && instruction.body.head == :(=) && !use_defined
            push!(defined_symbols, lhs)
        end
    end

    # get all loop domain symbols
    for domain in kernel.domains
        push!(defined_symbols, domain.iname)
    end

    # args are all symbols that are not defined symbols
    args = Set([s for s in all_symbols if !(s in defined_symbols)])

    append!(kernel.args, args)
    return args
end

"""
CONSTANT FINDING FUNCTIONS
"""

"""
determine all constants in kernel
"""
function set_kernel_consts(kernel::LoopKernel)
    const_symbols = Set{Symbol}()
    domain_inames = [domain.iname for domain in kernel.domains]

    # get all symbols in loop bounds (since must be constant)
    for domain in kernel.domains
        union!(const_symbols, get_kernel_consts(domain.lowerbound, domain_inames))
        union!(const_symbols, get_kernel_consts(domain.upperbound, domain_inames))
        union!(const_symbols, get_kernel_consts(domain.step, domain_inames))
    end

    append!(kernel.consts, const_symbols)
end

function get_kernel_consts(ex::Expr, domain_inames::Vector{Symbol})::Set{Symbol}
    const_symbols = Set{Symbol}()
    for arg in ex.args[2:end]
        union!(const_symbols, get_kernel_consts(arg, domain_inames))
    end
    return const_symbols
end

function get_kernel_consts(sym::Symbol, domain_inames::Vector{Symbol})::Set{Symbol}
    if !(sym in domain_inames)
        return Set([sym])
    end
    return Set()
end

function get_kernel_consts(any, domain_inames::Vector{Symbol})::Set{Symbol}
    return Set()
end

"""
COMPILATION
"""

"""
compile native julia code given a kernel
"""
function compile(kernel::LoopKernel; debug=false, verbose=0, tile=-1, thread=false)
    # kernel args
    args = get_kernel_args(kernel)
    set_kernel_consts(kernel)

    body = run_polyhedral_model(kernel, debug=debug, verbose=verbose, tile=tile, thread=thread)

    expr = quote
        function $(gensym(:Poly))(;$(args...))
            $(body)
        end
    end
    eval(expr)
end

"""
compile native julia code given a kernel to an expression
"""
function compile_expr(kernel::LoopKernel; debug=false, verbose=0, tile=-1, thread=false)::Expr
    # kernel "args" for isl
    get_kernel_args(kernel)
    set_kernel_consts(kernel)

    expr = run_polyhedral_model(kernel, debug=debug, verbose=verbose, tile=tile, thread=thread)

    return expr
end


text:
struct Fourier{T} <: Basis{T} end
Fourier() = Fourier{Float64}()

==(::Fourier, ::Fourier) = true

axes(F::Fourier) = (Inclusion(ℝ), _BlockedUnitRange(1:2:∞))

function getindex(F::Fourier{T}, x::Real, j::Int)::T where T
    isodd(j) && return cos((j÷2)*x)
    sin((j÷2)*x)
end

### transform
checkpoints(::Fourier{T}) where T = T[1.223972,3.14,5.83273484]

fouriergrid(T, n) = convert(T,π)*collect(0:2:2n-2)/n

function grid(Pn::SubQuasiArray{T,2,<:Fourier,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
    kr,jr = parentindices(Pn)
    n = maximum(jr)
    fouriergrid(T, n)
end

"""
Gives a shuffled version of the real FFT, with order
1,sin(θ),cos(θ),sin(2θ)…
"""
struct ShuffledRFFT{T,Plan} <: Factorization{T}
    plan::Plan
end

size(F::ShuffledRFFT, _) = size(F.plan,1)
size(F::ShuffledRFFT) = (size(F.plan,1),size(F.plan,1))

ShuffledRFFT{T}(p::Plan) where {T,Plan} = ShuffledRFFT{T,Plan}(p)
ShuffledRFFT{T}(n::Int) where T = ShuffledRFFT{T}(FFTW.plan_r2r!(Array{T}(undef, n), FFTW.R2HC))

function *(F::ShuffledRFFT{T}, b::AbstractVector) where T
    n = size(F,1)
    c = lmul!(convert(T,2)/n, F.plan * convert(Array, b))
    c[1] /= 2
    iseven(n) && (c[n÷2+1] /= 2)
    negateeven!(reverseeven!(interlace!(c,1)))
end

factorize(L::SubQuasiArray{T,2,<:Fourier,<:Tuple{<:Inclusion,<:OneTo}}) where T =
    TransformFactorization(grid(L), ShuffledRFFT{T}(size(L,2)))

import BlockBandedMatrices: _BlockSkylineMatrix

@simplify function *(A::QuasiAdjoint{<:Any,<:Fourier}, B::Fourier)
    TV = promote_type(eltype(A),eltype(B))
    PseudoBlockArray(Diagonal(Vcat(2convert(TV,π),Fill(convert(TV,π),∞))), (axes(A,1),axes(B,2)))
end

@simplify function *(D::Derivative, F::Fourier)
    TV = promote_type(eltype(D),eltype(F))
    Fourier{TV}()*_BlockArray(Diagonal(Vcat([reshape([0.0],1,1)], (1.0:∞) .* Fill([0 -one(TV); one(TV) 0], ∞))), (axes(F,2),axes(F,2)))
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), c::BroadcastQuasiVector{<:Any,typeof(cos),<:Tuple{<:Inclusion{<:Any,<:FullSpace}}}, F::Fourier)
    axes(c,1) == axes(F,1) || throw(DimensionMismatch())
    T = promote_type(eltype(c), eltype(F))
    F*mortar(Tridiagonal(Vcat([reshape([0; one(T)],2,1)], Fill(Matrix(one(T)/2*I,2,2),∞)),
                        Vcat([zeros(T,1,1)], Fill(Matrix(zero(T)I,2,2),∞)),
                        Vcat([[0 one(T)/2]], Fill(Matrix(one(T)/2*I,2,2),∞))))
end

function broadcasted(::LazyQuasiArrayStyle{2}, ::typeof(*), s::BroadcastQuasiVector{<:Any,typeof(sin),<:Tuple{<:Inclusion{<:Any,<:FullSpace}}}, F::Fourier)
    axes(s,1) == axes(F,1) || throw(DimensionMismatch())
    T = promote_type(eltype(s), eltype(F))
    F*mortar(Tridiagonal(Vcat([reshape([one(T); 0],2,1)], Fill([0 one(T)/2; -one(T)/2 0],∞)),
                        Vcat([zeros(T,1,1)], Fill(Matrix(zero(T)*I,2,2),∞)),
                        Vcat([[one(T)/2 0]], Fill([0 -one(T)/2; one(T)/2 0],∞))))
end


text:
using Plots

include("helper_functions.jl")

begin
    λ1 = 1.0
    λ2 = 1.0
    α = 1.0
    μ(α) = 1 / (1.5 * (α))
    ν(α) = 1 / (0.7 * (α))
    qs = 0:0.01:1
    parameters = (λ1, λ2, α, μ, ν)
end

# function pstar(q, λ1, λ2, α, μ, ν)
#     p = (λ1*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) - μ(α)*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) - ν(α)*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) - λ1*ν(α) + 2*μ(α)*ν(α) + λ2*q*Complex(μ(α)*(ν(α) - λ2*q))^(1/2) + λ1*λ2*q - 2*λ2*μ(α)*q)/(λ1*μ(α) - λ1*ν(α) + λ1*λ2*q) |> real
#     return clamp(p, 0, 1)
# end

# function best_responded_cost(p, λ1, λ2, α, μ, ν)
#     cost(p, pstar(p, λ2, λ1, α, μ, ν), parameters...)
# end

# function ∂best_responded_cost_∂p(p, λ1, λ2, α, μ, ν)
#     ForwardDiff.derivative(ϕ -> best_responded_cost(ϕ, λ1, λ2, α, μ, ν), p)
# end

# using BenchmarkTools
# brc = @benchmark best_responded_cost(0.5, λ1, λ2, α, μ, ν)
# abrc = @benchmark analytic_best_responded_cost(0.5, λ1, λ2, α, μ, ν)

# function analytic_best_responded_cost(p, λ1, λ2, α, μ, ν)
#     real(- (p - 1)/(μ(α) - λ1 + λ1*p) - p/(λ1*p - ν(α) + (λ2*(λ2*Complex(μ(α)*(ν(α) - λ1*p))^(1/2) - μ(α)*Complex(μ(α)*(ν(α) - λ1*p))^(1/2) - ν(α)*Complex(μ(α)*(ν(α) - λ1*p))^(1/2) - λ2*ν(α) + 2*μ(α)*ν(α) + λ1*p*Complex(μ(α)*(ν(α) - λ1*p))^(1/2) + λ1*λ2*p - 2*λ1*μ(α)*p))/(λ2*μ(α) - λ2*ν(α) + λ1*λ2*p)))
# end

# function diff_best_responded_cost(p, λ1, λ2, α, μ, ν)
#     squareroot_factor = real(Complex(μ(α)*(ν(α) - λ1*p))^(1/2))
#     A = (
#         (λ1*(p - 1)) / (μ(α) - λ1 + λ1*p)^2 -
#         (μ(α) - ν(α) + λ1*p)/
#         ((squareroot_factor - ν(α) + λ1*p)*(λ2 - μ(α) - ν(α) + λ1*p)) -
#         1/(μ(α) - λ1 + λ1*p)
#     )
#     B = (
#         p * λ1 *(
#             4*real(Complex(μ(α)*(ν(α) - λ1*p))^(3/2)) -
#             μ(α)^3 +
#             λ2*μ(α)^2 -
#             2*λ1^2*p^2*squareroot_factor -
#             4*μ(α)^2*(ν(α) - λ1*p) +
#             2*μ(α)^2*squareroot_factor -
#             2*ν(α)^2*squareroot_factor -
#             λ1*μ(α)*p*(ν(α) - λ1*p) +
#             4*λ1*ν(α)*p*squareroot_factor +
#             λ2*μ(α)*(ν(α) - λ1*p) +
#             μ(α)*ν(α)*(ν(α) - λ1*p) -
#             2*λ2*μ(α)*squareroot_factor
#         )
#     )

#     C = (
#         2*squareroot_factor *
#         (
#             squareroot_factor -
#             ν(α) +
#             λ1*p
#         )^2 *
#         (λ2 - μ(α) - ν(α) + λ1*p)^2
#     )
#     return real(A - B/C)
# end

# function equilibrium_probability(λ1, λ2, α, μ, ν; tolerance = 10e-9, first_actor = 1, returnall = false)
#     if first_actor == 1
#         parameters = (λ1, λ2, α, μ, ν)
#         swapped_parameters = (λ2, λ1, α, μ, ν)
#     else
#         parameters = (λ2, λ1, α, μ, ν)
#         swapped_parameters = (λ1, λ2, α, μ, ν)
#     end

#     root, ismin = classify_roots(p -> diff_best_responded_cost(p, parameters...), find_valid_roots(p -> diff_best_responded_cost(p, parameters...)))
#     if sum(ismin) == 0
#         val_at_zero = analytic_best_responded_cost(0,  parameters...)
#         val_at_one = analytic_best_responded_cost(0,  parameters...)
#         root = minimum((val_at_zero, val_at_one))
#     else
#         root = first(root[ismin])
#     end
#     first_pstar(p) = pstar(p, parameters...)
#     second_pstar(p) = pstar(p, swapped_parameters...)

#     diff = 100
#     response1 = 0
#     ps = Float64[]
#     qs = Float64[]
#     if tolerance < 10e-13
#         # println("tolerance of less than 10e-13 is dangerous")
#         tolerance = 10e-14
#     end
#     while diff > tolerance
#         response1 = second_pstar(root)
#         response2 = first_pstar(response1)
#         push!(ps, response2)
#         push!(qs, response2)
#         diff = abs(response2 - root)
#         root = response2
#     end
#     if returnall == true
#         return ps, qs
#     else
#         return root, response1
#     end
# end

# Static Graph
λ1s = 0.001:0.001:1.755

probability_tuples = equilibrium_probability.(λ1s, λ2, 1, μ, ν; tolerance = 10e-9)
response_probs = first.(probability_tuples)
response_qrobs = last.(probability_tuples)

plot(λ1s, response_probs, label = "p⋆", ylims = (0,1), xlims = (0,1.76))
plot!(λ1s, response_qrobs, label = "q⋆")
plot!(xlabel = "λ1", ylabel = "probability",
    title = "optimal probabilities for various λ1")

# Animation
λ2s = 0:0.01:2
animated_λs = @animate for λ2 in λ2s
    λ1s = 0:0.01:2
    probability_tuples = equilibrium_probability.(λ1s, λ2, 1, μ, ν; tolerance = 10e-9)
    response_probs = first.(probability_tuples)
    response_qrobs = last.(probability_tuples)
    plot(title = "Traffic at each server for varying λ1 for λ2 = $(dp_string(λ2, 3))",
        legend = :outerright,
        size = (800, 600))
    plot!(λ1s, response_probs, label = "p⋆", ylims = (0,1), xlims = (0,1.76))
    plot!(λ1s, response_qrobs, label = "q⋆")
    plot!(xlabel = "λ1", ylabel = "probability")

    p1_assumptions = calculate_assumptions(response_probs, response_qrobs, λ1s, λ2, α, μ, ν; which = "both")
    p2_assumptions = calculate_assumptions(response_probs, response_qrobs, λ2, λ1s, α, μ, ν; which = "both")
    μ1_shapes = shapes_from_intervals(shaded_regions(p1_assumptions[1], λ1s))
    ν1_shapes = shapes_from_intervals(shaded_regions(p1_assumptions[2], λ1s))
    μ2_shapes = shapes_from_intervals(shaded_regions(p2_assumptions[1], λ1s))
    ν2_shapes = shapes_from_intervals(shaded_regions(p2_assumptions[2], λ1s))
    vec_shapes = [μ1_shapes, ν1_shapes, μ2_shapes, ν2_shapes]
    resized_shapes = resize_interval_shapes(vec_shapes)

    plot!(resized_shapes[1], colour=:red, label = ["μ overload p1" fill(nothing, 1, length(vec_shapes)-1)], alpha=0.3)
    plot!(resized_shapes[2], colour=:blue, label = ["ν overload p1" fill(nothing, 1, length(vec_shapes)-1)], alpha=0.3)
    plot!(resized_shapes[3], colour=:green, label = ["μ overload p2" fill(nothing, 1, length(vec_shapes)-1)], alpha=0.3)
    plot!(resized_shapes[4], colour=:yellow, label = ["ν overload p2" fill(nothing, 1, length(vec_shapes)-1)], alpha=0.3)
end

gif(animated_λs, "altering_arrival_rates.gif", fps = 30)

plot(0.5:0.0001:0.8, p -> diff_best_responded_cost(p, 1.75, 1, 1, μ, ν), ylims = (-500, 500))
plot!(0.5:0.001:0.8, p -> analytic_best_responded_cost(p, 1.75, 1, 1, μ, ν))
hline!([0], label = "", colour = :black)

begin
    roots, ismin = classify_roots(p -> diff_best_responded_cost(p, 1.75, 1, 1, μ, ν), find_valid_roots(p -> diff_best_responded_cost(big(p), 1.75, 1, 1, μ, ν)))
    vline!(roots[ismin], colour = :black, label = "min")
    vline!(roots[.!ismin], colour = :red, label = "max")
end

roots[ismin]

difficult_params = (1.75, 1, 1, μ, ν)
difficult_p = equilibrium_probability(difficult_params...)[1]

plot(0:0.001:1, p -> best_responded_cost(p, difficult_params...), ylims = (-50, 50))
vline!([difficult_p])

difficult_br_costs = best_responded_cost.(0.5:0.001:1, difficult_params...)
min_index = findmin(difficult_br_costs)[2]
min_at = 0.5+0.001*318

plot(0.5:0.001:1, p -> best_responded_cost(p, difficult_params...))
vline!([min_at])


text:
include("../src/OteraEngine.jl")
using .OteraEngine
using Documenter

makedocs(;
    modules=[OteraEngine],
    authors="QGMW22 <ascendwatson@gmail.com> and contributors",
    sitename="OteraEngine.jl",
    pages=[
        "Home" => "index.md",
        "Tutorial" => "tutorial.md"
    ],
)


text:
"""
Unimodal Normal Distribution Crossover (UNDX).

See
    Kita, H., Ono, I., and Kobayashi, S., "Multi-parental Extension of the Unimodal Normal Distribution Crossover for Real-coded Genetic Algorithms," Proceedings of the 1999 Congress on Evolutionary Computation, pp. 1581-1588, 1999.
    Deb, K., Anand, A., and Joshi, D., "A Computationally Efficient Evolutionary Algorithm for Real-Parameter Optimization," Evolutionary Computation, vol. 10, no. 4, pp. 371-395, 2002.
"""
struct UnimodalNormalDistributionCrossover{NP} <: CrossoverOperator{NP,1}
	ζ::Float64 # sd for parent-defined orthogonal directions
    η::Float64 # sd for the directions orthogonal to the parent-defined ones, divided by √n

    function UnimodalNormalDistributionCrossover{NP}(ζ::Number, η::Number) where NP
        ζ > 0 || throw(ArgumentError("ζ must be positive"))
        η >= 0 || throw(ArgumentError("η must be non-negative"))
        new{NP}(ζ, η)
    end
    UnimodalNormalDistributionCrossover{NP}(params::Parameters) where NP =
        new{NP}(params[:UNDX_ζ], params[:UNDX_η])
end

const UNDX_DefaultOptions = ParamsDict(
    :UNDX_ζ => 0.5,
    :UNDX_η => 0.1
)

function apply!(xover::UnimodalNormalDistributionCrossover{NP},
                target::AbstractIndividual, targetIndex::Int,
                pop, parentIndices) where NP
    @assert length(parentIndices) == NP

    parents = pop[:, parentIndices]
    center = mean(parents, dims=2)
    parents .-= center
    sd = mean(map(sqrt, sum(abs2, parents, dims=2)))
    svdf = svd(parents, full=false)
    svals = svdf.S / norm(svdf.S)
    @inbounds for (i,s) in enumerate(svals)
        svals[i] = sd*randn() * (
            abs(s) > 1E-8 ? # check if parent subspace; FIXME is the threshold correct?
            s * xover.ζ : xover.η )
    end

    mul!(target, svdf.U, svals)
    @inbounds @simd for i in eachindex(target)
        target[i] += center[i,1]
    end
    return target
end


text:
using Base.Test
using LyapunovExponents
using LyapunovExponents: IntegratorError
using LyapunovExponents.TestTools: @test_nothrow

@testset "User interface" begin
    args = [
        (x...) -> nothing,  # f
        [0.0],              # u0
        nothing,            # p
    ]
    for lep in [DiscreteLEProblem, ContinuousLEProblem]
        @test_nothrow lep(args...; t_attr=1)
        try
            lep(args...)
            println("Calling $lep w/o t_attr didn't throw.")
            @test false
        catch err
            @test isa(err, ErrorException)
            @test contains(err.msg, "`t_attr` is required")
        end
    end
end

@testset "IntegratorError" begin
    prob = ContinuousLEProblem(
        (du, u, p, t) -> (@. du = exp(u)),  # f
        [0.0],    # u0 --- explosion time is 1.0
        nothing,  # p
        t_attr = 10.0,
        t_tran = 0,
    )
    err = try
        solve(prob, integrator_options=[:verbose => false])
        nothing
    catch err
        err
    end
    @test err isa IntegratorError
    @test err.retcode == :DtLessThanMin
end


text:
using DiffEqOperators, SparseArrays, LinearAlgebra

x = 0.0:0.01:π
x_ = x[2:end-1]
dx = diff(x)
y = sin.(2x)
y_ = y[2:(end-1)]
dy = [2cos.(2x_),-4sin.(2x_),-8cos.(2x_),16y_,32cos.(2x_)]

# test concretization against regular grid operator
for dor in 1:4, aor in 2:2:6
    Dr = CenteredDifference(dor,aor,dx[1],length(x)-2)
    Dir = CenteredDifference(dor,aor,dx,length(x)-2)

    @test sparse(Dr)≈sparse(Dir)
    @test Array(Dr)≈Array(Dir)

end

# test irregular grid operator with regular grid
for dor in 1:4, aor in 2:6

    D1 = CenteredDifference(dor,aor,dx,length(x)-2)

    #take derivative
    yprime1 = D1*y

    #test result
    @test yprime1 ≈ dy[dor] atol = 10.0^(1-aor)#2test operator with known derivative of exp(kx)

    #TODO: implement specific tests for the left and right boundary regions, waiting until after update
end

# test irregular grid

x = sin.(0.0:0.05:π)
x = cumsum(x)
x = x/x[end]*π
x_ = x[2:end-1]
dx = diff(x)
y  = sin.(2x)
y_ = y[2:(end-1)]
dy = [2cos.(2x_), -4sin.(2x_), -8cos.(2x_), 16y_, 32cos.(2x_)]

for dor in 1:4, aor in 4:10

    D1 = CenteredDifference(dor,aor,dx,length(x)-2)

    #take derivative
    yprime1 = D1*y

    #test result
    tol = 2*10.0^(2-aor)*maximum(dx)^(2-dor)
    # error estimate is fairly difficult for high order derivatives and small dx.

    # err = norm(yprime1.-dy[dor])
    # @show err
    # @show tol

    @test yprime1 ≈ dy[dor] atol = tol #2test operator with known derivative of exp(kx)

    #TODO: implement specific tests for the left and right boundary regions, waiting until after update
end


text:
# Compute Eccentric Anomally from Mean Anomally
# WARNING: Only valid for small e.  Not convergent for e>~0.66
# Therefore, this should only be used for testing purposes, not general science applications
function ecc_anom_bessel_series_approx{T}(mean_anom::T, ecc::T; tol::Real = 1.e-12)
  sinM = sin(mean_anom)
  sin2M = sin(2*mean_anom)
  sin3M = sin(3*mean_anom)
  sin4M = sin(4*mean_anom)
  E = mean_anom+ecc*(sinM+ecc*(0.5*sin2M+ecc*((3*sin3M-sinM)/8+ecc*((sin4M-0.5*sin2M)/3))))
end

#=  Compute Eccentric Anomally from Mean Anomally robustly
    Initial guess for eccentric anomaly for use by itterative solvers of Kepler's equation
    Based on "The Solution of Kepler's Equations - Part Three"
            Danby, J. M. A. (1987) Journal: Celestial Mechanics, Volume 40, Issue 3-4, pp. 303-312
            1987CeMec..40..303D
=#
function ecc_anom_init_guess_danby{T}(M::T, ecc::T)
    local k, E
    const k = convert(T,0.85)
    if(M<zero(T)) M += 2pi end
    E = (M<pi) ? M + k*ecc : M - k*ecc;
end

#= Update the current guess for solution to Kepler's equation
   Based on "An Improved Algorithm due to Laguerre for the Solution of Kepler's Equation"
   Conway, B. A.  (1986) Celestial Mechanics, Volume 39, Issue 2, pp.199-211
  1986CeMec..39..199C
=#
function update_ecc_anom_laguerre{T}(E::T, M::T, ecc::T)
  es = ecc*sin(E)
  ec = ecc*cos(E)
  F = (E-es)-M
  Fp = one(T)-ec
  Fpp = es
  const n = 5
  root = sqrt(abs((n-1)*((n-1)*Fp*Fp-n*F*Fpp)))
  denom = Fp>zero(T) ? Fp+root : Fp-root
  return E-n*F/denom
end

# Loop to update the current estimate of the solution to Kepler's equation
function calc_ecc_anom_itterative_laguerre{T}(mean_anom::T, ecc::T; tol::Real = 1.e-8)
    local M, E, E_old
    #M = atan2(sin(mean_anom),cos(mean_anom))  # WARNING: Kludge since missing autodif for mod2pi  # [-pi,pi]
    M = mod2pi(mean_anom)
    E = ecc_anom_init_guess_danby(M,ecc)
    const max_its_laguerre = 200
    for i in 1:max_its_laguerre
       E_old = E
       E = update_ecc_anom_laguerre(E_old, M, ecc)
       if abs(E-E_old)<tol break end
    end
    @assert abs(E-E_old)<tol
    return E
end

# Set the default Kepler solver
calc_ecc_anom = calc_ecc_anom_itterative_laguerre
#=
 function calc_ecc_anom{T}(mean_anom::T, ecc::T; tol::Real = 1.e-8)
  return calc_ecc_anom_itterative_laguerre(mean_anom, ecc, tol=tol)
end
=#

# export calc_ecc_anom


text:
#=
#define FMPR_RND_DOWN  0    RoundingMode{:Down}()
#define FMPR_RND_UP    1
#define FMPR_RND_FLOOR 2
#define FMPR_RND_CEIL  3
#define FMPR_RND_NEAR  4
=#

bits_to_rounded_digits(bits::Int32) = cld((bits * 3010%Int32), 10000%Int32)
bits_to_rounded_digits(bits::Int64) = bits_to_rounded_digits(bits%Int32)
digits_to_rounded_bits(digs::Int32) = fld((digs * 10000%Int32), 3010%Int32) + 1%Int32
digits_to_rounded_bits(digs::Int64) = digits_to_rounded_bits(digs%Int32)

integral_digits(x::ArbFloat{P}) where {P} = ceiled(Int, log10(1+floored(abs(x))))

function round(x::ArbFloat{P}, places::Int=integral_digits(P), base::Int=10) where {P}
    ((base==2) | (base==10)) || throw(ErrorException(string("Expecting base in (2,10), radix ",base," is not supported.")))
    sigbits = base==2 ? places+digits_to_rounded_bits(integral_digits(x)) : digits_to_rounded_bits(places+integral_digits(x))
    sigbits = max(1, sigbits) # library call chokes on a value of zero
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_set_round), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    return z
end

function ceiled(::Type{T}, x::ArbFloat{P}) where {P,T}
    y = ceiled(x)
    return convert(T,y)
end
function ceiled(x::ArbFloat{P}, places::Int=P, base::Int=2) where {P}
    ((base==2) | (base==10)) || throw(ErrorException(string("Expecting base in (2,10), radix ",base," is not supported.")))
    places = max(1,abs(places))
    sigbits = base==2 ? places : digits_to_rounded_bits(places)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_ceil), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    return z
end
function floored(::Type{T}, x::ArbFloat{P}) where {P,T}
    y = floored(x)
    return convert(T,y)
end
function floored(x::ArbFloat{P}, places::Int=P, base::Int=2) where {P}
    ((base==2) | (base==10)) || throw(ErrorException(string("Expecting base in (2,10), radix ",base," is not supported.")))
    places = max(1,abs(places))
    sigbits = base==2 ? places : digits_to_rounded_bits(places)
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_floor), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    return z
end

function ceil(x::ArbFloat{P}, places::Int=P, base::Int=2) where {P}
    ((base==2) | (base==10)) || throw(ErrorException(string("Expecting base in (2,10), radix ",base," is not supported.")))
    sigbits = base==2 ? places+digits_to_rounded_bits(integral_digits(x)) : digits_to_rounded_bits(places+integral_digits(x))
    sigbits = max(1, sigbits) # library call chokes on a value of zero
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_ceil), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    return z
end

function floor(x::ArbFloat{P}, places::Int=P, base::Int=2) where {P}
    ((base==2) | (base==10)) || throw(ErrorException(string("Expecting base in (2,10), radix ",base," is not supported.")))
    sigbits = base==2 ? places+digits_to_rounded_bits(integral_digits(x)) : digits_to_rounded_bits(places+integral_digits(x))
    sigbits = max(1, sigbits) # library call chokes on a value of zero
    z = initializer(ArbFloat{P})
    ccall(@libarb(arb_floor), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    return z
end

function trunc(x::ArbFloat{P}, places::Int=P, base::Int=2) where {P}
    ((base==2) | (base==10)) || throw(ErrorException(string("Expecting base in (2,10), radix ",base," is not supported.")))
    sigbits = base==2 ? places+digits_to_rounded_bits(integral_digits(x)) : digits_to_rounded_bits(places+integral_digits(x))
    sigbits = max(1, sigbits) # library call chokes on a value of zero
    z = initializer(ArbFloat{P})
    if signbit(x)
        ccall(@libarb(arb_ceil), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    else
        ccall(@libarb(arb_floor), Cvoid,  (Ref{ArbFloat{P}}, Ref{ArbFloat{P}}, Int), z, x, sigbits)
    end
    return z
end

function round(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=2) where {I <: Integer,P}
    z = round(x, digits=sig, base=base)
    return convert(I, z)
end
function ceil(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=2) where {I <: Integer,P}
    z = ceil(x, digits=sig, base=base)
    return convert(I, z)
end
function floor(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=2) where {I <: Integer,P}
    z = floor(x, digits=sig, base=base)
    return convert(I, z)
end
function trunc(::Type{I}, x::ArbFloat{P}, sig::Int=P, base::Int=2) where {I <: Integer,P}
    z = trunc(x, digits=sig, base=base)
    return convert(I, z)
end

fld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = convert(Int, floor(x/y))
cld(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = convert(Int, ceil(x/y))
div(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = convert(Int, trunc(x/y))

rem(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = convert(Int, x - div(x,y)*y)
mod(x::ArbFloat{P}, y::ArbFloat{P}) where {P} = convert(Int, x - fld(x,y)*y)

function divrem(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
   dv = div(x,y)
   r  = x - d*y
   rm = convert(Int, r)
   return dv,rm
end

function fldmod(x::ArbFloat{P}, y::ArbFloat{P}) where {P}
   fd = fld(x,y)
   m  = x - d*y
   md = convert(Int, m)
   return fd,md
end


text:
# ------------------------------------------------------------------------------------------
# # ModInt: a simple modular integer type
# ------------------------------------------------------------------------------------------

struct ModInt{n} <: Integer
    k::Int

    # Constructor definition...
    # note the left side looks like the call it defines
    ModInt{n}(k::Int) where {n} = new(mod(k,n))
end

a = ModInt{13}(1238279873492834)

b = ModInt{13}(9872349827349827)

a + b

# ------------------------------------------------------------------------------------------
# To extend standard functions we need to import them.
# ------------------------------------------------------------------------------------------

import Base: +

+(a::ModInt{n}, b::ModInt{n}) where {n} = ModInt{n}(a.k + b.k)

a + b

import Base: *, -

*(a::ModInt{n}, b::ModInt{n}) where{n} = ModInt{n}(a.k * b.k)
-(a::ModInt{n}, b::ModInt{n}) where {n} = ModInt{n}(a.k - b.k)
-(a::ModInt{n}) where {n} = ModInt{n}(-a.k)

a * b

a - b

-b

Base.show(io::IO, a::ModInt{n}) where {n} =
    get(io, :compact, false) ? show(io, a.k) : print(io, "$(a.k) mod $n")

a

b

a + 1

Base.promote_rule(::Type{ModInt{n}}, ::Type{Int}) where {n} = ModInt{n}

Base.convert(::Type{ModInt{n}}, x::Int) where {n} = ModInt{n}(x)

a + 1

1 + a

A = map(ModInt{13}, rand(1:100, 5, 5))

A^100000000

2A^100 .- 1

# ------------------------------------------------------------------------------------------
# ### Summary
#
# Here is all the code that defines the `ModInt` type:
# ```jl
# struct ModInt{n} <: Integer
#     k::Int
#
#     # Constructor definition...
#     # note the left side looks like the call it defines
#     ModInt{n}(k::Int) where {n} = new(mod(k,n))
# end
#
# import Base: +, *, -
#
# +(a::ModInt{n}, b::ModInt{n}) where {n} = ModInt{n}(a.k + b.k)
# *(a::ModInt{n}, b::ModInt{n}) where{n} = ModInt{n}(a.k * b.k)
# -(a::ModInt{n}, b::ModInt{n}) where {n} = ModInt{n}(a.k - b.k)
# -(a::ModInt{n}) where {n} = ModInt{n}(-a.k)
#
# Base.show(io::IO, a::ModInt{n}) where {n} =
#     get(io, :compact, false) ? show(io, a.k) : print(io, "$(a.k) mod $n")
#
# Base.promote_rule(::Type{ModInt{n}}, ::Type{Int}) where {n} = ModInt{n}
# Base.convert(::Type{ModInt{n}}, x::Int) where {n} = ModInt{n}(x)
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Exercise
#
# Add two methods that allows operations between modular integers with different modulus
# using the rule that they should combine in the modulus that is the `lcm` (least common
# multiple) of the moduli of the arguments.
#
# **Hint:** try something, see what fails, define something to make that work.
# ------------------------------------------------------------------------------------------

x = ModInt{12}(9)

y = ModInt{15}(13)

# two method definitions here...

@assert x + y == ModInt{60}(22)
@assert x * y == ModInt{60}(57)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
@testset "1104.path-in-zigzag-labelled-binary-tree.jl" begin
    @test path_in_zig_zag_tree(26) == [1, 2, 6, 10, 26]
    @test path_in_zig_zag_tree(14) == [1, 3, 4, 14]
    @test path_in_zig_zag_tree(7684) ==
        [1, 2, 7, 8, 30, 35, 120, 143, 480, 575, 1921, 2301, 7684]
end


text:
# This is an experimental build of SuiteSparse.

using BinaryBuilder, Pkg

name = "SuiteSparse_GPU"
version = v"5.8.1"

# Collection of sources required to build SuiteSparse
sources = [
    GitSource("https://github.com/DrTimothyAldenDavis/SuiteSparse.git",
              "1869379f464f0f8dac471edb4e6d010b2b0e639d"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SuiteSparse

# Apply Jameson's shlib patch
atomic_patch -p1 ${WORKSPACE}/srcdir/patches/SuiteSparse-shlib.patch

# Disable OpenMP as it will probably interfere with blas threads and Julia threads
FLAGS+=(INSTALL="${prefix}" INSTALL_LIB="${libdir}" INSTALL_INCLUDE="${prefix}/include" CFOPENMP=)

if [[ ${target} == *mingw32* ]]; then
    FLAGS+=(UNAME=Windows)
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir} -shared")
else
    FLAGS+=(UNAME="$(uname)")
    FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir}")
fi

# OpenBLAS
#if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
#    BLAS_64="-DSUN64 -DBLAS64 -DLONGBLAS='long long'"
#    BLAS_NAME=openblas64_
#else
#    BLAS_NAME=openblas
#fi
#FLAGS+=(BLAS="-l${BLAS_NAME}" LAPACK="-l${BLAS_NAME}")
#FLAGS+=(UMFPACK_CONFIG="$BLAS_64" CHOLMOD_CONFIG+="$BLAS_64" SPQR_CONFIG="$BLAS_64")

# MKL
FLAGS+=(MKLROOT="${prefix}")

# CUDA
FLAGS+=(CUDA_PATH="$prefix/cuda")

# METIS
FLAGS+=(MY_METIS_LIB="-lmetis" MY_METIS_INC="${prefix}/include")
# Disable METIS in CHOLMOD by passing -DNPARTITION and avoiding linking metis
#FLAGS+=(CHOLMOD_CONFIG+="-DNPARTITION")

make -j${nproc} -C SuiteSparse_config "${FLAGS[@]}" library config

for proj in SuiteSparse_config SuiteSparse_GPURuntime GPUQREngine AMD BTF CAMD CCOLAMD COLAMD CHOLMOD LDL KLU UMFPACK RBio SPQR SLIP_LU Mongoose; do
    make -j${nproc} -C $proj "${FLAGS[@]}" library CFOPENMP="$CFOPENMP"
    make -j${nproc} -C $proj "${FLAGS[@]}" install CFOPENMP="$CFOPENMP"
done

# For now, we'll have to adjust the name of the OpenBLAS library on macOS and FreeBSD.
# Eventually, this should be fixed upstream
if [[ ${target} == *-apple-* ]] || [[ ${target} == *freebsd* ]]; then
    echo "-- Modifying library name for OpenBLAS"

    for nm in libcholmod libspqr libumfpack; do
        # Figure out what version it probably latched on to:
        if [[ ${target} == *-apple-* ]]; then
            OPENBLAS_LINK=$(otool -L ${libdir}/${nm}.dylib | grep lib${BLAS_NAME} | awk '{ print $1 }')
            install_name_tool -change ${OPENBLAS_LINK} @rpath/lib${BLAS_NAME}.dylib ${libdir}/${nm}.dylib
        elif [[ ${target} == *freebsd* ]]; then
            OPENBLAS_LINK=$(readelf -d ${libdir}/${nm}.so | grep lib${BLAS_NAME} | sed -e 's/.*\[\(.*\)\].*/\1/')
            patchelf --replace-needed ${OPENBLAS_LINK} lib${BLAS_NAME}.so ${libdir}/${nm}.so
        fi
    done
fi

# Delete the extra soversion libraries built. https://github.com/JuliaPackaging/Yggdrasil/issues/7
if [[ "${target}" == *-mingw* ]]; then
    rm -f ${libdir}/lib*.*.${dlext}
    rm -f ${libdir}/lib*.*.*.${dlext}
fi

install_license LICENSE.txt

# Compile SuiteSparse_wrapper shim
cd $WORKSPACE/srcdir/SuiteSparse_wrapper
"${CC}" -O2 -shared -fPIC -I${prefix}/include SuiteSparse_wrapper.c -o ${libdir}/libsuitesparse_wrapper.${dlext} -L${libdir} -lcholmod
"""

platforms = [
    Platform("x86_64", "linux"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct("libsuitesparseconfig",      :libsuitesparseconfig),
    LibraryProduct("libamd",                    :libamd),
    LibraryProduct("libbtf",                    :libbtf),
    LibraryProduct("libcamd",                   :libcamd),
    LibraryProduct("libccolamd",                :libccolamd),
    LibraryProduct("libcolamd",                 :libcolamd),
    LibraryProduct("libcholmod",                :libcholmod),
    LibraryProduct("libldl",                    :libldl),
    LibraryProduct("libklu",                    :libklu),
    LibraryProduct("libumfpack",                :libumfpack),
    LibraryProduct("librbio",                   :librbio),
    LibraryProduct("libspqr",                   :libspqr),
    LibraryProduct("libsliplu",                 :libsliplu),
#    LibraryProduct("libmongoose",               :libmongoose),
    LibraryProduct("libGPUQREngine",            :libGPUQREngine),
    LibraryProduct("libSuiteSparse_GPURuntime", :libSuiteSparse_GPURuntime),
    LibraryProduct("libsuitesparse_wrapper",    :libsuitesparse_wrapper),
]

# Dependencies that must be installed before this package can be built
cuda_version = v"9.0.176"
dependencies = [
    Dependency("MKL_jll"),
    Dependency("METIS_jll"),
    Dependency("MPFR_jll"),
    Dependency("GMP_jll"),
    BuildDependency(PackageSpec(name="CUDA_full_jll", version=cuda_version))
]

# Note: we explicitly lie about this because we don't have the new
# versioning APIs worked out in BB yet.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6")

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests of back-end interface.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Back-end auto selection" begin
    data = rand(3,3)

    # Text
    # ==========================================================================

    auto   = pretty_table(String, data, tf = tf_unicode)
    manual = pretty_table(String, data, tf = tf_unicode, backend = :text)

    @test auto == manual

    # HTML
    # ==========================================================================

    auto   = pretty_table(String, data, tf = tf_html_simple)
    manual = pretty_table(String, data, tf = tf_html_simple, backend = :html)

    @test auto == manual

    # LaTeX
    # ==========================================================================

    auto   = pretty_table(String, data, tf = tf_latex_default)
    manual = pretty_table(String, data, tf = tf_latex_default, backend = :latex)

    @test auto == manual

    # Error
    # ==========================================================================

    @test_throws TypeError pretty_table(data, tf = [])
end


text:
@testset "pathof_noload" begin
    pnl = CompileBot.pathof_noload("Example")
    import Example
    p = GoodPath(pathof(Example))
    @test p == pnl
end

################################################################

@testset "detectOS" begin
    if Base.Sys.iswindows()
        @test ("windows", Base.Sys.iswindows) == CompileBot.detectOS()
    elseif Base.Sys.islinux()
        @test ("linux", Base.Sys.islinux) == CompileBot.detectOS()
    elseif Base.Sys.isapple()
        @test ("apple", Base.Sys.isapple) == CompileBot.detectOS()
    end
end
################################################################
# JuliaVersionNumber

# https://github.com/JuliaLang/julia/pull/36223:
# @test CompileBot.JuliaVersionNumber("nightly") ==
      # VersionNumber(replace(Base.read("VERSION", String), "\n" => ""))
# @test thispatch(CompileBot.JuliaVersionNumber("nightly")) == thispatch(VERSION)
@test CompileBot.JuliaVersionNumber("1.2.3") == v"1.2.3"
@test CompileBot.JuliaVersionNumber(v"1.2.3") == v"1.2.3"