text:
using BenchmarkTools,Knet

const N=1000

const sizes = (1,10,100,128,512,1000,1024,2048)

function f01(a,b)
    c=similar(a)
    n=length(a)
    for i=1:N; ccall(("add_32_01",Knet.libknet8),Nothing,(Cint,Float32,Ptr{Float32},Ptr{Float32}),n,b,a,c); end
    Knet.cudaDeviceSynchronize()
end

function f12(x,y)

    (dz,sx,nx,sy,ny,xlast,ylast,xdims,ydims,multi) = Knet.vbroadcast_shape(x,y)
    z = similar(x,dz)
    nz = length(z)
    for i=1:N; ccall(("add_32_12",Knet.libknet8),Nothing,(Cint,Ptr{Float32},Cint,Cint,Ptr{Float32},Cint,Cint,Ptr{Float32}),nz,x,sx,nx,y,sy,ny,z); end
    Knet.cudaDeviceSynchronize()
end

# other than first dim broadcast
# y is vector to be broadcasted, then ylast is broadcasted dim
function f13_x_y(x,y)
    (dz,sx,nx,sy,ny,xlast,ylast,xdims,ydims,multi) = Knet.vbroadcast_shape(x,y)
    z = similar(x,dz)
    brdcastdimstride = strides(x)[ylast]
    # if broadcast last dimension, nextstride is zero
    brdcastnextstride = ((ylast+1) > ndims(x) ? 0 : strides(x)[ylast+1])
    multidimsize = prod(size(x)[ylast+1:end])
    for i=1:N; ccall(("add_32_13_x_y",Knet.libknet8),Nothing,(Ptr{Float32},Ptr{Float32},Ptr{Float32},Cint,Cint,Cint,Cint,Cint),x,y,z,brdcastdimstride,brdcastnextstride,multidimsize,length(x),length(y)) end
    Knet.cudaDeviceSynchronize()
end

#x is N-dim y is vector
function f14_x_y(x,y)
    (dz,sx,nx,sy,ny,xlast,ylast,xdims,ydims,multi) = Knet.vbroadcast_shape(x,y)
    z = similar(x,dz)
    flat_dimsize=(length(x)/length(y))
    for i=1:N; ccall(("add_32_14_x_y",Knet.libknet8),Nothing,(Ptr{Float32},Ptr{Float32},Ptr{Float32},Cint,Cint,Cint),x,y,z,length(y),length(x),flat_dimsize) end
    Knet.cudaDeviceSynchronize()
end

#r=1 f13_x_y, r=2 f14
# m=nrows n=ncols
for r in (0,1,2)
    println(r==0 ? "a[m,n].+b" : r==1 ? "a[m,n].+b[1,n]" : r==2 ? "a[m,n].+b[m,1]" : error())
    for s in sizes; print("\t$s"); end; println()
    for nrows in sizes
        print(nrows)
        for ncols in sizes
            a = KnetArray(rand(Float32,nrows,ncols))
            b = (r==0 ? rand(Float32) : r==1 ? KnetArray(rand(Float32,1,ncols)) : KnetArray(rand(Float32,nrows,1)))
            # bm = (r==0 ? (@benchmark f01($a,$b) seconds=1) : (@benchmark f12($a,$b) seconds=1))
            yvectorFlag= (nrows>=2048 && (100<=ncols<128 )) || (nrows>=512 && (128<=ncols<512 )) || (nrows>=100 && (512<=ncols ))
            if r==0
              bm=(@benchmark f01($a,$b) seconds=1)
            elseif (r==1 && (ncols<704 || nrows<512)) || (r==2)
            # elseif (r==1 && (ncols<704 || nrows<512)) || (r==2 && (!yvectorFlag))
              bm=(@benchmark f12($a,$b) seconds=1)
            elseif r==1
              bm=(@benchmark f13_x_y($a,$b) seconds=1)
            # else
            #   bm=(@benchmark f14_x_y($a,$b) seconds=1)
            end
            m = (round(Int, minimum(bm.times)/N))
            #
            # m = (ncols*nrows*4)/(round(Int, minimum(bm.times)/N))
            print("\t$m")
            a=b=nothing; Knet.gc()
        end
        println()
    end
end

# commit 359d3646 2018-09-22, julia 1.0.0                               commit 4aa5f92f 2018-08-14, julia 0.6.4
# GPU:V100, CPU:Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz               GPU:V100, CPU:Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz
#
# a[m,n].+b                                                     	a[m,n].+b
# 	1	10	100	128	512	1000	1024	2048    	1	10	100	128	512	1000	1024	2048
# 1	3168	3007	3055	3321	3317	3321	3334	3316    1	3362	3349	3367	3353	3358	3350	3372	3645
# 10	3349	3290	3201	3186	3188	3199	3063	3081    10	3634	3629	3636	3622	3645	3633	3631	3571
# 100	3070	3067	3064	3074	3062	3236	3217	3232    100	3527	3365	3369	3348	3363	3379	3372	3384
# 128	3223	3342	3333	3352	3347	3338	3332	3329    128	3367	3370	3344	3372	3643	3631	3637	3656
# 512	3333	3294	3192	3348	3231	3776	3819	14590   512	3625	3622	3645	3653	3575	3714	3768	14632
# 1000	3219	3226	3210	3027	3760	13998	14180	26367   1000	3486	3488	3375	3360	3731	13972	14175	26436
# 1024	3312	3322	3413	3522	3807	14191	14671	26865   1024	3486	3632	3631	3624	3774	14196	14571	26947
# 2048	3225	3064	3219	3217	14637	26230	27076	51566   2048	3370	3540	3484	3336	14600	26448	27023	51247
# a[m,n].+b[1,n]                                                        a[m,n].+b[1,n]
# 	1	10	100	128	512	1000	1024	2048    	1	10	100	128	512	1000	1024	2048
# 1	3361	3548	3561	3369	3361	3533	3292	3245    1	3262	3279	3555	3551	3573	3369	3373	3380
# 10	3232	3245	3220	3250	3221	3060	3222	3231    10	3368	3376	3382	3274	3086	3093	3082	3082
# 100	3242	3236	3245	3253	3550	3539	3558	4156    100	3096	3076	3106	3092	3396	3359	3381	4081
# 128	3554	3515	3491	3466	3378	3375	3378	4548    128	3475	3376	3378	3398	3380	3369	3315	4507
# 512	3239	3254	3237	3253	4546	4105	4106	12231   512	3298	3322	3266	3072	4502	4056	4054	12364
# 1000	3241	3552	3554	3561	6783	15565	15766	25255   1000	3366	3388	3372	3370	6757	15579	15845	25196
# 1024	3544	3390	3245	3297	7168	14399	14576	24069   1024	3101	3090	3084	3252	7094	14324	14512	24043
# 2048	3244	3261	4134	4536	20508	26210	26678	45702   2048	3098	3090	4098	4524	20345	26623	26452	45965
# a[m,n].+b[m,1]                                                        a[m,n].+b[m,1]
# 	1	10	100	128	512	1000	1024	2048    	1	10	100	128	512	1000	1024	2048
# 1	3568	3569	3548	3544	3503	3396	3391	3390    1	3355	3355	3350	3353	3352	3351	3289	3293
# 10	3388	3234	3249	3256	3230	3237	3248	3235    10	3240	3254	3298	3071	3065	3071	3084	3075
# 100	3237	3236	3237	3250	3247	3546	3556	3783    100	3085	3350	3382	3380	3558	3583	3581	3743
# 128	3549	3555	3500	3484	3473	3372	3244	4016    128	3347	3311	3261	3244	3234	3234	3251	3968
# 512	3258	3237	3240	3237	4017	5433	5616	16587   512	3134	3150	3127	3116	3959	5406	5592	16779
# 1000	3247	3247	3552	3567	5990	16724	16946	30295   1000	3616	3592	3301	3287	6003	16639	16962	30431
# 1024	3554	3563	3059	3234	5418	16085	16515	30598   1024	3303	3285	3292	3300	5400	16201	16470	30657
# 2048	3244	3240	3624	3814	16551	29967	30636	57885   2048	3314	3606	3617	3746	16658	30043	30680	58019

# COMMIT 3d32e16 2017-05-17
#
# a[m,n].+b
# 	1	10	100	128	512	1000	1024	2048
# 1	4264	4287	4295	4292	4302	4305	4306	4337
# 10	4288	4293	4305	4315	4365	4444	4459	4681
# 100	4295	4305	4457	4469	5209	6192	6238	11079
# 128	4299	4314	4472	4560	5600	6701	6787	17908
# 512	4303	4375	5211	5600	17957	30842	31749	59719
# 1000	4313	4436	6201	6699	30841	57163	58332	113552
# 1024	4312	4459	6245	6797	31668	58113	59591	116099
# 2048	4341	4683	11695	18025	59673	113487	116187	232430
# a[m,n].+b[1,n]
# 	1	10	100	128	512	1000	1024	2048
# 1	4304	4308	4317	4316	4330	4356	4356	4361
# 10	4307	4487	4576	4576	4706	5018	5104	5523
# 100	4315	4575	5006	5204	7389	10042	10129	16564
# 128	4316	4576	5202	5317	8237	11773	12162	22033
# 512	4325	4707	7414	8240	21068	33924	34016	60985
# 1000	4357	5006	10047	11783	36528	72782	73926	137811
# 1024	4355	5100	10153	12170	38988	61619	62400	115017
# 2048	4361	5531	16586	20558	76857	118232	119759	227061
# a[m,n].+b[m,1]
# 	1	10	100	128	512	1000	1024	2048
# 1	4311	4311	4316	4316	4325	4355	4354	4359
# 10	4307	4314	4366	4371	4696	5065	5075	5873
# 100	4318	4406	4897	5040	7304	10312	10450	17811
# 128	4318	4366	4861	5018	7302	10195	10441	18789
# 512	4334	4484	6512	7154	18362	31899	32853	61838
# 1000	4356	4729	8474	9675	32089	59411	60838	119014
# 1024	4356	4732	8459	9856	32833	60652	62001	121412
# 2048	4356	5196	13393	18371	61866	118420	121254	242820

# COMMIT: ccad0cb 2017-04-15
# a[m,n].+b
# 	1	10	100	1000
# 1	4264	4272	4281	4299
# 10	4268	4277	4298	4424
# 100	4275	4299	4422	6173
# 1000	4300	4412	6173	57050
# a[m,n].+b[1,n]
# 	1	10	100	1000
# 1	4296	4295	4312	4354
# 10	4304	4481	4573	5001
# 100	4311	4573	4992	10034
# 1000	4345	4998	10047	69640
# a[m,n].+b[m,1]
# 	1	10	100	1000
# 1	4304	4306	4314	4349
# 10	4300	4313	4359	5050
# 100	4314	4402	4865	10295
# 1000	4354	4718	8433	59328


text:
module TestShow
    using DataFrames
    using Compat
    using Base.Test
    import Compat.String
    df = DataFrame(A = 1:3, B = ["x", "y", "z"])

    io = IOBuffer()
    show(io, df)
    show(io, df, true)
    showall(io, df)
    showall(io, df, true)

    subdf = sub(df, [2, 3]) # df[df[:A] .> 1.0, :]
    show(io, subdf)
    show(io, subdf, true)
    showall(io, subdf)
    showall(io, subdf, true)

    dfvec = DataFrame[df for _=1:3]
    show(io, dfvec)
    showall(io, dfvec)

    gd = groupby(df, :A)
    show(io, gd)
    showall(io, gd)

    dfr = DataFrameRow(df, 1)
    show(io, dfr)

    df = DataFrame(A = Array(String, 3))

    A = DataFrames.StackedVector(Any[[1, 2, 3], [4, 5, 6], [7, 8, 9]])
    show(io, A)
    A = DataFrames.RepeatedVector([1, 2, 3], 5, 1)
    show(io, A)
    A = DataFrames.RepeatedVector([1, 2, 3], 1, 5)
    show(io, A)

    #Test show output for REPL and similar
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, Nullable()])
    io = IOBuffer()
    show(io, df)
    str = takebuf_string(io)
    @test str == """
2×2 DataFrames.DataFrame
│ Row │ Fish │ Mass  │
├─────┼──────┼───────┤
│ 1   │ Suzy │ 1.5   │
│ 2   │ Amir │ #NULL │"""
end

ext: jl


text:
"""
    pass2(project)

 - Remove `pass2/`
 - Copy files from `pass1/` into `pass2`
 - Remove comments (if indicated)
"""
function pass2(project)
    println("Performing pass 2...")
    files = prep_pass(project, 2)
    if !project.preserve_comments
        for f in files
            remove_comments(f, project.language_from)
        end
    end
end


text:
# ============================================================
# Geometric Methods
# ============================================================

include("geo_transform.jl")
include("geo_points.jl")
include("geo_jacobi.jl")
include("geo_neighbor.jl")


text:
using ArrayLayouts, LinearAlgebra, FillArrays, Test
import ArrayLayouts: MemoryLayout, DenseRowMajor, DenseColumnMajor, StridedLayout,
                        ConjLayout, RowMajor, ColumnMajor, UnitStride,
                        SymmetricLayout, HermitianLayout, UpperTriangularLayout,
                        UnitUpperTriangularLayout, LowerTriangularLayout,
                        UnitLowerTriangularLayout, ScalarLayout, UnknownLayout,
                        hermitiandata, symmetricdata, FillLayout, ZerosLayout,
                        DiagonalLayout, TridiagonalLayout, SymTridiagonalLayout, colsupport, rowsupport,
                        diagonaldata, subdiagonaldata, supdiagonaldata, BidiagonalLayout, bidiagonaluplo

struct FooBar end
struct FooNumber <: Number end

@testset "MemoryLayout" begin
    @testset "Trivial" begin
        @test MemoryLayout(3.141) == MemoryLayout(Float64) == MemoryLayout(Int) == MemoryLayout(FooNumber) == ScalarLayout()
        @test MemoryLayout(FooBar()) == MemoryLayout(FooBar) == UnknownLayout()

        A = randn(6)
        @test MemoryLayout(A) == MemoryLayout(Base.ReshapedArray(A,(2,3),())) ==
            MemoryLayout(reinterpret(Float32,A)) == MemoryLayout(A) == DenseColumnMajor()

        @test MemoryLayout(view(A,1:3)) == DenseColumnMajor()
        @test MemoryLayout(view(A,Base.OneTo(3))) == DenseColumnMajor()
        @test MemoryLayout(view(A,:)) == DenseColumnMajor()
        @test MemoryLayout(view(A,CartesianIndex(1,1))) == DenseColumnMajor()
        @test MemoryLayout(view(A,1:2:4)) == StridedLayout()

        A = randn(6,6)
        V = view(A, 1:3,:)
        @test MemoryLayout(V) == MemoryLayout(V) == ColumnMajor()
    end

    @testset "adjoint and transpose MemoryLayout" begin
        A = [1.0 2; 3 4]
        @test MemoryLayout(A') == DenseRowMajor()
        @test MemoryLayout(transpose(A)) == DenseRowMajor()
        B = [1.0+im 2; 3 4]
        @test MemoryLayout(B') == ConjLayout{DenseRowMajor}()
        @test MemoryLayout(transpose(B)) == DenseRowMajor()
        VA = view(A, 1:1, 1:1)
        @test MemoryLayout(VA') == RowMajor()
        @test MemoryLayout(transpose(VA)) == RowMajor()
        VB = view(B, 1:1, 1:1)
        @test MemoryLayout(VB') == ConjLayout{RowMajor}()
        @test MemoryLayout(transpose(VB)) == RowMajor()
        VA = view(A, 1:2:2, 1:2:2)
        @test MemoryLayout(VA') == StridedLayout()
        @test MemoryLayout(transpose(VA)) == StridedLayout()
        VB = view(B, 1:2:2, 1:2:2)
        @test MemoryLayout(VB') == ConjLayout{StridedLayout}()
        @test MemoryLayout(transpose(VB)) == StridedLayout()
        VA2 = view(A, [1,2], :)
        @test MemoryLayout(VA2') == UnknownLayout()
        @test MemoryLayout(transpose(VA2)) == UnknownLayout()
        VB2 = view(B, [1,2], :)
        @test MemoryLayout(VB2') == UnknownLayout()
        @test MemoryLayout(transpose(VB2)) == UnknownLayout()
        VA2 = view(A, 1:2, :)
        @test MemoryLayout(VA2') == RowMajor()
        @test MemoryLayout(transpose(VA2)) == RowMajor()
        VB2 = view(B, 1:2, :)
        @test MemoryLayout(VB2') == ConjLayout{RowMajor}()
        @test MemoryLayout(transpose(VB2)) == RowMajor()
        VA2 = view(A, :, 1:2)
        @test MemoryLayout(VA2') == DenseRowMajor()
        @test MemoryLayout(transpose(VA2)) == DenseRowMajor()
        VB2 = view(B, :, 1:2)
        @test MemoryLayout(VB2') == ConjLayout{DenseRowMajor}()
        @test MemoryLayout(transpose(VB2)) == DenseRowMajor()
        VAc = view(A', 1:1, 1:1)
        @test MemoryLayout(VAc) == RowMajor()
        VAt = view(transpose(A), 1:1, 1:1)
        @test MemoryLayout(VAt) == RowMajor()
        VBc = view(B', 1:1, 1:1)
        @test MemoryLayout(VBc) == ConjLayout{RowMajor}()
        VBt = view(transpose(B), 1:1, 1:1)
        @test MemoryLayout(VBt) == RowMajor()

        @test MemoryLayout(view(randn(5)',[1,3])) == UnknownLayout()
    end

    @testset "Bi/Tridiagonal" begin
        T = Tridiagonal(randn(5),randn(6),randn(5))
        S = SymTridiagonal(T.d, T.du)
        Bl = Bidiagonal(T.d, T.dl, :L)
        Bu = Bidiagonal(T.d, T.du, :U)

        @test MemoryLayout(T) isa TridiagonalLayout
        @test MemoryLayout(Adjoint(T)) isa TridiagonalLayout
        @test MemoryLayout(Transpose(T)) isa TridiagonalLayout
        @test MemoryLayout(S) isa SymTridiagonalLayout
        @test MemoryLayout(Adjoint(S)) isa SymTridiagonalLayout
        @test MemoryLayout(Transpose(S)) isa SymTridiagonalLayout
        @test MemoryLayout(Bl) isa BidiagonalLayout
        @test MemoryLayout(Adjoint(Bl)) isa BidiagonalLayout
        @test MemoryLayout(Transpose(Bl)) isa BidiagonalLayout
        @test MemoryLayout(Bu) isa BidiagonalLayout
        @test MemoryLayout(Adjoint(Bu)) isa BidiagonalLayout
        @test MemoryLayout(Transpose(Bu)) isa BidiagonalLayout

        @test bidiagonaluplo(Bl) == bidiagonaluplo(Adjoint(Bu)) == 'L'
        @test bidiagonaluplo(Bu) == bidiagonaluplo(Adjoint(Bl)) == 'U'

        @test diagonaldata(T) == diagonaldata(T') == diagonaldata(S) == diagonaldata(Bl) == diagonaldata(Bu)
        @test supdiagonaldata(T) == subdiagonaldata(Adjoint(T)) == subdiagonaldata(Transpose(T)) ==
                    supdiagonaldata(S) == subdiagonaldata(S) ==
                    supdiagonaldata(Bu) == subdiagonaldata(Adjoint(Bu)) == subdiagonaldata(Transpose(Bu))
        @test subdiagonaldata(T) == supdiagonaldata(Adjoint(T)) == supdiagonaldata(Transpose(T)) ==
                subdiagonaldata(Bl) == supdiagonaldata(Adjoint(Bl)) == supdiagonaldata(Transpose(Bl)) ==
                T.dl

        @test colsupport(T,3) == rowsupport(T,3) == colsupport(S,3) == rowsupport(S,3) == 2:4
        @test colsupport(T,3:6) == rowsupport(T,3:6) == colsupport(S,3:6) == rowsupport(S,3:6) == 2:6
        @test colsupport(Bl,3) == rowsupport(Bu,3) == rowsupport(Adjoint(Bl),3) == 3:4
        @test rowsupport(Bl,3) == colsupport(Bu,3) == colsupport(Adjoint(Bl),3) == 2:3
        @test colsupport(Bl,3:6) == rowsupport(Bu,3:6) == 3:6
        @test colsupport(Bu,3:6) == rowsupport(Bl,3:6) == 2:6
    end

    @testset "Symmetric/Hermitian" begin
        A = [1.0 2; 3 4]
        @test MemoryLayout(Symmetric(A)) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Hermitian(A)) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Transpose(Symmetric(A))) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Transpose(Hermitian(A))) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Adjoint(Symmetric(A))) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Adjoint(Hermitian(A))) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(view(Symmetric(A),:,:)) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(view(Hermitian(A),:,:)) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Symmetric(A')) == SymmetricLayout{DenseRowMajor}()
        @test MemoryLayout(Hermitian(A')) == SymmetricLayout{DenseRowMajor}()
        @test MemoryLayout(Symmetric(transpose(A))) == SymmetricLayout{DenseRowMajor}()
        @test MemoryLayout(Hermitian(transpose(A))) == SymmetricLayout{DenseRowMajor}()

        @test symmetricdata(Symmetric(A)) ≡ A
        @test symmetricdata(Hermitian(A)) ≡ A
        @test symmetricdata(Transpose(Symmetric(A))) ≡ A
        @test symmetricdata(Transpose(Hermitian(A))) ≡ A
        @test symmetricdata(Adjoint(Symmetric(A))) ≡ A
        @test symmetricdata(Adjoint(Hermitian(A))) ≡ A
        @test symmetricdata(view(Symmetric(A),:,:)) ≡ A
        @test symmetricdata(view(Hermitian(A),:,:)) ≡ A
        @test symmetricdata(Symmetric(A')) ≡ A'
        @test symmetricdata(Hermitian(A')) ≡ A'
        @test symmetricdata(Symmetric(transpose(A))) ≡ transpose(A)
        @test symmetricdata(Hermitian(transpose(A))) ≡ transpose(A)

        @test colsupport(Symmetric(A),2) ≡ colsupport(Symmetric(A),1:2) ≡
                rowsupport(Symmetric(A),2) ≡ rowsupport(Symmetric(A),1:2) ≡ 1:2
                @test colsupport(Hermitian(A),2) ≡ colsupport(Hermitian(A),1:2) ≡
                rowsupport(Hermitian(A),2) ≡ rowsupport(Hermitian(A),1:2) ≡ 1:2

        B = [1.0+im 2; 3 4]
        @test MemoryLayout(Symmetric(B)) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Hermitian(B)) == HermitianLayout{DenseColumnMajor}()
        @test MemoryLayout(Transpose(Symmetric(B))) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(Transpose(Hermitian(B))) == UnknownLayout()
        @test MemoryLayout(Adjoint(Symmetric(B))) == UnknownLayout()
        @test MemoryLayout(Adjoint(Hermitian(B))) == HermitianLayout{DenseColumnMajor}()
        @test MemoryLayout(view(Symmetric(B),:,:)) == SymmetricLayout{DenseColumnMajor}()
        @test MemoryLayout(view(Hermitian(B),:,:)) == HermitianLayout{DenseColumnMajor}()
        @test MemoryLayout(Symmetric(B')) == UnknownLayout()
        @test MemoryLayout(Hermitian(B')) == UnknownLayout()
        @test MemoryLayout(Symmetric(transpose(B))) == SymmetricLayout{DenseRowMajor}()
        @test MemoryLayout(Hermitian(transpose(B))) == HermitianLayout{DenseRowMajor}()

        @test symmetricdata(Symmetric(B)) ≡ B
        @test hermitiandata(Hermitian(B)) ≡ B
        @test symmetricdata(Transpose(Symmetric(B))) ≡ B
        @test hermitiandata(Adjoint(Hermitian(B))) ≡ B
        @test symmetricdata(view(Symmetric(B),:,:)) ≡ B
        @test hermitiandata(view(Hermitian(B),:,:)) ≡ B
        @test symmetricdata(Symmetric(B')) ≡ B'
        @test hermitiandata(Hermitian(B')) ≡ B'
        @test symmetricdata(Symmetric(transpose(B))) ≡ transpose(B)
        @test hermitiandata(Hermitian(transpose(B))) ≡ transpose(B)

        @testset "Bidiagonal" begin
            B = Bidiagonal(randn(6),randn(5),:U)
            Bc = Bidiagonal(randn(6) .+ 0im,randn(5) .+ 1im,:U)
            S = Symmetric(B)
            H = Hermitian(B)
            Sc = Symmetric(Bc)
            Hc = Hermitian(Bc)

            @test MemoryLayout(S) isa SymTridiagonalLayout
            @test MemoryLayout(H) isa SymTridiagonalLayout
            @test MemoryLayout(Sc) isa SymTridiagonalLayout
            @test MemoryLayout(Hc) isa HermitianLayout

            @test diagonaldata(S) == diagonaldata(B)
            @test subdiagonaldata(S) == supdiagonaldata(S) == supdiagonaldata(B)

            @test colsupport(S,3) == colsupport(H,3) == colsupport(Sc,3) == colsupport(Hc,3) == 2:4
            @test rowsupport(S,3) == rowsupport(H,3) == rowsupport(Sc,3) == rowsupport(Hc,3) == 2:4
        end
     end

    @testset "triangular MemoryLayout" begin
        A = [1.0 2; 3 4]
        B = [1.0+im 2; 3 4]
        for (TriType, TriLayout, TriLayoutTrans) in ((UpperTriangular, UpperTriangularLayout, LowerTriangularLayout),
                              (UnitUpperTriangular, UnitUpperTriangularLayout, UnitLowerTriangularLayout),
                              (LowerTriangular, LowerTriangularLayout, UpperTriangularLayout),
                              (UnitLowerTriangular, UnitLowerTriangularLayout, UnitUpperTriangularLayout))
            @test MemoryLayout(TriType(A)) == TriLayout{DenseColumnMajor}()
            @test MemoryLayout(TriType(transpose(A))) == TriLayout{DenseRowMajor}()
            @test MemoryLayout(TriType(A')) == TriLayout{DenseRowMajor}()
            @test MemoryLayout(transpose(TriType(A))) == TriLayoutTrans{DenseRowMajor}()
            @test MemoryLayout(TriType(A)') == TriLayoutTrans{DenseRowMajor}()

            @test MemoryLayout(TriType(B)) == TriLayout{DenseColumnMajor}()
            @test MemoryLayout(TriType(transpose(B))) == TriLayout{DenseRowMajor}()
            @test MemoryLayout(TriType(B')) == TriLayout{ConjLayout{DenseRowMajor}}()
            @test MemoryLayout(transpose(TriType(B))) == TriLayoutTrans{DenseRowMajor}()
            @test MemoryLayout(TriType(B)') == TriLayoutTrans{ConjLayout{DenseRowMajor}}()
        end

        @test MemoryLayout(UpperTriangular(B)') == MemoryLayout(LowerTriangular(B'))
    end

    @testset "Reinterpreted/Reshaped" begin
       @test MemoryLayout(reinterpret(Float32, UInt32[1 2 3 4 5])) == DenseColumnMajor()
       @test MemoryLayout(reinterpret(Float32, UInt32[1 2 3 4 5]')) == UnknownLayout()
       @test MemoryLayout(Base.__reshape((randn(6),IndexLinear()),(2,3))) == DenseColumnMajor()
       @test MemoryLayout(Base.__reshape((1:6,IndexLinear()),(2,3))) == UnknownLayout()
    end

    @testset "Fill" begin
        @test MemoryLayout(Fill(1,10)) == FillLayout()
        @test MemoryLayout(Ones(10)) == FillLayout()
        @test MemoryLayout(Zeros(10)) == ZerosLayout()
        @test MemoryLayout(view(Fill(1,10),1:3)) == FillLayout()
        @test MemoryLayout(view(Fill(1,10),1:3,1)) == FillLayout()
    end

    @testset "Triangular col/rowsupport" begin
        A = randn(5,5)
        @test colsupport(UpperTriangular(A),3) ≡ Base.OneTo(3)
        @test rowsupport(UpperTriangular(A),3) ≡ 3:5
        @test colsupport(LowerTriangular(A),3) ≡ 3:5
        @test rowsupport(LowerTriangular(A),3) ≡ Base.OneTo(3)
    end

    @testset "PermutedDimsArray" begin
        A = [1.0 2; 3 4]
        @test MemoryLayout(PermutedDimsArray(A, (1,2))) == DenseColumnMajor()
        @test MemoryLayout(PermutedDimsArray(A, (2,1))) == DenseRowMajor()
        @test MemoryLayout(transpose(PermutedDimsArray(A, (2,1)))) == DenseColumnMajor()
        @test MemoryLayout(adjoint(PermutedDimsArray(A, (2,1)))) == DenseColumnMajor()
        B = [1.0+im 2; 3 4]
        @test MemoryLayout(PermutedDimsArray(B, (2,1))) == DenseRowMajor()
        @test MemoryLayout(transpose(PermutedDimsArray(B, (2,1)))) == DenseColumnMajor()
        @test MemoryLayout(adjoint(PermutedDimsArray(B, (2,1)))) == ConjLayout{DenseColumnMajor}()

        C = view(ones(10,20,30), 2:9, 3:18, 4:27);
        @test MemoryLayout(C) == ColumnMajor()
        @test MemoryLayout(PermutedDimsArray(C, (1,2,3))) == ColumnMajor()
        @test MemoryLayout(PermutedDimsArray(C, (1,3,2))) == UnitStride{1}()

        @test MemoryLayout(PermutedDimsArray(C, (3,1,2))) == UnitStride{2}()
        @test MemoryLayout(PermutedDimsArray(C, (2,1,3))) == UnitStride{2}()

        @test MemoryLayout(PermutedDimsArray(C, (3,2,1))) == RowMajor()
        @test MemoryLayout(PermutedDimsArray(C, (2,3,1))) == UnitStride{3}()

        revC = PermutedDimsArray(C, (3,2,1));
        @test MemoryLayout(PermutedDimsArray(revC, (3,2,1))) == ColumnMajor()
        @test MemoryLayout(PermutedDimsArray(revC, (3,1,2))) == UnitStride{1}()

        D = ones(10,20,30,40);
        @test MemoryLayout(D) == DenseColumnMajor()
        @test MemoryLayout(PermutedDimsArray(D, (1,2,3,4))) == DenseColumnMajor()
        @test MemoryLayout(PermutedDimsArray(D, (1,4,3,2))) == UnitStride{1}()

        @test MemoryLayout(PermutedDimsArray(D, (4,1,3,2))) == UnitStride{2}()
        @test MemoryLayout(PermutedDimsArray(D, (2,1,4,3))) == UnitStride{2}()

        @test MemoryLayout(PermutedDimsArray(D, (4,3,2,1))) == DenseRowMajor()
        @test MemoryLayout(PermutedDimsArray(D, (4,2,1,3))) == UnitStride{3}()

        twoD = PermutedDimsArray(D, (3,1,2,4));
        MemoryLayout(PermutedDimsArray(twoD, (2,1,4,3))) == UnitStride{1}()

        revD = PermutedDimsArray(D, (4,3,2,1));
        MemoryLayout(PermutedDimsArray(revD, (4,3,2,1))) == DenseColumnMajor()
        MemoryLayout(PermutedDimsArray(revD, (4,2,3,1))) == UnitStride{1}()

        issorted((1,2,3,4))
        # Fails on Julia 1.4, in tests. Could use BenchmarkTools.@ballocated instead.
        @test_skip 0 == @allocated issorted((1,2,3,4))
        reverse((1,2,3,4))
        @test_skip 0 == @allocated reverse((1,2,3,4))
        MemoryLayout(revD)
        @test 0 == @allocated MemoryLayout(revD)
    end
end


text:
abstract type AbstractGravity <: AbstractPhysicsComponent end

struct ThinShellGravityFromPotential <: AbstractGravity end
struct TotalEnergyGravityFromPotential <: AbstractGravity end

@Base.kwdef struct DeepShellGravity{FT} <: AbstractGravity
    g::FT # ms⁻²
    a::FT # m
end

@Base.kwdef struct ThinShellGravity{FT} <: AbstractGravity
    g::FT # ms⁻²
end

@Base.kwdef struct Buoyancy{FT} <: AbstractGravity
    α::FT # K⁻¹
    g::FT # ms⁻²
end

@inline calc_force!(source, ::Nothing, state, _...) = nothing

@inline function calc_force!(source, gravity::DeepShellGravity, state, aux, orientation, _...)
    g = gravity.g
    a = gravity.a
    ρ = state.ρ
    r⃗ = @SVector [aux.x, aux.y, aux.z]
    k̂ = vertical_unit_vector(orientation, aux)

    r = norm(r⃗)

    # TODO!: Need numerical gradient of geopotential
    source.ρu -= g * (a / r)^2 * k̂ * ρ

    return nothing
end

@inline function calc_force!(source, gravity::ThinShellGravity, state, aux, orientation, _...)
    g = gravity.g
    ρ = state.ρ
    k̂ = vertical_unit_vector(orientation, aux)

    source.ρu -= g * k̂ * ρ

    return nothing
end

@inline function calc_force!(source, gravity::Buoyancy, state, aux, orientation, _...)
    α = gravity.α
    g = gravity.g
    ρθ = state.ρθ
    k̂ = vertical_unit_vector(orientation, aux)

    source.ρu -= -α * g * k̂ * ρθ

    return nothing
end

@inline function calc_force!(source, ::ThinShellGravityFromPotential, state, aux, _...)
    ∇Φ = aux.∇Φ
    source.ρu -= state.ρ * ∇Φ
end

@inline function calc_force!(source, ::TotalEnergyGravityFromPotential, state, aux, _...)
    ∇Φ = aux.∇Φ
    source.ρe -= state.ρu' * ∇Φ
end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine

max_forks_repo_name: CliMA/CouplerMachine


text:
# ---
# title: 825. Friends Of Appropriate Ages
# id: problem825
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/friends-of-appropriate-ages/description/>
# hidden: true
# ---
#
# Some people will make friend requests. The list of their ages is given and
# `ages[i]` is the age of the ith person.
#
# Person A will NOT friend request person B (B != A) if any of the following
# conditions are true:
#
#   * `age[B] <= 0.5 * age[A] + 7`
#   * `age[B] > age[A]`
#   * `age[B] > 100 && age[A] < 100`
#
# Otherwise, A will friend request B.
#
# Note that if A requests B, B does not necessarily request A.  Also, people
# will not friend request themselves.
#
# How many total friend requests are made?
#
# **Example 1:**
#
#
#
#     Input: [16,16]
#     Output: 2
#     Explanation: 2 people friend request each other.
#
#
# **Example 2:**
#
#
#
#     Input: [16,17,18]
#     Output: 2
#     Explanation: Friend requests are made 17 -> 16, 18 -> 17.
#
# **Example 3:**
#
#
#
#     Input: [20,30,100,110,120]
#     Output: 3
#     Explanation: Friend requests are made 110 -> 100, 120 -> 110, 120 -> 100.
#
#
#
#
# Notes:
#
#   * `1 <= ages.length <= 20000`.
#   * `1 <= ages[i] <= 120`.
#
#
## @lc code=start
using LeetCode

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


text:
module TestReleaseFlow

using Test

@testset "$path" for path in [
    "test_replace_commits.jl"
]
    include(path)
end

end  # module


text:
## %% Global Settings Structure (to be set by modeller)

mutable struct STABLEglobalSettings
      Hours::Integer       # The number of hours.
      Timestep::Integer    # Time interval size as multiple of 30 minutes.
      DiscountRate::Float64    # The annual interest rate used to compute the net present value.
      # ReferenceYear::Integer # Year of data set to use for renewable traces
      CostScaling::Float64
end

## %% Scenario Structure

mutable struct STABLEscenario
      ScenarioSettings::STABLEglobalSettings
      ScenarioIdentifier::String
      ScenarioName::String
      ScenarioNumber::Integer     # Number identifier for data input/lookup.
      TGScenarioIdentifier::String    # Bridge to pre-populated scenarios.
      TGScenarioNumber::String        # Bridge to pre-populated scenarios.
      WeatherYear::Integer        # Historical year as basis for renewable traces
      TraceYear::Integer          # Which year to use from the ReferenceYear trace dataset
      BatteryImproveYear::Integer # Projected year of improved battery response speed
      # GridFollowFraction::Float64 # Fraction of installed renewables that have grid following inverters in BatteryImproveYear
      MinimumRenewShare::Integer  # Default value on the minimum yearly renewables requirement for the whole of system.
      TechCostScenario::String    # Scenario identifier for the technology cost data
      FuelCostScenario::String    # Scenario identifier for the fuel cost data
      VPPbatteryScenario          # Scenario identifier for the virtual power plant battery uptake data.
      HydrogenScenario            # Scenario identifier for hydrogen uptake.
      HydrogenGasPeaker::Integer  # Scenario identifier for allowing hydrogen gas peaking plant.
      FixExistingCap::Bool        # Flag to indicate whether to fix the variables for existing capacity to the given values if true.
      NewGasAllowed::Bool         # Flag to indicate whether new gas plant capacity is built if true.
      NewDistillateAllowed::Bool  # Flag to indicate whether new liquid distillate technology capacity is built if true.
      NoNewFossilFuelsYear::Integer # Projected year at which no new fossil-fuel technologies are built.
      TxScenario::Array{String}   # List of committed transmission upgrades.
end

## %% Instance Structure

mutable struct STABLEinstance
      InstanceScenario::STABLEscenario
      InstanceIdentifier::String
      InstanceYear::Integer # The year (projected or historical) to simulate within the given Scenario.
end
# e.g # STinst.Flags["FixExistingCapFlag"] = true, NoNewGas, NoNewDistillateAllowed

## %% Run Settings Structure for specific configuration of runs

mutable struct STABLErunSettings
      RunSettingsIdentifier::String #
      Flags::OrderedDict{String,Bool}  #
      Parameters::OrderedDict{String,Any}
end

## Transmission parameters:
# Tx_Loss_Factor::Float64  # Transmission loss factor multiplier on demand e.g. 1.06 represents 6% losses.
# Tx_Lifetime::Integer       # Transmission economic lifetime. Units: years
# Tx_Cost_Scaling::Float64   # Transmission cost scaling factor
# Tx_Exp_UB_Factor::Float64   # Factor for setting upper bound on size of transmission expansion as multiple of existing capacity
## Synchronous condensers parameters:
# SynCon_Parameters: Lifetime::Integer # Synchronous condenser plant economic lifetime. Units: years
# SynCon_Parameters: CapCost::Float64  # Cost of synchronous condenser type. Units: $/MW

# Expected Flags:
#   Flag H2::Bool  # :h2 - include H2 technology constraints if true.
#   Flag EV::Bool  # :ev - include EV technology constraints if true.
#   Flag HeatLoads::Bool # :heat - include Building Heat Load constraints if true.
#   Flag FixExistingCap::Bool  # FixExistingCapFlag - Fix the variables for existing capacity to the given values if true
#   Flag NewGasAllowed::Bool   # NoNewGas -  New gas technology built if true.
#   Flag NewDistillateAllowed::Bool   # NoNewDistillateAllowed - New liquid distillate technology built if true.

# Run Structure

mutable struct STABLErun
      RunInstance::STABLEinstance
      RunSettings::STABLErunSettings
      RunIdentifier::String
      RunTimestampStart::String  # run_timestamp
      RunTimestampFinish::String
      RunStatus::String
      RunNote::String    # A text annotation in order to aid identification of runs
      ResultsIdentifier::String  # Identifier for Results (e.g. filenames)
      ResultsPath::String        # The path to directory where Results will be written - resultsdir
end
# Construct from STrun::STABLErun the RunDescriptor::String # scenario_timestamp

# Examples of Tests:
# @assert STBL_scen.ScenarioName in Allowed_Scenarios
# @assert ScenarioName in Allowed_Scenarios


text:
#!/usr/bin/env julia

using Pkg
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--base", "-b"
        help = "Baseline commit id or branch"
        arg_type = String
        default = "master"
        "--target", "-t"
        help = "Target commit id or branch"
        arg_type = String
        required = true
        "--pr"
        help = "ID of the PR to comment on"
        arg_type = Int
        "--comment", "-c"
        help = "ID of the comment you want to update"
        arg_type = Int
        "--file", "-f"
        help = "Markdown filename"
        arg_type = String
    end

    return parse_args(s)
end

if isinteractive() == false
    Pkg.activate(".")
    Pkg.instantiate()

    args = parse_commandline()
    using PkgBenchmark
    using ConstraintSolver, Cbc
    using GitHub, JSON, Statistics

    github_auth = GitHub.authenticate(ENV["GITHUB_AUTH"])
    target = args["target"]
    base = args["base"]

    baseline_config = BenchmarkConfig(id = base, juliacmd = `julia -O3`)
    target_config = BenchmarkConfig(id = target, juliacmd = `julia -O3`)

    judged = judge("ConstraintSolver", target_config, baseline_config; f = median)

    markdown = sprint(export_markdown, judged)
    if args["file"] !== nothing
        export_markdown(args["file"], judged)
    end
    if args["comment"] !== nothing
        comment = edit_comment(
            "Wikunia/ConstraintSolver.jl",
            Comment(args["comment"]),
            :pr;
            params = Dict(:body => markdown),
            auth = github_auth,
        )
        println("Updated comment: $(comment.html_url)")
    elseif args["pr"] !== nothing
        comment = create_comment(
            "Wikunia/ConstraintSolver.jl",
            PullRequest(args["pr"]),
            markdown;
            auth = github_auth,
        )
        println("New comment: $(comment.html_url)")
    end
end


text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

## Powers

# CRlibm does not contain a correctly-rounded ^ function for Float64
# Use the BigFloat version from MPFR instead, which is correctly-rounded:

# Write explicitly like this to avoid ambiguity warnings:
for T in (:Integer, :Rational, :Float64, :BigFloat, :Interval)
    @eval ^(a::Interval{Float64}, x::$T) = atomic(Interval{Float64}, big53(a)^x)
end

# Integer power:

# overwrite new behaviour for small integer powers:
# ^{p}(x::IntervalArithmetic.Interval, ::Type{Val{p}}) = x^p

function ^(a::Interval{BigFloat}, n::Integer)
    isempty(a) && return a
    n == 0 && return one(a)
    n == 1 && return a
    # n == 2 && return sqr(a)
    n < 0 && a == zero(a) && return emptyinterval(a)

    T = BigFloat

    if isodd(n) # odd power
        isentire(a) && return a
        if n > 0
            a.lo == 0 && return @round(0, a.hi^n)
            a.hi == 0 && return @round(a.lo^n, 0)
            return @round(a.lo^n, a.hi^n)
        else
            if a.lo ≥ 0
                a.lo == 0 && return @round(a.hi^n, Inf)
                return @round(a.hi^n, a.lo^n)

            elseif a.hi ≤ 0
                a.hi == 0 && return @round(-Inf, a.lo^n)
                return @round(a.hi^n, a.lo^n)
            else
                return entireinterval(a)
            end
        end

    else # even power
        if n > 0
            if a.lo ≥ 0
                return @round(a.lo^n, a.hi^n)
            elseif a.hi ≤ 0
                return @round(a.hi^n, a.lo^n)
            else
                return @round(mig(a)^n, mag(a)^n)
            end

        else
            if a.lo ≥ 0
                return @round(a.hi^n, a.lo^n)
            elseif a.hi ≤ 0
                return @round(a.lo^n, a.hi^n)
            else
                return @round(mag(a)^n, mig(a)^n)
            end
        end
    end
end

function sqr(a::Interval{T}) where T<:Real
    return a^2
    # isempty(a) && return a
    # if a.lo ≥ zero(T)
    #     return @round(a.lo^2, a.hi^2)
    #
    # elseif a.hi ≤ zero(T)
    #     return @round(a.hi^2, a.lo^2)
    # end
    #
    # return @round(mig(a)^2, mag(a)^2)
end

^(a::Interval{BigFloat}, x::AbstractFloat) = a^big(x)

# Floating-point power of a BigFloat interval:
function ^(a::Interval{BigFloat}, x::BigFloat)

    domain = Interval{BigFloat}(0, Inf)

    if a == zero(a)
        a = a ∩ domain
        x > zero(x) && return zero(a)
        return emptyinterval(a)
    end

    isinteger(x) && return a^(round(Int, x))
    x == 0.5 && return sqrt(a)

    a = a ∩ domain
    (isempty(x) || isempty(a)) && return emptyinterval(a)

    xx = atomic(Interval{BigFloat}, x)

    # @round() can't be used directly, because both arguments may
    # Inf or -Inf, which throws an error
    # lo = @round(a.lo^xx.lo, a.lo^xx.lo)
    lolod = @round_down(a.lo^xx.lo)
    lolou = @round_up(a.lo^xx.lo)
    lo = (lolod == Inf || lolou == -Inf) ?
        wideinterval(lolod) : Interval(lolod, lolou)

    # lo1 = @round(a.lo^xx.hi, a.lo^xx.hi)
    lohid = @round_down(a.lo^xx.hi)
    lohiu = @round_up(a.lo^xx.hi)
    lo1 = (lohid == Inf || lohiu == -Inf) ?
        wideinterval(lohid) : Interval(lohid, lohiu)

    # hi = @round(a.hi^xx.lo, a.hi^xx.lo)
    hilod = @round_down(a.hi^xx.lo)
    hilou = @round_up(a.hi^xx.lo)
    hi = (hilod == Inf || hilou == -Inf) ?
        wideinterval(hilod) : Interval(hilod, hilou)

    # hi1 = @round(a.hi^xx.hi, a.hi^xx.hi)
    hihid = @round_down(a.hi^xx.hi)
    hihiu = @round_up(a.hi^xx.hi)
    hi1 = (hihid == Inf || hihiu == -Inf) ?
        wideinterval(hihid) : Interval(hihid, hihiu)

    lo = hull(lo, lo1)
    hi = hull(hi, hi1)

    return hull(lo, hi)
end

function ^(a::Interval{Rational{T}}, x::AbstractFloat) where T<:Integer
    a = Interval(a.lo.num/a.lo.den, a.hi.num/a.hi.den)
    a = a^x
    atomic(Interval{Rational{T}}, a)
end

# Rational power
function ^(a::Interval{BigFloat}, r::Rational{S}) where S<:Integer
    T = BigFloat
    domain = Interval{T}(0, Inf)

    if a == zero(a)
        a = a ∩ domain
        r > zero(r) && return zero(a)
        return emptyinterval(a)
    end

    isinteger(r) && return atomic(Interval{T}, a^round(S,r))
    r == one(S)//2 && return sqrt(a)

    a = a ∩ domain
    (isempty(r) || isempty(a)) && return emptyinterval(a)

    y = atomic(Interval{BigFloat}, r)

    a^y
end

# Interval power of an interval:
function ^(a::Interval{BigFloat}, x::Interval)
    T = BigFloat
    domain = Interval{T}(0, Inf)

    a = a ∩ domain

    (isempty(x) || isempty(a)) && return emptyinterval(a)

    lo1 = a^x.lo
    lo2 = a^x.hi
    lo1 = hull(lo1, lo2)

    hi1 = a^x.lo
    hi2 = a^x.hi
    hi1 = hull(hi1, hi2)

    hull(lo1, hi1)
end

function sqrt(a::Interval{T}) where T
    domain = Interval{T}(0, Inf)
    a = a ∩ domain

    isempty(a) && return a

    @round(sqrt(a.lo), sqrt(a.hi))  # `sqrt` is correctly-rounded
end

"""
    pow(x::Interval, n::Integer)

A faster implementation of `x^n`, currently using `power_by_squaring`.
`pow(x, n)` will usually return an interval that is slightly larger than that calculated by `x^n`, but is guaranteed to be a correct
enclosure when using multiplication with correct rounding.
"""
function pow(x::Interval, n::Integer)  # fast integer power

    isempty(x) && return x

    if iseven(n) && 0 ∈ x

        return hull(zero(x),
                    hull(Base.power_by_squaring(Interval(mig(x)), n),
                        Base.power_by_squaring(Interval(mag(x)), n))
            )

    else

      return hull( Base.power_by_squaring(Interval(x.lo), n),
                    Base.power_by_squaring(Interval(x.hi), n) )

    end

end

function pow(x::Interval, y::Real)  # fast real power, including for y an Interval

    isempty(x) && return x

    return exp(y * log(x))

end

for f in (:exp, :expm1)
    @eval begin
        function ($f)(a::Interval{T}) where T
            isempty(a) && return a
            @round( ($f)(a.lo), ($f)(a.hi) )
        end
    end
end

for f in (:exp2, :exp10)

    @eval function ($f)(x::BigFloat, r::RoundingMode)  # add BigFloat functions with rounding:
            setrounding(BigFloat, r) do
                ($f)(x)
            end
        end

    @eval ($f)(a::Interval{Float64}) = atomic(Interval{Float64}, $f(big53(a)))  # no CRlibm version

    @eval function ($f)(a::Interval{BigFloat})
            isempty(a) && return a
            @round( ($f)(a.lo), ($f)(a.hi) )
        end
end

for f in (:log, :log2, :log10, :log1p)

    @eval function ($f)(a::Interval{T}) where T
            domain = Interval{T}(0, Inf)
            a = a ∩ domain

            (isempty(a) || a.hi ≤ zero(T)) && return emptyinterval(a)

            @round( ($f)(a.lo), ($f)(a.hi) )

        end
end


text:
using ADCME
using LinearAlgebra
using Random

npoints = parse(Int64, ARGS[1])
SEED = parse(Int64, ARGS[2])

Random.seed!(SEED)

x0 = LinRange(0, 1, 100)|>collect
y0 = sin.(π*x0)

θ = Variable(fc_init([1,20,20,20,1]))
y = fc(x0, [20,20,20,1], θ)

idx = Int64.(round.(collect(LinRange(1, 100, npoints+2))))[2:end-1]

loss = sum((y[idx] - y0[idx])^2)

H = hessian(loss, θ)

sess = Session(); init(sess)
BFGS!(sess, loss)
H0 = run(sess, H)
λ = real.(eigvals(H0))
nlambda = length(findall(λ .> λ[end] * 1e-6))

db = Database("hessian.db")
execute(db, """
CREATE TABLE IF NOT EXISTS eigvals (
    seed integer,
    npoints integer,
    nlambda integer,
    PRIMARY KEY (seed, npoints)
)""")
execute(db, """
INSERT OR REPLACE INTO eigvals VALUES ($SEED, $npoints, $nlambda)
""")
close(db)


text:
# Run this script from the command line:
# $ julia sekira ARGS...

# To see available and required ARGS, run
# $ julia sekira --help
# or see
# julia> ?Sekira.plateau
using Sekira
add = BoseFS(ntuple(i -> ifelse(i==10, 10, 0), 20))
H = HubbardMom1D(add; u=0.5)

Sekira.plateau(H, ARGS; post_step=ProjectedEnergy(H, DVec(add => 1.0)))


text:
#### 03c - BART predictions ####
include("required.jl")

## Conditional arguments
# save_data = true

## Load predictions from CSV
summaries = CSV.read(joinpath("data", "proc", "bart_summaries.csv"), DataFrame)
varimps = CSV.read(joinpath("data", "proc", "bart_varimps.csv"), DataFrame)
pred_df = CSV.read(
    joinpath("data", "proc", "bart_predictions_prob.csv"), DataFrame; missingstrings=["NA"]
)
lower_df = CSV.read(
    joinpath("data", "proc", "bart_predictions_lower.csv"), DataFrame; missingstrings=["NA"]
)
upper_df = CSV.read(
    joinpath("data", "proc", "bart_predictions_upper.csv"), DataFrame; missingstrings=["NA"]
)
pres_df = CSV.read(
    joinpath("data", "proc", "bart_predictions_pres.csv"), DataFrame; missingstrings=["NA"]
)

## Create Y matrices

# Get matrix Y
Y = replace(Array(pres_df), missing => nothing) |> Array{Union{Nothing,Float32}}
Yprob = replace(Array(pred_df), missing => nothing) |> Array{Union{Nothing,Float32}}
Ylower = replace(Array(lower_df), missing => nothing) |> Array{Union{Nothing,Float32}}
Yupper = replace(Array(upper_df), missing => nothing) |> Array{Union{Nothing,Float32}}

# Set values to nothing if no species present
inds_zeros = _indsnotobs(Y)
Y[inds_zeros, :] .= nothing

## Create distributions

# Load raw distributions (for grid size)
glossary = CSV.read(joinpath("data", "proc", "glossary.csv"), DataFrame)
spenames = filter(:type => ==("species"), glossary).full_name
distributions = [
    geotiff(SimpleSDMPredictor, joinpath("data", "raster", "distributions_raw.tif"), i) for
    i in eachindex(spenames)
]
raw_distributions = copy(distributions)

# Cut to Quebec coordinates (optional)
# coords_qc = (left = -80.0, right = -55.0, bottom = 45.0, top = 63.0)
# raw_distributions = [d[coords_qc] for d in raw_distributions]

# Get layer dimensions & limits
dims = size(raw_distributions[1])
lims = boundingbox(raw_distributions[1])

# Create distribution layers
layers = []
for Y in (Y, Yprob, Ylower, Yupper)
    Ydistrib = replace(Y, 0.0 => nothing)
    Ygrids = [Ydistrib[:, col] for col in 1:size(Ydistrib, 2)]
    Ygrids = reshape.(Ygrids, dims...) .|> Array
    local distributions = SimpleSDMResponse.(Ygrids, lims...)
    push!(layers, distributions)
end
distributions, prob_distrib, lower_distrib, upper_distrib = layers;
distributions

## Export results
# save_data = true
if (@isdefined save_data) && save_data == true
    # Export BART distributions
    geotiff(joinpath("data", "raster", "distributions_bart.tif"), distributions)

    # Extras
    geotiff(joinpath("data", "raster", "bart_xtras_prob-distrib.tif"), prob_distrib)
    geotiff(joinpath("data", "raster", "bart_xtras_lower-distrib.tif"), lower_distrib)
    geotiff(joinpath("data", "raster", "bart_xtras_upper-distrib.tif"), upper_distrib)

    # Update placeholder files (as files are too big for version control)
    placeholder_paths = [
        joinpath("data", "raster", "bart_xtras_prob-distrib_placeholder.tif"),
        joinpath("data", "raster", "bart_xtras_lower-distrib_placeholder.tif"),
        joinpath("data", "raster", "bart_xtras_upper-distrib_placeholder.tif"),
    ]
    for p in placeholder_paths
        open(p, "w") do io
            write(io, string(Dates.now()))
        end
    end
end

## Map uncertainty

# Get uncertainty per species
uncertainty = upper_distrib .- lower_distrib
plot_layer(uncertainty[1]; c=:viridis)

# Uncertainty sum
uncertainty_sum = sum(uncertainty)
plot_layer(uncertainty_sum; c=:viridis)

# Uncertainty mean
uncertainty_mean = mean(uncertainty)
plot_layer(uncertainty_mean; c=:viridis)
histogram(uncertainty_mean)

# Plot uncertainty & histogram
function uncertainty_plot(layer; title="", kw...)
    unc_map = plot_layer(
        layer;
        c=:viridis,
        clim=extrema(layer),
        title="Uncertainty",
        colorbar_title="Uncertainty mean",
    )
    unc_hist = histogram(
        layer;
        legend=:none,
        # ylim = maxrange, # xlabel = "Difference",
        title="Distribution of uncertainty values",
        orientation=:horizontal,
    )
    unc_title = plot(; annotation=(0.5, 0.5, "$(title)"), framestyle=:none)
    l = @layout [t{0.01h}; a{0.6w} b{0.38w}]
    unc_plot = plot(unc_title, unc_map, unc_hist; size=(800, 400), layout=l)
    return unc_plot
end
unc_plot = uncertainty_plot(uncertainty_mean; title="Uncertainty mean (BART SDMs)", dpi=200)

# Export uncertainty plot
if (@isdefined save_figures) && save_figures == true
    savefig(unc_plot, joinpath("fig", "bart", "x_bart_uncertainty.png"))
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
module ZXCalculus

using YaoHIR
using MLStyle
using YaoLocations
using YaoHIR: X, Y, Z, H, S, T, Rx, Ry, Rz, UGate, shift
using YaoLocations: plain

export convert_to_chain, convert_to_zxd

include("phase.jl")
include("scalar.jl")
include("abstract_zx_diagram.jl")
include("zx_layout.jl")
include("zx_diagram.jl")
include("zx_graph.jl")
include("rules.jl")
include("simplify.jl")

include("ir.jl")

include("circuit_extraction.jl")
include("phase_teleportation.jl")

include("deprecations.jl")

end # module


text:
# For low level cudnn functions that require a pointer to a number
cptr(x,a::CuArray{Float64})=Float64[x]
cptr(x,a::CuArray{Float32})=Float32[x]
cptr(x,a::CuArray{Float16})=Float32[x]

# Conversion between Julia and CUDNN datatypes
cudnnDataType(::Type{Float16})=CUDNN_DATA_HALF
cudnnDataType(::Type{Float32})=CUDNN_DATA_FLOAT
cudnnDataType(::Type{Float64})=CUDNN_DATA_DOUBLE
juliaDataType(a)=(a==CUDNN_DATA_HALF ? Float16 :
                  a==CUDNN_DATA_FLOAT ? Float32 :
                  a==CUDNN_DATA_DOUBLE ? Float64 : error())

tuple_strides(A::Tuple) = _strides((1,), A)
_strides(out::Tuple{Int}, A::Tuple{}) = ()
_strides(out::NTuple{N,Int}, A::NTuple{N}) where {N} = out
function _strides(out::NTuple{M,Int}, A::Tuple) where M
    Base.@_inline_meta
    _strides((out..., out[M]*A[M]), A)
end

# Descriptors

mutable struct TensorDesc; ptr; end
free(td::TensorDesc) = cudnnDestroyTensorDescriptor(td.ptr)
Base.unsafe_convert(::Type{cudnnTensorDescriptor_t}, td::TensorDesc) = td.ptr
Base.unsafe_convert(::Type{Ptr{Nothing}}, td::TensorDesc) = convert(Ptr{Nothing}, td.ptr)

function TensorDesc(T::Type, size::NTuple{N,Integer}, strides::NTuple{N,Integer} = tuple_strides(size)) where N
    sz = Cint.(size) |> reverse |> collect
    st = Cint.(strides) |> reverse |> collect
    d = Ref{cudnnTensorDescriptor_t}()
    cudnnCreateTensorDescriptor(d)
    cudnnSetTensorNdDescriptor(d[], cudnnDataType(T), length(sz), sz, st)
    this = TensorDesc(d[])
    finalizer(free, this)
    return this
end

TensorDesc(a::CuArray) = TensorDesc(eltype(a), size(a), strides(a))

mutable struct FilterDesc
  ptr
end
free(fd::FilterDesc)=cudnnDestroyFilterDescriptor(fd.ptr)
Base.unsafe_convert(::Type{cudnnFilterDescriptor_t}, fd::FilterDesc)=fd.ptr
Base.unsafe_convert(::Type{Ptr{Nothing}}, fd::FilterDesc)=fd.ptr

function createFilterDesc()
  d = Ref{cudnnFilterDescriptor_t}()
  @check cudnnCreateFilterDescriptor(d)
  return d[]
end

function FilterDesc(T::Type, size::Tuple; format = CUDNN_TENSOR_NCHW)
    # The only difference of a FilterDescriptor is no strides.
    sz = Cint.(size) |> reverse |> collect
    d = createFilterDesc()
    version() >= v"5" ?
        cudnnSetFilterNdDescriptor(d, cudnnDataType(T), format, length(sz), sz) :
    version() >= v"4" ?
        cudnnSetFilterNdDescriptor_v4(d, cudnnDataType(T), format, length(sz), sz) :
        cudnnSetFilterNdDescriptor(d, cudnnDataType(T), length(sz), sz)
    this = FilterDesc(d)
    finalizer(free, this)
    return this
end

FilterDesc(a::CuArray; format = CUDNN_TENSOR_NCHW) = FilterDesc(eltype(a), size(a), format = format)

function Base.size(f::FilterDesc)
  typ = Ref{Cuint}()
  format = Ref{Cuint}()
  ndims = Ref{Cint}()
  dims = Vector{Cint}(undef, 8)
  @check ccall((:cudnnGetFilterNdDescriptor,libcudnn), cudnnStatus_t, (Ptr{Nothing}, Cint, Ptr{UInt32}, Ptr{UInt32}, Ptr{Cint}, Ptr{Cint}),
    f, 8, typ, format, ndims, dims)
  @assert ndims[] ≤ 8
  return (dims[1:ndims[]]...,) |> reverse
end

mutable struct ConvDesc; ptr; end
free(cd::ConvDesc) = cudnnDestroyConvolutionDescriptor(cd.ptr)
Base.unsafe_convert(::Type{cudnnConvolutionDescriptor_t}, cd::ConvDesc)=cd.ptr

function cdsize(w, nd)
    isa(w, Integer) ? Cint[fill(w,nd)...] :
    length(w)!=nd ? error("Dimension mismatch") :
    Cint[reverse(w)...]
end

pdsize(w, nd)=Cint[reverse(psize(w,nd))...]
psize(w, nd)=(isa(w,Integer)  ? fill(w,nd) : length(w) != nd ? error("Dimension mismatch") : w)

function ConvDesc(T, N, padding, stride, dilation, mode)
    cd = Ref{cudnnConvolutionDescriptor_t}()
    cudnnCreateConvolutionDescriptor(cd)
    version() >= v"4" ? cudnnSetConvolutionNdDescriptor(cd[],N,cdsize(padding,N),cdsize(stride,N),cdsize(dilation,N),mode,cudnnDataType(T)) :
    version() >= v"3" ? cudnnSetConvolutionNdDescriptor_v3(cd[],N,cdsize(padding,N),cdsize(stride,N),cdsize(dilation,N),mode,cudnnDataType(T)) :
    cudnnSetConvolutionNdDescriptor(cd[],N,cdsize(padding,N),cdsize(stride,N),cdsize(dilation,N),mode)
    this = ConvDesc(cd[])
    finalizer(free, this)
    return this
end

mutable struct PoolDesc; ptr; end
free(pd::PoolDesc)=cudnnDestroyPoolingDescriptor(pd.ptr)
Base.unsafe_convert(::Type{cudnnPoolingDescriptor_t}, pd::PoolDesc)=pd.ptr

function PoolDesc(nd, window, padding, stride, mode, maxpoolingNanOpt=CUDNN_NOT_PROPAGATE_NAN)
    pd = Ref{cudnnPoolingDescriptor_t}()
    cudnnCreatePoolingDescriptor(pd)
    cudnnSetPoolingNdDescriptor(pd[],mode,maxpoolingNanOpt,nd,pdsize(window,nd),pdsize(padding,nd),pdsize(stride,nd))
    this = PoolDesc(pd[])
    finalizer(free, this)
    return this
end

mutable struct ActivationDesc; ptr; end
free(ad::ActivationDesc)=cudnnDestroyActivationDescriptor(ad.ptr)
Base.unsafe_convert(::Type{cudnnActivationDescriptor_t}, ad::ActivationDesc)=ad.ptr

function ActivationDesc(mode, coeff, reluNanOpt=CUDNN_NOT_PROPAGATE_NAN)
    ad = Ref{cudnnActivationDescriptor_t}()
    cudnnCreateActivationDescriptor(ad)
    cudnnSetActivationDescriptor(ad[],mode,reluNanOpt,coeff)
    this = ActivationDesc(ad[])
    finalizer(free, this)
    return this
end


text:
# When the previous transformation didnt modify the table inplace
# then applying f to rows would ignore the previous transformation.
# One solution is to create a deepcopy of the table, apply the
# previsous transformation inplace on that copy and finally apply f
# to its rows. Copying the whole table is wasteful, so ByRow and
# ByRow! structs can store the columns needed.
for T in (:ByRow, :ByRow!)
    F = Symbol(lowercase(String(T)))
    @eval begin
        struct $(T){Bool,S} <: AbstractTransformation{Bool} where S
            f::S
            cols::Tuple{Vararg{Symbol}}
        end
        $(F)(x::Callable) = $(T){false, typeof(x)}(x, Tuple{Vararg{Symbol,0}}())
        $(F)(x::Callable, cols) = $(T){false, typeof(x)}(x, cols)
        $(F)(x::Pair{Tuple{Vararg{Symbol}},Callable}) = $(T){false, typeof(x)}(last(x), first(x))
        $(F)(x::Pair{Symbol,Callable}) = $(T){false, typeof(x)}(last(x), (first(x),))

        @inline function Broadcast.broadcasted(::typeof($(F)), x::Callable)
            $(T){true, typeof(x)}(x, Tuple{Vararg{Symbol,0}}())
        end
        @inline function Broadcast.broadcasted(::typeof($(F)), x::Callable, cols)
            $(T){true, typeof(x)}(x, cols)
        end
        @inline function Broadcast.broadcasted(::typeof($(F)), x::Pair{Tuple{Vararg{Symbol}},Callable})
            $(T){true, typeof(x)}(last(x), first(x))
        end
        @inline function Broadcast.broadcasted(::typeof($(F)), x::Pair{Symbol,Callable})
            $(T){true, typeof(x)}(last(x), (first(x),))
        end
        @inline (t::$(T){true})(x, colname) = broadcasted(t.f, x, colname)
        @inline (t::$(T){false})(x, colname) = t.f(x, colname)
        Base.show(io::IO, f::$(T){false}) = print(io, $(F), "(", f.f, ")")
        Base.show(io::IO, f::$(T){true}) = print(io, $(F), ".(", f.f, ")")
    end
end

for T in (:ByTab, :ByTab!)
    F = Symbol(lowercase(String(T)))
    @eval begin
        struct $(T){Bool,S} <: AbstractTransformation{Bool} where S
            f::S
            cols::Tuple{Vararg{Symbol}}
        end
        $(F)(x::Callable) = $(T){false, typeof(x)}(x, Tuple{Vararg{Symbol,0}}())
        $(F)(x::Callable, cols) = $(T){false, typeof(x)}(x, cols)
        $(F)(x::Pair{Tuple{Vararg{Symbol}},Callable}) = $(T){false, typeof(x)}(last(x), first(x))
        $(F)(x::Pair{Symbol,Callable}) = $(T){false, typeof(x)}(last(x), (first(x),))

        @inline (t::$(T){false})(x, colname) = t.f(x, colname)
        Base.show(io::IO, f::$(T){false}) = print(io, $(F), "(", f.f, ")")
    end
end

for T in (:ByCol, :ByCol!)
    F = Symbol(lowercase(String(T)))
    @eval begin
        struct $(T){Bool,S} <: AbstractTransformation{Bool} where S
            f::S
        end
        $(F)(x::Callable) = $(T){false, typeof(x)}(x)
        @inline function Broadcast.broadcasted(::typeof($(F)), x::Callable)
            $(T){true, typeof(x)}(x)
        end
        @inline (t::$(T){true})(x) = broadcasted(t.f, x)
        @inline (t::$(T){false})(x) = t.f(x)
        Base.show(io::IO, f::$(T){false}) = print(io, $(F), "(", f.f, ")")
        Base.show(io::IO, f::$(T){true}) = print(io, $(F), ".(", f.f, ")")
    end
end

const RowwiseTrans{T} = Union{ByRow{T}, ByRow!{T}} where T
const ColwiseTrans{T} = Union{ByCol{T}, ByCol!{T}} where T
const TabwiseTrans{T} = Union{ByTab{T}, ByTab!{T}} where T
const InplaceTrans{T} = Union{ByCol!{T}, ByRow!{T}, ByTab!{T}} where T
const OnCopyTrans{T} = Union{ByCol{T}, ByRow{T}, ByTab{T}} where T

####################################################################################################
# Nested ###########################################################################################
####################################################################################################

apply_trans_nested(::Nothing, colname, tab) = columntable(tab)[colname]

function apply_trans_nested(c::Composition, colname, tab)
    apply_trans_nested(c.funcs, colname, tab, columntable(tab)[colname], nothing)
end

function apply_trans_nested(fs::Tuple, colname, tab, x, f)
    apply_trans_nested(
        Base.tail(fs),
        colname,
        tab,
        _apply_trans_nested(first(fs), f, x, tab, colname),
        first(fs))
end

apply_trans_nested(::Tuple{}, colname, tab, x, f) = x

apply_trans_nested(::Tuple{}, colname, tab, x, f::OnCopyTrans{true}) = materialize(x)

function apply_trans_nested(::Tuple{}, colname, tab, x, f::InplaceTrans{true})
    materialize!(columntable(tab)[colname], x)
end

function apply_trans_nested(::Tuple{}, colname, tab, x, f::InplaceTrans{false})
    copyto!(columntable(tab)[colname], x)
end

# RowwiseTrans #####################################################################################

function _apply_trans_nested(f::RowwiseTrans, prev_f::OnCopyTrans{true}, x, tab, colname)
    if isempty(f.cols)
        _df = deepcopy(tab)
    else
        _df = deepcopy(subset_cols(tab, f.cols))
    end
    materialize!(columntable(_df)[colname], x)
    f(rowtable(_df), colname)
end

function _apply_trans_nested(
        f::RowwiseTrans,
        prev_f::AbstractTransformation{false},
        x,
        tab,
        colname)
    if isempty(f.cols)
        _df = deepcopy(tab)
    else
        _df = deepcopy(subset_cols(tab, f.cols))
    end
    copyto!(columntable(_df)[colname], x)
    f(rowtable(_df), colname)
end

function _apply_trans_nested(f::RowwiseTrans, prev_f::InplaceTrans{true}, x, tab, colname)
    materialize!(columntable(tab)[colname], x)
    f(rowtable(tab), colname)
end

_apply_trans_nested(f::RowwiseTrans, prev_f::Nothing, x, tab, colname) = f(rowtable(tab), colname)

_apply_trans_nested(f::T, prev_f::T, x, tab, colname) where T<:RowwiseTrans{true} = f(rowtable(tab), colname)

# ColwiseTrans #####################################################################################

function _apply_trans_nested(
        f::ColwiseTrans,
        prev_f::Union{AbstractTransformation{false}, Nothing},
        x,
        tab,
        colname)
    f(x)
end

function _apply_trans_nested(f::ColwiseTrans, prev_f::InplaceTrans{true}, x, tab, colname)
    f(materialize!(columntable(tab)[colname], x))
end

function _apply_trans_nested(f::ColwiseTrans, prev_f::OnCopyTrans{true}, x, tab, colname)
    f(materialize(x))
end

_apply_trans_nested(f::T, prev_f::T, x, tab, colname) where T<:ColwiseTrans{true} = f(x)

# TabwiseTrans #####################################################################################

_apply_trans_nested(f::TabwiseTrans, prev_f::Nothing, x, tab, colname) = f(columntable(tab), colname)

function _apply_trans_nested(f::TabwiseTrans, prev_f::AbstractTransformation{false}, x, tab, colname)
    if isempty(f.cols)
        _df = columntable(deepcopy(tab))
    else
        _df = subset_cols(columntable(deepcopy(tab)), f.cols)
    end
    copyto!(_df[colname], x)
    f(_df, colname)
end

function _apply_trans_nested(f::TabwiseTrans, prev_f::InplaceTrans{true}, x, tab, colname)
    _df = columntable(tab)
    materialize!(_df[colname], x)
    f(_df, colname)
end

function _apply_trans_nested(f::TabwiseTrans, prev_f::OnCopyTrans{true}, x, tab, colname)
    if isempty(f.cols)
        _df = columntable(deepcopy(tab))
    else
        _df = subset_cols(columntable(deepcopy(tab)), f.cols)
    end
    materialize!(_df[colname], x)
    f(_df, colname)
end

####################################################################################################
# Flat #############################################################################################
####################################################################################################
# TODO: inplace byrow and bytab

function apply_trans_flat(f::ByRow, selected_colnames, newcolname, nt)
    _selected_colnames = ensure_tuple(length(f.cols) > 1 ? f.cols : selected_colnames)
    materialize(f(rowtable(NamedTuple{_selected_colnames}(nt)), newcolname))
end

function apply_trans_flat(f::ByRow!{false}, selected_colnames, newcolname, nt)
    !(newcolname in colnames(nt)) && throw(KeyError(newcolname))
    _selected_colnames = ensure_tuple(length(f.cols) > 1 ? f.cols : selected_colnames)
    copyto!(nt[newcolname], f(rowtable(NamedTuple{_selected_colnames}(nt)), newcolname))
end

function apply_trans_flat(f::ByRow!{true}, selected_colnames, newcolname, nt)
    !(newcolname in colnames(nt)) && throw(KeyError(newcolname))
    _selected_colnames = ensure_tuple(length(f.cols) > 1 ? f.cols : selected_colnames)
    materialize!(nt[newcolname], f(rowtable(NamedTuple{_selected_colnames}(nt)), newcolname))
end

function apply_trans_flat(f::ByTab{false}, selected_colnames, newcolname, nt)
    _selected_colnames = ensure_tuple(length(f.cols) > 1 ? f.cols : selected_colnames)
    materialize(f(NamedTuple{_selected_colnames}(nt), newcolname))
end

function apply_trans_flat(f::ByTab!{false}, selected_colnames, newcolname, nt)
    !(newcolname in colnames(nt)) && throw(KeyError(newcolname))
    _selected_colnames = ensure_tuple(length(f.cols) > 1 ? f.cols : selected_colnames)
    copyto!(nt[newcolname], f(NamedTuple{_selected_colnames}(nt), newcolname))
end

function apply_trans_flat(f::ByCol, selected_colnames, newcolname, nt)
    if length(selected_colnames) > 1
        throw(error("Multiple columns as an input for $(typeof(f)) is not implemented."))
    else
        materialize(f(nt[@inbounds(selected_colnames[1])]))
    end
end

function apply_trans_flat(f::ByCol!{true}, selected_colnames, newcolname, nt)
    if !(newcolname in colnames(nt))
        throw(KeyError(newcolname))
    elseif length(selected_colnames) > 1
        throw(error("Multiple columns as an input for $(typeof(f)) is not implemented."))
    else
        x = nt[selected_colnames[1]]
        materialize!(x, f(x))
    end
end

function apply_trans_flat(f::ByCol!{false}, selected_colnames, newcolname, nt)
    if !(newcolname in colnames(nt))
        throw(KeyError(newcolname))
    elseif length(selected_colnames) > 1
        throw(error("Multiple columns as an input for $(typeof(f)) is not implemented."))
    else
        x = nt[selected_colnames[1]]
        copyto!(x, f(x))
    end
end

transformation(f::Callable) = bycol(f)
transformation(t::AbstractTransformation) = t

"""
```
bycol(f::Callable)           # ~ f(column)
bycol!(f::Callable)          # ~ f(column)
byrow(f::Callable[, cols])   # ~ f(rowtable, name::Symbol))
byrow!(f::Callable[, cols])  # ~ f(rowtable, name::Symbol))
bytab(f::Callable[, cols])   # ~ f(coltable, name::Symbol))
bytab!(f::Callable[, cols])  # ~ f(coltable, name::Symbol))

bycol.(f::Callable)           # ~ f(element)
bycol!.(f::Callable)          # ~ f(element)
byrow.(f::Callable[, cols])   # ~ f(row, name::Symbol))
byrow!.(f::Callable[, cols])  # ~ f(row, name::Symbol))
```

Wrappers around transformations `f` to be applied to selected columns in `select(tab, args...; kwargs...)`.
`cols` are the source columns to make available inside the function.

As their name suggest, functions wrapped in `bycol[!.]` are applied to individual selected columns,
`byrow[!.]` to rows and `bytab[!]` to a whole table (a `NamedTuple` of columns).
`byrow[!.]` and `bytab[!]` will supply the name of the currently selected column as the second argument to `f`

Functions ending with `!` apply the transformation inplace. This requires the columns to exist
(i.e. you won't be able to create a new column in `kwargs...` this way) and to have the same input and output eltype.

Broadcasting of these wrappers is a signal to interpret the wrapped function `f` as an elementwise function.
You you chain multiple broadcasted transformation together and if they are of same type (e.g. they are all `bycol.`)
then the multiple wrapped function would be fused together.

See also: [`select`](@ref), [`selection`](@ref), [`renaming`](@ref)
"""
bycol, bycol!, byrow, byrow!, bytab, bytab!, transformation


text:
#### Frank wolfe for matrix completion
using FirstOrderOptimization
import Base: axpy!, scale!

## Generate data

# make it deterministic
srand(10)

m, n, r, rtrue = 100, 100, 2, 2
Xtrue = randn(m, rtrue) * randn(rtrue, n)
tau = nucnorm(Xtrue)

pobs = .5
obs = sprand(m, n, pobs).>0
iobs = find(obs)
A = IndexingOperator(m, n, iobs)
b = A*Xtrue

nonsketched = true

params = FrankWolfeParams(maxiters = 200, reltol = 1e-1, stepsize = BacktrackingStepSize(1,.5,.3)) # DecreasingStepSize(2,1))

srand(10)

objective_sketched(z) = .5*sum((z - b).^2)
grad_objective_sketched(z) = LowRankOperator(A, z-b, transpose = Symbol[:T, :N]) # concretely, should be A'*(z-b)
function min_lin_st_nucnorm_sketched(G, tau)
	u,s,v = svds(Array(G), nsv=1) # for this case, g is a sparse matrix so representing it is O(m)
	return LowRankOperator(-tau*u, v')
end
const_nucnorm(z) = tau # we'll always saturate the constraint, don't bother computing it
# I can't think of a pretty way to compute
# <g, tilde_x - x> = <A'(z-b), u tau v' - x> = <z-b, A(u tau v') - A(x)> = <z-b, A(u tau v') - z>
# dot_g_w_tx_minus_x(g, tilde_x, x) = dot(g.factors[2], g.factors[1]*tilde_x - x)
dot_g_w_tx_minus_x(g, Delta, z) = dot(vec(g.factors[2]), g.factors[1]*Delta - z)
# we end up computing A*Delta twice with this scheme
update_var!(z, Delta, a) = (scale!(z, 1-a); axpy!(a, A*Delta, z))

z = zeros(length(b))
X_sketched = frank_wolfe_sketched(
	z,
	objective_sketched, grad_objective_sketched,
	const_nucnorm,
	tau,
	min_lin_st_nucnorm_sketched,
	AsymmetricSketch(m,n,r),
	params,
#	LB = 0,
	verbose=true
	)

@show vecnorm(Array(X_sketched) - Xtrue) / vecnorm(Xtrue)

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
using IterTools
const gb = IterTools.groupby

function sumby_itgb{T,S}(by::AbstractArray{T}, v1::AbstractArray{S})
  res = Dict{T, S}()
  szero = zero(S)
  for i in IterTools.groupby(x -> x[1], zip(by,v1))
    byi = i[1][1]
    vali = mapreduce(x->x[2],+,i)
    index = ht_keyindex(res, byi)
    if index > 0
      @inbounds vw = res.vals[index]
      new_vw = vw + vali
      @inbounds res.vals[index] = new_vw
    else
      @inbounds res[byi] = vali
    end

  end
  return res
end

@time sumby_itgb(id4, v1)
@time sumby_itgb(id4, v1)


text:
using Random

using WordTokenizers

"""
    recursive_readdir(path::AbstractString="./")

recursive read all file from a dir. return a list of filenames.
"""
function recursive_readdir(path::AbstractString="./")
  ret = String[]
  for (root, dirs, files) in walkdir(path)
    append!(ret, map(file->joinpath(root, file), files))
  end
  ret
end

"""
    bert_pretrain_task(datachn::Channel, wordpiece::WordPiece;
                       buffer_size = 100,
                       channel_size = 100
                       wordpiece::WordPiece,
                       sentences_pool = sentences;
                       start_token = "[CLS]",
                       sep_token = "[SEP]",
                       mask_token = "[MASK]",
                       mask_ratio = 0.15,
                       real_token_ratio = 0.1,
                       random_token_ratio = 0.1,
                       whole_word_mask = false,
                       next_sentence_ratio = 0.5,
                       next_sentence = true,
                       tokenizer = tokenize,
                       istokenized = false,
                       return_real_sentence = false)

helper function to generate bert mask language modeling and next sentence prediction data. `datachn` is a `Channel` with input documents line by line.
"""
function bert_pretrain_task(datachn::Channel,
                       wordpiece::WordPiece;
                       buffer_size = 100,
                       channel_size = 100,
                       kwargs...
                       )
  outchn = Channel(channel_size)
  bert_pretrain_task(outchn, datachn, wordpiece; buffer_size = buffer_size, kwargs...)
  outchn
end

function bert_pretrain_task(outchn::Channel,
                       datachn::Channel,
                       wordpiece::WordPiece;
                       buffer_size = 100,
                       kwargs...
                       )
  task = @async begin
    buffer = Vector(undef, buffer_size)
    while isopen(datachn)
      i = 1
      eod = false
      while i <= buffer_size
        try
          sentence = take!(datachn)
          if isempty(sentence)
            continue
          else
            buffer[i] = sentence
            i+=1
          end
        catch e
          if isa(e, InvalidStateException) && e.state==:closed
            eod = true
            break
          else
            rethrow()
          end
        end
      end

      i -= 1

      if eod || i == buffer_size
        bert_pretrain_task(outchn, @view(buffer[1:(eod ? i - 1 : i)]), wordpiece; kwargs...)
      end
    end
  end
  bind(outchn, task)
end

function bert_pretrain_task(sentences,
                       wordpiece::WordPiece,
                       sentences_pool = sentences;
                       channel_size = 100,
                       kwargs...
                       )
  chn = Channel(channel_size)
  task = @async bert_pretrain_task(chn, sentences, wordpiece, sentences_pool; kwargs...)
  bind(chn, task)
  chn
end

function bert_pretrain_task(chn::Channel,
                       sentences,
                       wordpiece::WordPiece,
                       sentences_pool = sentences;
                       start_token = "[CLS]",
                       sep_token = "[SEP]",
                       mask_token = "[MASK]",
                       mask_ratio = 0.15,
                       real_token_ratio = 0.1,
                       random_token_ratio = 0.1,
                       whole_word_mask = false,
                       next_sentence_ratio = 0.5,
                       next_sentence = true,
                       tokenizer = tokenize,
                       istokenized = false,
                       return_real_sentence = false)

  foreach(enumerate(sentences)) do (i, sentence)
    sentenceA = masksentence(
      istokenized ? sentence : tokenizer(sentence),
      wordpiece;
      mask_token = mask_token,
      mask_ratio = mask_ratio,
      real_token_ratio = real_token_ratio,
      random_token_ratio = random_token_ratio,
      whole_word_mask = whole_word_mask
    )

    if next_sentence
      if rand() <= next_sentence_ratio && i != length(sentences)
        sentenceB = sentences[i+1]
        isnext = true
      else
        sentenceB = rand(sentences_pool)
        isnext = false
      end

      sentenceB = masksentence(
        istokenized ? sentenceB : tokenizer(sentenceB),
        wordpiece;
        mask_token = mask_token,
        mask_ratio = mask_ratio,
        real_token_ratio = real_token_ratio,
        random_token_ratio = random_token_ratio,
        whole_word_mask = whole_word_mask
      )

      masked_sentence = _wrap_sentence(sentenceA[1],
                                       sentenceB[1];
                                       start_token = start_token,
                                       sep_token = sep_token)

      sentence = _wrap_sentence(sentenceA[2],
                                sentenceB[2];
                                start_token = start_token,
                                sep_token = sep_token)

      mask_idx = _wrap_idx(sentenceA[3],
                           sentenceB[3],
                           length(sentenceA[1]))
    else
      masked_sentence = _wrap_sentence(sentenceA[1];
                                       start_token = start_token,
                                       sep_token = sep_token)

      sentence = _wrap_sentence(sentenceA[2];
                                start_token = start_token,
                                sep_token = sep_token)

      mask_idx = _wrap_idx(sentenceA[3])
    end

    masked_token = sentence[mask_idx]

    if return_real_sentence
      if next_sentence
        put!(chn, (masked_sentence, mask_idx, masked_token, isnext, sentence))
      else
        put!(chn, (masked_sentence, mask_idx, masked_token, sentence))
      end
    else
      if next_sentence
        put!(chn, (masked_sentence, mask_idx, masked_token, isnext))
      else
        put!(chn, (masked_sentence, mask_idx, masked_token))
      end
    end
  end
end

function _wrap_sentence(sentence1, sentence2...; start_token = "[CLS]", sep_token = "[SEP]")
  pushfirst!(sentence1, start_token)
  push!(sentence1, sep_token)
  map(s->push!(s, sep_token), sentence2)
  vcat(sentence1, sentence2...)
end

_wrap_idx(sentence1_idx, pre_len = 1) = sentence1_idx .+= pre_len
function _wrap_idx(sentence1_idx, sentence2_idx, len1)
  _wrap_idx(sentence1_idx)
  _wrap_idx(sentence2_idx, len1)
  vcat(sentence1_idx, sentence2_idx)
end

masksentence(sentence::String,
             tokenizer::Tf, wordpiece::WordPiece;
             mask_token = "[MASK]",
             mask_ratio = 0.15,
             real_token_ratio = 0.1,
             random_token_ratio = 0.1,
             whole_word_mask = false) where Tf = masksentence(
               tokenizer(sentence),
               wordpiece;
               mask_token = mask_token,
               mask_ratio = mask_ratio,
               real_token_ratio = real_token_ratio,
               random_token_ratio = random_token_ratio,
               whole_word_mask = whole_word_mask)

function masksentence(words,
                      wordpiece;
                      mask_token = "[MASK]",
                      mask_ratio = 0.15,
                      real_token_ratio = 0.1,
                      random_token_ratio = 0.1,
                      whole_word_mask = false)

  if whole_word_mask
    masked_word_idx = randsubseq(1:length(words), mask_ratio)

    is_masked = fill(false, length(words))
    map(masked_word_idx) do i
      is_masked[i] = true
    end

    tokens = Vector{String}()
    masked_idx = Vector{Int}()

    sizehint!(tokens, length(words))
    sizehint!(masked_idx, length(masked_word_idx))

    len = 0
    for (word, im) ∈ zip(words, is_masked)
      wps = wordpiece(word)
      wpl = length(wps)
      append!(tokens, wps)
      im && append!(masked_idx, map(i->i+len, 1:wpl))
      len += wpl
    end
  else
    tokens = wordpiece(words)
    masked_idx = randsubseq(1:length(tokens), mask_ratio)
  end

  masked_tokens = copy(tokens)

  for idx ∈ masked_idx
    r = rand()
    if r <= random_token_ratio
      masked_tokens[idx] = rand(wordpiece.vocab)
    elseif r > real_token_ratio + random_token_ratio
      masked_tokens[idx] = mask_token
    end
  end

  return masked_tokens, tokens, masked_idx
end


text:
"""
Strip from the right end of the `bytes` array for any byte that matches the ones
specified in the `remove` argument.  See Python's bytes.rstrip function.
"""
function brstrip(bytes::AbstractVector{UInt8}, remove::AbstractVector{UInt8})
    for i in length(bytes):-1:1
        x = bytes[i]
        found = false
        for m in remove
            if x == m
                found = true
                break
            end
        end
        if !found
            return bytes[1:i]
        end
    end
    return Vector{UInt8}()
end

# """
# Faster version of rstrip (slightly modified version from julia master branch)
# """
# function rstrip2(s::String)
#     i = endof(s)
#     while 1 ≤ i
#         c = s[i]
#         j = prevind(s, i)
#         c == ' ' || return s[1:i]
#         i = j
#     end
#     EMPTY_STRING
# end

"""
Find needle in the haystack with both `Vector{UInt8}` type arguments.
"""
function contains(haystack::AbstractVector{UInt8}, needle::AbstractVector{UInt8})
    hlen = length(haystack)
    nlen = length(needle)
    if hlen >= nlen
        for i in 1:hlen-nlen+1
            if haystack[i:i+nlen-1] == needle
                return true
            end
        end
    end
    return false
end

# Fast implementation to `reinterpret` int/floats
# See https://discourse.julialang.org/t/newbie-question-convert-two-8-byte-values-into-a-single-16-byte-value/7662/5

# Version a.  Original implementation... slow.
# """
# Byte swap is needed only if file the array represent a different endianness
# than the system.  This function does not make any assumption and the caller
# is expected to pass `true` to the `swap` argument when needed.
# """
# function convertfloat64a(bytes::Vector{UInt8}, swap::Bool)
#     # global count_a
#     # count_a::Int64 += 1
#     values = [bytes[i:i+8-1] for i in 1:8:length(bytes)]
#     values = map(x -> reinterpret(Float64, x)[1], values)
#     swap ? bswap.(values) : values
# end

# Version b.  Should be a lot faster.
# julia> @btime convertfloat64b(r, :LittleEndian);
#   370.677 μs (98 allocations: 395.41 KiB)
# """
# It turns out that `reinterpret` consider a single UInt64 as BigEndian
# Hence it's necessary to swap bytes if the array is in LittleEndian convention.
# This function does not make any assumption and the caller
# is expected to pass `true` to the `swap` argument when needed.
# """
# function convertfloat64b(bytes::Vector{UInt8}, endianess::Symbol)
#     # global count_b
#     # count_b::Int64 += 1
#     v = endianess == :LittleEndian ? reverse(bytes) : bytes
#     c = convertint64.(v[1:8:end],v[2:8:end],v[3:8:end],v[4:8:end],
#             v[5:8:end], v[6:8:end], v[7:8:end], v[8:8:end])
#     r = reinterpret.(Float64, c)
#     endianess == :LittleEndian ? reverse(r) : r
# end

# Version c
# julia> @btime convertfloat64c(r, :LittleEndian);
#   75.835 μs (2 allocations: 78.20 KiB)
# function convertfloat64c(bytes::Vector{UInt8}, endianess::Symbol)
#     L = length(bytes)
#     n = div(L, 8)               # numbers to convert
#     r = zeros(Float64, n)       # results
#     j = 1                       # result index
#     @inbounds for i in 1:8:L
#         if endianess == :LittleEndian
#             r[j] = reinterpret(Float64, convertint64(
#                         bytes[i+7], bytes[i+6], bytes[i+5], bytes[i+4],
#                         bytes[i+3], bytes[i+2], bytes[i+1], bytes[i]))
#         else
#             r[j] = reinterpret(Float64, convertint64(
#                         bytes[i],   bytes[i+1], bytes[i+2], bytes[i+3],
#                         bytes[i+4], bytes[i+5], bytes[i+6], bytes[i+7]))
#         end
#         j += 1
#     end
#     r
# end

# Version d
# julia> @btime convertfloat64d(r, :LittleEndian);
#   184.463 μs (4 allocations: 156.47 KiB)
# function convertfloat64d(bytes::Vector{UInt8}, endianess::Symbol)
#     if endianess == :LittleEndian
#         v = reverse(bytes)
#     else
#         v = bytes
#     end
#     L = length(bytes)
#     n = div(L, 8)               # numbers to convert
#     r = zeros(Float64, n)       # results
#     j = n                       # result index
#     for i in 1:8:L
#         r[j] = reinterpret(Float64, convertint64(
#             v[i],   v[i+1], v[i+2], v[i+3],
#             v[i+4], v[i+5], v[i+6], v[i+7]))
#         j -= 1
#     end
#     r
# end

# Version e - same as version c except that it does bit-shifting
# inline here in a new way.
# julia> @btime convertfloat64e(r, :LittleEndian);
#   174.685 μs (2 allocations: 78.20 KiB)
# function convertfloat64e(bytes::Vector{UInt8}, endianess::Symbol)
#     L = length(bytes)
#     n = div(L, 8)               # numbers to convert
#     r = zeros(Float64, n)       # results
#     j = 1                       # result index
#     for i in 1:8:L
#         v = UInt64(0)
#         if endianess == :LittleEndian
#             for k in 7:-1:1
#                 v = (v | bytes[i + k]) << 8
#             end
#             v |= bytes[i]
#         else
#             for k in 0:6
#                 v = (v | bytes[i + k]) << 8
#             end
#             v |= bytes[i + 7]
#         end
#         r[j] = reinterpret(Float64, v)
#         j += 1
#     end
#     r
# end

# Version f.  Best one so far!
# julia> @btime convertfloat64f(r, :LittleEndian)
#   35.132 μs (2 allocations: 78.20 KiB)
#
# results will be updated directly in the provided array `r`
function convertfloat64f!(r::AbstractVector{Float64}, bytes::Vector{UInt8}, endianess::Symbol)
    L = length(bytes)
    n = div(L, 8)               # numbers to convert
    j = 1                       # result index
    @inbounds for i in 1:8:L
        if endianess == :LittleEndian
            r[j] = reinterpret(Float64,
                        UInt64(bytes[i+7]) << 56 |
                        UInt64(bytes[i+6]) << 48 |
                        UInt64(bytes[i+5]) << 40 |
                        UInt64(bytes[i+4]) << 32 |
                        UInt64(bytes[i+3]) << 24 |
                        UInt64(bytes[i+2]) << 16 |
                        UInt64(bytes[i+1]) << 8  |
                        UInt64(bytes[i]))
        else
            r[j] = reinterpret(Float64,
                        UInt64(bytes[i]) << 56 |
                        UInt64(bytes[i+1]) << 48 |
                        UInt64(bytes[i+2]) << 40 |
                        UInt64(bytes[i+3]) << 32 |
                        UInt64(bytes[i+4]) << 24 |
                        UInt64(bytes[i+5]) << 16 |
                        UInt64(bytes[i+6]) << 8  |
                        UInt64(bytes[i+7]))
        end
        j += 1
    end
    r
end

# Conversion routines for 1,2,4,8-byte words into a single 64-bit integer
@inline function convertint64B(a::UInt8,b::UInt8,c::UInt8,d::UInt8,e::UInt8,f::UInt8,g::UInt8,h::UInt8)
    (Int64(a) << 56) | (Int64(b) << 48) | (Int64(c) << 40) | (Int64(d) << 32) |
    (Int64(e) << 24) | (Int64(f) << 16) | (Int64(g) << 8)  |  Int64(h)
end
@inline function convertint64L(a::UInt8,b::UInt8,c::UInt8,d::UInt8,e::UInt8,f::UInt8,g::UInt8,h::UInt8)
    (Int64(h) << 56) | (Int64(g) << 48) | (Int64(f) << 40) | (Int64(e) << 32) |
    (Int64(d) << 24) | (Int64(c) << 16) | (Int64(b) << 8)  |  Int64(a)
end
@inline function convertint64B(a::UInt8,b::UInt8,c::UInt8,d::UInt8)
    (Int64(a) << 24) | (Int64(b) << 16) | (Int64(c) << 8)  |  Int64(d)
end
@inline function convertint64L(a::UInt8,b::UInt8,c::UInt8,d::UInt8)
    (Int64(d) << 24) | (Int64(c) << 16) | (Int64(b) << 8)  |  Int64(a)
end
@inline function convertint64B(a::UInt8,b::UInt8)
    (Int64(a) << 8)  |  Int64(b)
end
@inline function convertint64L(a::UInt8,b::UInt8)
    (Int64(b) << 8)  | Int64(a)
end

# this version is slightly slower
# function convertint64b(a::UInt8,b::UInt8,c::UInt8,d::UInt8,e::UInt8,f::UInt8,g::UInt8,h::UInt8)
#     v = UInt64(a) << 8
#     v |= b ; v <<= 8
#     v |= c ; v <<= 8
#     v |= d ; v <<= 8
#     v |= e ; v <<= 8
#     v |= f ; v <<= 8
#     v |= g ; v <<= 8
#     v |= h
#     v
# end

# TODO cannot use AbstractString for some reasons
# """
# Concatenate an array of strings to a single string
# """
# concatenate(strArray::Vector{T} where T <: AbstractString, separator=",") =
#     foldl((x, y) -> *(x, y, separator), "", strArray)[1:end-length(separator)]

# """
# Convert a dictionary to an array of k=>v strings
# """
# stringarray(dict::Dict) =
#     ["$x => $y" for (x, y) in dict]


text:
# Autogenerated wrapper script for Rmath_jll for x86_64-unknown-freebsd
export libRmath

JLLWrappers.@generate_wrapper_header("Rmath")
JLLWrappers.@declare_library_product(libRmath, "libRmath-julia.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libRmath,
        "lib/libRmath-julia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Logging
io = open("test_log.txt", "w+")
logger = SimpleLogger(io)
with_logger(logger) do
    testdirectory = dirname(@__FILE__)
    cd(testdirectory)
    directories = readdir()
    for directory in directories
        if isdir(directory)
            include(joinpath(directory, "runtest.jl"))
        end
        cd(testdirectory)
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using Test
using Random
using LinearAlgebra
using SparseArrays
using MatrixDepot
using JLD2, FileIO
using TraceEstimation

@testset "SLQ" begin
    Random.seed!(123543)
    @testset "Dense SPD Matrices" begin
        @testset "a[i,j] = exp(-2 * abS(i - j)) (small size)" begin
            println("Executing Test 01: Dense SPD - Small Size - Inverse")
            A = rand(610, 610)
            for i in 1:610
                for j in 1:610
                    A[i, j] = exp(-2 * abs(i - j))
                end
            end
            obv = slq(A, fn = inv, ctol = 0.01, eps = 0.05, mtol = 0.01)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=0.01)
        end
        @testset "a[i,j] = exp(-2 * abS(i - j)) (large size)" begin
            println("Executing Test 02: Dense SPD - Large Size - Inverse")
            A = rand(4610, 4610)
            for i in 1:4610
                for j in 1:4610
                    A[i, j] = exp(-2 * abs(i - j))
                end
            end
            obv = slq(A, fn = inv, ctol = 0.01, eps = 0.05, mtol = 0.01)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=0.01)
        end
        @testset "Random SPD Matrix (large size)" begin
            println("Executing Test 03: Random Dense SPD - Large Size - Log")
            A = rand(4000, 4000)
            A = A + A'
            while !isposdef(A)
                A = A + 30I
            end
            A = A + 30I
            obv = slq(A, fn = log, ctol = 0.01)
            acv = tr(log(A))
            @test isapprox(obv, acv, rtol=0.01)
        end
        @testset "Hilbert Matrix" begin
            println("Executing Test 04: Dense Hilbert SPD - Inverse")
            A = matrixdepot("hilb", 3000)
            A = A + I
            obv = slq(A, fn = inv, ctol = 0.01)
            acv = tr(inv(A))
            @test isapprox(obv, acv, rtol=0.01)
        end
    end
    @testset "Sparse SPD Matrices" begin
        @testset "Sparse Matrix" begin
            println("Executing Test 05: Random Sparse Matrix - SQRT")
            A = sprand(6000, 6000, 0.001)
            A = A + A'
            while !isposdef(A)
                A = A + 60I
            end
            obv = slq(A, fn = sqrt, ctol = 0.01)
            M = Matrix(A)
            acv = tr(sqrt(M))
            @test isapprox(obv, acv, rtol = 0.01)
        end
        @testset "Wathen Sparse Matrix" begin
            println("Executing Test 06: Finite Element Matrix - Wathen")
            A = matrixdepot("wathen", 40) #7601x7601
            obv = slq(A, fn = inv, eps = 0.05)
            M = Matrix(A)
            acv = tr(inv(M))
            @test isapprox(obv, acv, rtol = 0.1)
        end
        @testset "TopOpt Shift +1 Matrix" begin
            println("Executing Test 07: TopOpt Shift +1")
            println("Loading from JLD2 file")
            #=
            s = (40, 10) # increase to increase the matrix size
            xmin = 0.9 # decrease to increase the condition number
            problem = HalfMBB(Val{:Linear}, s, (1.0, 1.0), 1.0, 0.3, 1.0)
            solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
            n = length(solver.vars)
            solver.vars[rand(1:n, n÷2)] .= 0
            solver()
            K = solver.globalinfo.K
            K = K+1*I
            =#
            @load "topoptfile.jld2" K
            obv = slq(K, fn = inv, eps = 0.005)
            M = Matrix(K)
            acv = tr(inv(M))
            @test isapprox(obv, acv, rtol=0.01)
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using Test

@testset "test integrate" begin
    element = Element(Seg2, [1, 2])
    ips1 = get_integration_points(element)
    ips2 = get_integration_points(element, 1)
    @test length(ips1) == 2
    @test length(ips2) == 3
    element = Element(Poi1, [1])
    ips3 = get_integration_points(element)
    @test length(ips3) == 1
end


text:
import Plots: heatmap, contour, contourf

"""
    heatmap(fld::MeshArray; args...)

Calls `heatmap` for one elementary array per subplot
"""
function heatmap(x::MeshArray; args...)
    n=x.grid.nFaces
    isa(x,MeshArrays.gcmarray) ? n=length(x.fIndex) : nothing

    m=Int(sqrt(n))
    j=collect(1:n)
    j=reshape(collect(j),(m,m))
    j=reverse(permutedims(j,(2,1)); dims=2)[:]

    p=()
    for i=1:n; p=(p...,heatmap(x.f[j[i]]; args...)); end
    plot(p...)
end

"""
    contour(fld::MeshArray; args...)

Calls `contour` for one elementary array per subplot
"""
function contour(x::MeshArray; args...)
    n=x.grid.nFaces
    isa(x,MeshArrays.gcmarray) ? n=length(x.fIndex) : nothing
    p=()
    for i=1:n; p=(p...,contour(x.f[i]; args...)); end
    plot(p...)
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

MyGame.Example.eval(quote

@enum Race::Int8 begin
    RaceNone = -1
    RaceHuman = 0
    RaceDwarf = 1
    RaceElf = 2
end

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia