text:
using NLPModels
using LinearOperators
using OptimizationProblems
import CSV
import DataFrames

# Convex problems
# Problems = [arglina, arglinb, arglinc, arwhead, bdqrtic, brownden, cliff, clplatea, clplateb, clplatec, dixon3dq, dqdrtic, dqrtic, engval1, nasty, nondquar,
# palmer1c, palmer1d, palmer2c, palmer3c, palmer4c, palmer5c, palmer5d, palmer6c, palmer7c, palmer8c, power, quartc, tridia, vardim]

# Convex problems var ≥ 10 (17 pbs)
# Problems = [arglina, arglinb, arglinc, arwhead, bdqrtic, clplatea, clplateb, clplatec,
# dixon3dq, dqdrtic, dqrtic, engval1, nondquar, power, quartc, tridia, vardim]

# Nonconvex problems
# Problems = [AMPGO02, AMPGO03, AMPGO04, AMPGO05, AMPGO06, AMPGO08, AMPGO09,
# AMPGO10, AMPGO11, AMPGO12, AMPGO13, AMPGO14, AMPGO15, AMPGO18,
# AMPGO20, AMPGO21, AMPGO22,
# beale, brownbs, broydn7d, brybnd,
# chainwoo, chnrosnb_mod, cosine, cragglvy, curly, dixmaane, dixmaanf, dixmaang,
# dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl, dixmaanm, dixmaann, dixmaano, dixmaanp,
# Dus2_1, Dus2_3, Dus2_9, Duscube,
# edensch, eg2, errinros_mod, extrosnb, fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose, genrose_nash,
# indef_mod, liarwhd, meyer3, morebv, ncb20, ncb20b, noncvxu2, noncvxun, nondia, penalty2, powellsg, schmvett,
# scosine, Shpak1, Shpak2, Shpak3, Shpak4, Shpak5, Shpak6, sinquad, sparsine, sparsqur, srosenbr, tointgss, tquartic, woods];

# Nonconvex problems of dimension > 1
# Problems = [beale, brownbs, broydn7d, brybnd,
# chainwoo, chnrosnb_mod, cosine, cragglvy, curly, dixmaane, dixmaanf, dixmaang,
# dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl, dixmaanm, dixmaann, dixmaano, dixmaanp,
# edensch, eg2, errinros_mod, extrosnb, fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose, genrose_nash,
# indef_mod, liarwhd, meyer3, morebv, ncb20, ncb20b, noncvxu2, noncvxun, nondia, penalty2, powellsg, schmvett,
# scosine, sinquad, sparsine, sparsqur, srosenbr, tointgss, tquartic, woods]

# Nonconvex problems var ≥ 10 (50 pbs)
Problems = [broydn7d, brybnd, chainwoo, chnrosnb_mod, cosine, cragglvy, curly,
dixmaane, dixmaanf, dixmaang, dixmaanh, dixmaani, dixmaanj, dixmaank, dixmaanl,
dixmaanm, dixmaann, dixmaano, dixmaanp, edensch, eg2, errinros_mod, extrosnb,
fletcbv2, fletcbv3_mod, fletchcr, fminsrf2, freuroth, genhumps, genrose,
genrose_nash, indef_mod, liarwhd, morebv, ncb20, ncb20b, noncvxu2, noncvxun,
nondia, penalty2, powellsg, schmvett, scosine, sinquad, sparsine, sparsqur,
srosenbr, tointgss, tquartic, woods]

# All pbs var ≥ 10 and not in CUTEst
# Problems = [clplatea, clplateb, clplatec, chnrosnb_mod, errinros_mod,
# fletcbv3_mod, genrose_nash, scosine]

algo = CRlin
# csvname = string("noncvxlin",string(algo)[1:2],".csv")

D = ["model     algo nvar   f(x)    f(x0)   ‖g(x)‖  ‖g(x0)‖   #f  #g #Hv  #it optimal"]

MiniLogging.@info(loggerlaunch, D[1])
n = 100
for problem in Problems
    model = MathProgNLPModel(problem(n), name=string(problem))
    L = Linesearch(model, algo)
    S = @sprintf("%5s %5s %4d %5s %5s %5s %5s %4d %4d %4d %4d %5s", L[1], L[2], L[3], L[4], L[5], L[6], L[7], L[8], L[9], L[10], L[11], L[12])
    MiniLogging.@info(loggerlaunch, S)
    D = vcat(D, S)
    reset!(model)
end
# DF = DataFrames.DataFrame(D)
# CSV.write(csvname, DF)
# writedlm("ResultsCGlin.txt", D)
writedlm("nonconvpbsCRlin.txt", D)


text:
module TestWrappedFunctions

# using Revise
using Test
using MLJBase
using ..Models
using CategoricalArrays
import Random.seed!
seed!(1234)

struct YourTransformer <: Static
    ftr::Symbol
end

MLJBase.transform(transf::YourTransformer, verbosity, X) =
    selectcols(X, transf.ftr)

@testset "nodal machine constructor for static transformers" begin
    X = (x1=rand(3), x2=[1, 2, 3]);
    mach = machine(YourTransformer(:x2))
    fit!(mach, verbosity=0)
    @test transform(mach, X) == [1, 2, 3]
end

x1 = rand(30)
x2 = rand(30)
x3 = rand(30)
y = exp.(x1 - x2 -2x3 + 0.1*rand(30))
X = (x1=x1, x2=x2, x3=x3)

f(X) = (a=selectcols(X, :x1), b=selectcols(X, :x2))

knn = KNNRegressor()

# 1. function in a pipeline:
comp1 = @pipeline f Standardizer knn target=UnivariateBoxCoxTransformer
e = evaluate(comp1, X, y, measure=mae, resampling=Holdout(), verbosity=0)

# 2. function with parameters in a pipeline:
mutable struct GreatTransformer <: Static
    ftr::Symbol
end
MLJBase.transform(transf::GreatTransformer, verbosity, X) =
    (a=selectcols(X, transf.ftr), b=selectcols(X, :x2))
comp2 = @pipeline(GreatTransformer(:x3), Standardizer, knn,
                      target=UnivariateBoxCoxTransformer)

comp2.great_transformer.ftr = :x1 # change the parameter
e2 =  evaluate(comp2, X, y, measure=mae, resampling=Holdout(), verbosity=0)
@test e2.measurement[1] ≈ e.measurement[1]

# 3. function in an exported learning network:
Xs = source(X)
ys = source(y)
f(X::AbstractNode) = node(f, X)
# or, without arrow syntax:
# W = Xs |> f |> Standardizer()
# z = ys |> UnivariateBoxCoxTransformer()
# zhat = (W, z) |> knn
# yhat = zhat |> inverse_transform(z)
X2 = f(Xs)
W = transform(machine(Standardizer(), X2), X2)
box_mach = machine(UnivariateBoxCoxTransformer(), ys)
z = transform(box_mach, ys)
knn_mach = machine(knn, W, z)
zhat = predict(knn_mach, W)
yhat = inverse_transform(box_mach, zhat)

@from_network machine(Deterministic(), Xs, ys; predict=yhat) begin
    mutable struct Comp3
        rgs=knn
    end
end

comp3 = Comp3()
e3 =  evaluate(comp3, X, y, measure=mae, resampling=Holdout(), verbosity=0)
@test e2.measurement[1] ≈ e.measurement[1]

# 4. function with parameters in exported learning network:
inserter = GreatTransformer(:x3)
# W = Xs |> inserter |> Standardizer()
# z = ys |> UnivariateBoxCoxTransformer()
# zhat = (W, z) |> knn
# yhat = zhat |> inverse_transform(z)
inserter_mach = machine(inserter)
X2 = transform(inserter_mach, Xs)
W = transform(machine(Standardizer(), X2), X2)
box_mach = machine(UnivariateBoxCoxTransformer(), ys)
z = transform(box_mach, ys)
knn_mach = machine(knn, W, z)
zhat = predict(knn_mach, W)
yhat = inverse_transform(box_mach, zhat)
@from_network machine(Deterministic(), Xs, ys; predict=yhat) begin
    mutable struct Comp4
        transf=inserter
        rgs=knn
    end
end
comp4 = Comp4()
comp4.transf.ftr = :x1 # change the parameter
e4 =  evaluate(comp4, X, y, measure=mae, resampling=Holdout(), verbosity=0)
@test e4.measurement[1] ≈ e.measurement[1]

end

true

ext: jl


text:
@lexer MatlabLexer let
    keywords = [
        "break", "case", "catch", "classdef", "continue", "else", "elseif",
        "end", "enumerated", "events", "for", "function", "global", "if",
        "methods", "otherwise", "parfor", "persistent", "properties",
        "return", "spmd", "switch", "try", "while",
    ]
    Dict(
        :name => "Matlab",
        :description => "A lexer for Matlab source files.",
        :comments => "Based on the lexer from Pygments.",
        :aliases => ["matlab"],
        :tokens => Dict(
            :root => [
                (r"^!.*"m, STRING_OTHER),
                (r"%\{\s*$"ms, COMMENT_MULTILINE, :blockcomment),
                (r"%.*$"m, COMMENT_SINGLE),
                (r"^\s*function", KEYWORD, :deffunc),

                (words(keywords, suffix = "\\b"), KEYWORD),

                (r"\.\.\..*$"m, COMMENT),

                (r"(-|==|~=|<|>|<=|>=|&&|&|~|\|\|?|\^)", OPERATOR),
                (r"(\.\*|\*|\+|\.\^|\.\\|\.\/|\/|\\)", OPERATOR),

                (r"(\[|\]|\(|\)|\{|\}|:|@|\.|,)", PUNCTUATION),
                (r"(=|:|;)", PUNCTUATION),

                (r"(?<=[\w)\].])'+", OPERATOR),

                (r"(\d+\.\d*|\d*\.\d+)([eEf][+-]?[0-9]+)?", NUMBER_FLOAT),
                (r"\d+[eEf][+-]?[0-9]+", NUMBER_FLOAT),
                (r"\d+", NUMBER_INTEGER),

                (r"(?<![\w)\].])'", STRING, :string),
                (r"([a-zA-Z_]\w*)(\()", (NAME_FUNCTION, PUNCTUATION)),
                (r"[a-zA-Z_]\w*", NAME),
                (r".", TEXT),
            ],
            :blockcomment => Any[
                (r"^\s*%\}"ms, COMMENT_MULTILINE, :__pop__),
                (r"."sm, COMMENT_MULTILINE),
            ],
            :deffunc => [
                (
                    r"(\s*)(?:(.+)(\s*)(=)(\s*))?(.+)(\()(.*)(\))(\s*)",
                    (
                        WHITESPACE, TEXT, WHITESPACE, PUNCTUATION, WHITESPACE,
                        NAME_FUNCTION, PUNCTUATION, TEXT, PUNCTUATION, WHITESPACE,
                    ),
                    :__pop__,
                ),
                (r"(\s*)([a-zA-Z_]\w*)", (TEXT, NAME_FUNCTION), :__pop__),
            ],
            :string => [
                (r"[^']*'", STRING, :__pop__),
            ],
        ),
    )
end


text:
using Printf
using Statistics
using LinearAlgebra
using MDBM
using PyPlot
pygui(true);
using Interpolations
using Plots

function evplot(a,k)
U=a[:,k]
eigRe1=[real(U[n,1]) for n in 1:size(U,1)]
eigIm1=[imag(U[n,1]) for n in 1:size(U,1)]
plot1=Plots.plot(eigRe1,eigIm1,seriestype=:scatter)
Re_circle=[cos(n) for n in 0:0.05:2*pi]
Im_circle=[sin(n) for n in 0:0.05:2*pi]
return(Plots.plot!(plot1,Re_circle,Im_circle))
end

function check(a)
    v=false
    s=maximum(size(a))
    dist=zeros(s)
    aRe=real(a)
    aIm=broadcast(abs,imag(a))
    for j=1:s
        if abs(aIm[j])<mar
        else
            for k=1:s
                if k==j
                dist[k]=1
                else
                dist[k]=norm((aRe[j]+im*aIm[j])-(aRe[k]+im*aIm[k]))
                end
            end
            if minimum(dist)<mar
            else v=true
            end
        end
    end
    return(v)
end

function it(A)
    matrdim=size(A,2)-1
    scaleitp=real(A[:,1])
    knots=(scaleitp,)
    ARe=real(A)
    AIm=imag(A)

    imRe=[interpolate(knots, ARe[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imRetemp=[interpolate(knots, ARe[:,j], Gridded(Linear()))]
        imRe=vcat(imRe,imRetemp)
    end
    imIm=[interpolate(knots, AIm[:,2],Gridded(Linear()))]
    for j = 3:matrdim+1
        imImtemp=[interpolate(knots, AIm[:,j], Gridded(Linear()))]
        imIm=vcat(imIm,imImtemp)
    end
    return(hcat(imRe,imIm))
end

#substitution in interpolating function
function sub(it,t)
    subdim=size(it,1)
    out=zeros(ComplexF64,subdim)
    for j = 1:subdim
        out[j]=it[j,1](t)+it[j,2](t)*im
    end
    return(out)
end

function findnormmax(a)
    s=maximum(size(a))
    norma=zeros(s)
    for j = 1:s
        norma[j]=norm(a[j])
    end
    return(a[findmax(norma)[2]])
end

function control(a,r)
    s=maximum(size(a))
    norma=zeros(s)
    for j = 1:s
        norma[j]=norm(a[j])
    end
        max1=findmax(norma)[2]
        norma[max1]=0
        max2=findmax(norma)[2]
        norma[max2]=0
    var=true
    for j=1:s
        if norma[j]<r
        else var=false
        end
    end
    if (abs(imag(a[max1]))<mar && abs(imag(a[max2]))<mar) || abs((real(a[max1])-real(a[max2])))<mar
    else var=false
    end
    return(var)
end

dim=2 #DoF of the system (after Cauchy transcription)
mult=16 #Ns=dim*mult
n=40 #timestep number
mar=0.0001
#choosing numerical simulation method
#RK4: 4th order Runge-Kutta
#EI: Implicit Euler
method="RK4"
i=50 #number of iteration
r=n-1

#parameters of the delayed Mathieu equation
tau=2*pi
omega=1
kappa=0.1
delta=0.5
epsilon=1
b=1.2
T=2*pi/omega

#timestep
dt=tau/(n-1)

#matrix coefficients
function A(t,vv)
    Amatr=zeros(dim,dim)
    #vv[1]=kappa
    #vv[2]=delta
    #vv[3]=epsilon
    #vv[4]=omega
    Amatr=[-vv[1] -(vv[2]+vv[3]*cos(vv[4]*t)); 1 0]
    return(Amatr)
end

function B(bv)
    Bmatr=zeros(dim*dim)
    Bmatr=[0 bv; 0 0]
    return(Bmatr)
end

function simulation(kappa1,delta1,epsilon1,omega1,b1)
        v=[kappa1, delta1, epsilon1, omega1]
        solreturn=zeros(ComplexF64,1,1+mult*dim) #matrix for time integration solution (test)
        eval=zeros(ComplexF64,mult,1)
        S1=zeros(ComplexF64,n,mult*dim)
        V1=zeros(ComplexF64,n,mult*dim)
        S=zeros(ComplexF64,n*dim,mult)
        V=zeros(ComplexF64,n*dim,mult)
        Vj=zeros(ComplexF64,n*dim,mult)
        Vjnorm=zeros(ComplexF64,n*dim,mult)
        H=zeros(ComplexF64,mult,mult)
        tvect0=collect(-tau:dt:tau+dt/2) #discretized time vector
        sol=zeros(ComplexF64,(n-1),mult*dim) #empty solution vector
        sol00=rand(n,mult*dim)*(-2)+2*ones(n,mult*dim) #S0 inital random array
        g=1
        #iteration loop
        for g = 1:i
        sol0=hcat(tvect0,vcat(sol00,sol)) #constructing solution array
            #time integration
            for m = 2:dim:(dim*mult-(dim-2)) #choosing "one physical system"
            sol0m=hcat(sol0[:,1],sol0[:,m:m+dim-1])
                #time integration steps
            if method == "EE"
                for j = 0:n-2

                    Y=sol0m[n+j,2:2+(dim-1)]+dt*(A(sol0m[n+j,1],v)*sol0m[n+j,2:2+(dim-1)]+B(b1)*sol0m[1+j,2:2+(dim-1)])
                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y) #filling solution array of the physical system

                end

            elseif method == "RK4_int"
                for j = 0:n-2
                    sol0mit=it(sol0m)

                    Y1=sub(sol0mit,j*dt)
                    Y2=Y1+(dt/2)*(A(j*dt,v)*Y1+B(b1)*sub(sol0mit,-tau+j*dt))
                    Y3=Y1+(dt/2)*(A((j+0.5)*dt,v)*Y2+B(b1)*sub(sol0mit,-tau+(j+0.5)*dt))
                    Y4=Y1+dt*(A((j+0.5)*dt,v)*Y3+B(b1)*sub(sol0mit,-tau+(j+0.5)*dt))

                    Y=Y1+(dt/6)*((A(j*dt,v)*Y1+B(b1)*sub(sol0mit,-tau+j*dt))+2*(A((j+0.5)*dt,v)*Y2+B(b1)*sub(sol0mit,-tau+(j+0.5)*dt))+2*(A((j+0.5)*dt,v)*Y3+B(b1)*sub(sol0mit,-tau+(j+0.5)*dt))+(A((j+1)*dt,v)*Y4+B(b1)*sub(sol0mit,-tau+(j+1)*dt)))

                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y) #filling solution array of the physical system

                end

            elseif method == "RK4"
                for j = 0:n-2

                    Ydtau=0.5*(sol0m[j+1,2:2+(dim-1)]+sol0m[j+2,2:2+(dim-1)])

                    Y1=sol0m[n+j,2:2+(dim-1)]
                    Y2=Y1+(dt/2)*(A(j*dt,v)*Y1+B(b1)*sol0m[j+1,2:2+(dim-1)])
                    Y3=Y1+(dt/2)*(A((j+0.5)*dt,v)*Y2+B(b1)*Ydtau)
                    Y4=Y1+dt*(A((j+0.5)*dt,v)*Y3+B(b1)*Ydtau)

                    Y=Y1+(dt/6)*((A(j*dt,v)*Y1+B(b1)*sol0m[j+1,2:2+(dim-1)])+2*(A((j+0.5)*dt,v)*Y2+B(b1)*Ydtau)+2*(A((j+0.5)*dt,v)*Y3+B(b1)*sol0m[j+1,2:2+(dim-1)])+(A((j+1)*dt,v)*Y4+B(b1)*sol0m[j+2,2:2+(dim-1)]))

                    sol0m[n+j+1,2:2+(dim-1)]=transpose(Y) #filling solution array of the physical system

                end
            end
                sol0[:,m:m+dim-1]=sol0m[:,2:2+(dim-1)] #filling iterational solution array
            end
            #solreturn=vcat(solreturn,sol0) #filling global solution array (test)
            S1=sol0[1:n,2:mult*dim+1]
            V1=sol0[n:2*n-1,2:mult*dim+1]
                for p=1:n
                    for s=1:mult
                        for q=1:dim
                            S[1+(p-1)*dim+(q-1),s]=S1[p,1+(s-1)*dim+(q-1)]
                            V[1+(p-1)*dim+(q-1),s]=V1[p,1+(s-1)*dim+(q-1)]
                        end
                    end
                end
            H=pinv(S)*V
            Am=eigvecs(H) #calculation eigenvectors of Hj, creating matrix Am (Ns x Ns)
            for o = 1:n*dim
                Vj[o,:]=Am*V[o,:]
            end
            for h = 1:mult
                Vjnorm[:,h]=normalize(Vj[:,h])
            end
            sol00=zeros(ComplexF64,n,mult*dim)
            for p=1:n
                for s=1:mult
                    for q=1:dim
                             sol00[p,1+(s-1)*dim+(q-1)]=Vjnorm[1+dim*(p-1)+(q-1),s]
                    end
                end
            end
            #sol00=sol1
            evm=eigvals(H) #obtaining dominant eigenvalues in the iteration step
            eval=hcat(eval,evm)      #filling up the eigenvalue array
            #if (control(evm,mar) && g>3) break  end
            if check(evm) print(g-1); break end
        end
        #evm=eigvals(H)
         return(norm(findnormmax(eval[:,end-1]))-1)
        #return(norm(findnormmax(evm))-1)
        #return(eval[:,2:end]) #returning the global eigenvalue array
        #return(solreturn)
end

Sim=simulation(kappa,delta,epsilon,omega,b)

evplot(Sim,1)
evplot(Sim,2)
evplot(Sim,3)
evplot(Sim,4)
evplot(Sim,5)
evplot(Sim,6)
evplot(Sim,7)
evplot(Sim,8)
evplot(Sim,10)
evplot(Sim,20)
evplot(Sim,15)
evplot(Sim,20)
evplot(Sim,30)
evplot(Sim,100)
evplot(Sim,200)
evplot(Sim,400)

check(Sim[:,4])

function foo(x,y)
    v=[kappa, x, epsilon, omega]
    return(simulation(v[1],v[2],v[3],v[4],y))
end

ax1=Axis([-1.0,0,1,2,3,4,5],"delta") # initial grid in x direction
ax2=Axis(-1.5:0.5:1.5,"b") # initial grid in y direction

mymdbm=MDBM_Problem(foo,[ax1,ax2])
iteration=4 #number of refinements (resolution doubling)
@time solve!(mymdbm,iteration)

x_eval,y_eval=getevaluatedpoints(mymdbm)
x_sol,y_sol=getinterpolatedsolution(mymdbm)
fig = figure(1);clf()
PyPlot.scatter(x_eval,y_eval,s=5)
PyPlot.scatter(x_sol,y_sol,s=5)

[1,2,5,-4,5][end-1]


text:
module FunctionOperators
using Parameters, MacroTools, Printf, LinearAlgebra
import Base: show, *, +, -, adjoint, ==, eltype, reshape

export FunctionOperator, FunOp, reshape, setPlan, @♻, @recycle, FunctionOperators_global_settings

# Structure to hold global settings
mutable struct Settings
    verbose::Bool
    macro_verbose::Bool
    auto_reshape::Bool
end
"""
Object that holds global settings for `FunctionOperators` library

Fields:
 - `verbose::Bool` If set to true, then allocation information and calculated plan function will be displayed upon creation (i.e., when a composite operator is first used). Default: `false`
 - `macro_verbose::Bool` If set to true, then recycling macros (@♻ and @recycle) will print the transformed code. Default: `false`
 - `auto_reshape::Bool` If set to true, then input and output is reshaped according to the inDims and outDims values of the FunctionOperator before and after any multiplication. Default: `false`
"""
const FunctionOperators_global_settings = Settings(false, false, false)

include("StructDefs.jl")
include("Helpers.jl")
include("Show.jl")
include("BuildCompTree.jl")
include("getPlan.jl")
include("mul.jl")
include("Auxiliary.jl")
include("recycle.jl")

function __init__()
    empty!(bufferPool)
    for _ in 1:Threads.nthreads()
        push!(bufferPool, WeakKeyDict{Buffer, Int}())
    end
end

end # module


text:
checkDir()
cwd = pwd()
#fix the seed
Random.seed!(1)

# only use lineageMarkers for clustering
(lineageMarkers,)= getMarkers(panel)
cc = map(Symbol, lineageMarkers)
dfSom = daf.fcstable[:,cc]

som2 = initGigaSOM(dfSom, 10, 10)

@testset "Dimensions - batch" begin
    @test size(som2.codes) == (100,10)
    @test som2.xdim == 10
    @test som2.ydim == 10
    @test som2.numCodes == 100
end

# using batch som with epochs
som2 = trainGigaSOM(som2, dfSom, epochs = 1)

winners = mapToGigaSOM(som2, dfSom)

embed = embedGigaSOM(som2, dfSom, k=10, smooth=0.0, adjust=0.5)

#test batch
@testset "Batch" begin
    codes = som2.codes
    @test size(codes) == (100,10)

    dfCodes = DataFrame(codes)
    names!(dfCodes, Symbol.(som2.colNames))
    dfEmbed = DataFrame(embed)
    CSV.write(genDataPath*"/batchDfCodes.csv", dfCodes)
    CSV.write(genDataPath*"/batchWinners.csv", winners)
    CSV.write(genDataPath*"/batchEmbedded.csv", dfEmbed)

    #load the ref data
    refBatchDfCodes = CSV.File(refDataPath*"/refBatchDfCodes.csv") |> DataFrame
    refBatchWinners = CSV.File(refDataPath*"/refBatchWinners.csv") |> DataFrame
    refBatchEmbedded = CSV.File(refDataPath*"/refBatchEmbedded.csv") |> DataFrame

    #load the generated data
    batchDfCodes = CSV.File(genDataPath*"/batchDfCodes.csv") |> DataFrame
    batchDfCodesTest = first(batchDfCodes, 10)
    batchWinners = CSV.File(genDataPath*"/batchWinners.csv") |> DataFrame
    batchWinnersTest = first(batchWinners, 10)
    batchEmbedded = CSV.File(genDataPath*"/batchEmbedded.csv") |> DataFrame
    batchEmbeddedTest = first(batchEmbedded, 10)

    # test the generated data against the reference data
    @test refBatchWinners == batchWinnersTest
    @test refBatchDfCodes == batchDfCodesTest
    @test Array{Float64,2}(refBatchEmbedded) ≈ Array{Float64,2}(batchEmbeddedTest) atol=1e-4

    for (i, j) in zip(batchDfCodesTest[:,1], refBatchDfCodes[:,1])
        @test isapprox(i, j; atol = 0.001)
    end
end


text:
using Conda

ENV["PYTHON"]=" "  # Configure PyCall.jl to use Conda.jl's Python

# We need install a specific version of llvmlite which supports the same minor version of
# llvm used with Julia. https://github.com/numba/llvmlite#compatibility
const LLVM_VERSION = VersionNumber(Base.libllvm_version)

const LLVMLITE_VERSION = if LLVM_VERSION > v"4.0.0"
    v"0.17.0"
elseif LLVM_VERSION > v"3.9.0"
    v"0.16.0"
elseif LLVM_VERSION > v"3.8.0"
    v"0.15.0"
elseif LLVM_VERSION > v"3.7.0"
    v"0.12.1"
elseif LLVM_VERSION > v"3.6.0"
    v"0.8.0"
else
    error("Can't find a version of llvmlite that supports $LLVM_VERSION")
end

Conda.add_channel("numba")  # channel containing llvmlite
Conda.add("llvmlite==$LLVMLITE_VERSION")  # https://github.com/numba/llvmlite

Conda.add_channel("invenia")
Conda.update()
Conda.add("pandapower")
Conda.add("numpy==1.12.1")

info("Verifying pandapower install.")
python_bin = joinpath(Conda.PYTHONDIR, "python")
run(`$python_bin -c "import pandapower, pandapower.networks, pandapower.topology, pandapower.converter, pandapower.estimation"`)
run(`$python_bin -c "from pandapower.networks import case9; case9()"`)


text:
mutable struct Chunk
	blks::SMatrix{64,64,Block,4096}
	ligs::SVector{UInt8,2048}
	tm::Int
end
function get_lig(c::Chunk,x::Int,y::Int)::UInt8
	id=(x<<5)|(y>>1)+1
	return (y&1==0) ? c.ligs[id]>>4 : c.ligs[id]&0xf
end
function set_lig(c::Chunk,x::Int,y::Int,v::UInt8)::UInt8
	id=(x<<5)|(y>>1)+1
	(y&1==0) ? c.ligs[id]=c.ligs[id]&0xf|(v<<4) : c.ligs[id]=c.ligs[id]&0xf0|v
end
function generate(c::Chunk,x::Int,y::Int;map_generator::MapGenerator)
	use(map_generator,c,x,y)
	if haskey(sav,Pair(x,y))
		t=@inbounds sav[Pair(x,y)]
		for i in t
			x0=i.first>>6
			y0=i.first&63
			c.blks[x0,y0]=i.second
		end
	end
end
const visible_width=480
const visible_height=480
const visible_size=32
# 不应使用特性改动
"将(lx,ly)左上角显示到画板(px,py)"
function draw(c::Chunk,lx::Int,ly::Int,px::Int,py::Int;envl::UInt8)
	rx=min(lx+15-px>>5,64)
	ry=min(ly+15-py>>5,64)
	for i in lx:rx
		for j in ly:ry
			tl=get_lig(c,i,j)
			if tl<envl
				tl=envl
			end
			x0=px+(i-lx)<<5
			y0=py+(j-ly)<<5
			if tl==0x0
				fill_rect(dcon,x0,y0,32,32;color=RGB(0,0,0))
				continue
			end
			b=c.blks[i,j]
			if hole(b) fill_rect(dcon,x0,y0,32,32;color=RGB(15/16,15/16,15/16))
			end
			b_show(b,dcon,x0,y0)
			if tl!=0x15
				fill_rect(dcon,x0,y0,32,32;color=RGBA(0,0,0,(15-tl)/16))
			end
		end
	end
end

mutable struct Dimension
	cs::Dict{Pair{Int,Int},Chunk}
end
"画板(0,0)显示坐标(x,y)"
function draw(d::Dimension,x::Float,y::Float)
#=
区块坐标
(x,y)------+---------+
  |    A   |    B    |
  +--(lx+64,ly+64)---+
  |    C   |    D    |
  +--------+---(x+15,y+15)
lx+64-x=64-tx
lx+64<x+15 <=> tx>49 <=> 额外区块需渲染
||显示的区块左上角坐标|取整后|不取整显示坐标/32|显示坐标/32|
|:-:|:-:|:-:|:-:|:-:|
|A|(tx,ty)|(tx1,ty1)|(0,0)|(tx1-tx,ty1-ty)|
|B|(tx,0)|(tx1,0)|(0,64-ty)|(tx1-tx,64-ty)|
|C|(0,ty)|(0,ty1)|(64-tx,0)|(64-tx,ty1-ty)|
|D|(0,0)|(0,0)|(64-tx,64-ty)|(64-tx,64-ty)|
=#
	envl=envlight()
	tx=mod(x,64)
	ty=mod(y,64)
	tx1=floor(Int,tx)
	ty1=floor(Int,ty)
	lx=floor(Int,x)>>5<<5
	ly=floor(Int,y)>>5<<5
	draw(d.cs[Pair(lx,ly)],tx1+1,ty1+1,floor((tx1-tx)*32),floor((ty1-ty)*32);envl=envl) # A
	if ty>49
		draw(d.cs[Pair(lx,ly+64)],tx1+1,1,floor((tx1-tx)*32),floor((64-ty)*32);envl=envl) # B
	end
	if tx>49
		lx+=64
		draw(d.cs[Pair(lx,ly)],1,ty1+1,floor((64-tx)*32),floor((ty1-ty)*32);envl=envl) # C
		if ty>49
			draw(d.cs[Pair(lx,ly+64)],1,1,floor((64-tx)*32),floor((64-ty)*32);envl=envl) # D
		end
	end
end

mutable struct Game
	ds::Dict{String,Dimension}
	ps::Vector{玩家}
	set::Dict{String,Any}
end
function Game()
	return Game(Dict{String,Dimension}(),Vector{玩家}(),Dict{String,Any}())
end
function pushsave(lx::Int,ly::Int,tx::Int,ty::Int,b::Block)
	p=Pair(lx,ly)
	if !haskey(sav,p)
		sav[p]=Dict{Int,Block}()
	end
	sav[p][tx>>6|ty]=b
end
function getlig(x::Int,y::Int)
	tx=x&63;ty=y&63
	lx=x>>6;ly=y>>6
	c=chk(dim,lx,ly)
	return get_lig(c,tx,ty)
end
function setlig(x::Int,y::Int,v::UInt8)
	tx=x&63;ty=y&63
	lx=x>>6;ly=y>>6
	c=chk(dim,lx,ly)
	set_lig(c,tx,ty,v)
end
function getblk(x::Int,y::Int)
	tx=x&63+1;ty=y&63+1
	lx=x>>6;ly=y>>6
	c=chk(dim,lx,ly)
	return c.blks[tx,ty]
end
function setblk(x::Int,y::Int,b::Block;player_put::Bool=false)
	tx=x&63+1;ty=y&63+1
	lx=x>>6;ly=y>>6
	c=chk(dim,lx,ly)
	e_exist(c.blks[tx,ty],x,y)
	c.blks[tx,ty]=b
	if player_put
		e_put(b,x,y)
	else
		exist(b,x,y)
	end
	pushsave(lx,ly,tx,ty,b)
end

function calc_facing(x::Int,y::Int;to::Entity=ply)
	x0=x-to.x;y0=y-to.y
	if x0>0
		if x0<y0 return 0x2
		elseif x0<-y0 return 0x0
		else return 0x1
		end
	else
		if x0>y0 return 0x0
		elseif x0>-y0 return 0x2
		else return 0x3
		end
	end
end


text:
module LegendrePolynomials
using HyperDualNumbers,OffsetArrays

export Pl,Pl!,
dPl,dPl!,
d2Pl,d2Pl!,
Pl_dPl,Pl_dPl!,
Pl_d2Pl,Pl_d2Pl!,
dPl_d2Pl,dPl_d2Pl!,
Pl_dPl_d2Pl, Pl_dPl_d2Pl!

##########################################################################################
# Legendre Polynomials and derivatives up to order 2

@inline function Pl_recursion(ℓ,Plm1,Plm2,x)
	((2*(ℓ-1)+1)* x * Plm1 - (ℓ-1) * Plm2)/ℓ
end

function Pl_hyperdual_allmodes(x::Real,lmax::Integer)
	T = typeof(float(x))
	arr = OffsetArray{Hyper{T}}(undef,0:lmax)
	Pl_hyperdual_allmodes!(arr,x,lmax)
	return arr
end

function Pl_hyperdual_allmodes!(arr::OffsetVector{Hyper{TR},Vector{Hyper{TR}}},x::Real,
	lmax::Integer=maximum(axes(arr,1))) where {TR<:Real}

	# Check if x is within limits
	abs(x) > 1 && throw(DomainError(x,
		"Legendre Polynomials are defined for arguments lying in -1 ⩽ x ⩽ 1"))

	if (lmax≥0)
		arr[0] = hyper(one(x))
	end
	if (lmax≥1)
		hdx = hyper(x,one(x),one(x),zero(x))
		arr[1] = hdx
	end
	@inbounds for ℓ in 2:lmax
		arr[ℓ] = Pl_recursion(ℓ,arr[ℓ-1],arr[ℓ-2],hdx)
	end
	arr
end

# Compute Pl without storing intermediate values
function Pl_hyperdual(x::Real,l::Integer)
	# Check if x is within limits
	abs(x) > 1 && throw(DomainError(x,
		"Legendre Polynomials are defined for arguments lying in -1 ⩽ x ⩽ 1"))

	hd1 = hyper(one(x))
	hdx = hyper(x,one(x),one(x),zero(x))

	if l==0
		return hd1
	elseif l==1
		return hdx
	end

	Plm1 = hdx; Plm2 = hd1
	Pl = Plm1

	for ℓ in 2:l
		Pl = Pl_recursion(ℓ,Plm1,Plm2,hdx)
		Plm2,Plm1 = Plm1,Pl
	end

	Pl
end

function Pl_derivatives_allmodes!(arr::AbstractVector{<:Real},x::Real,
	lmax::Integer=maximum(axes(arr,1)),deriv=0)

	# Only compute Pl
	arr_hyperdual = view(Pl_hyperdual_allmodes(x,lmax),0:lmax)

	@inbounds @. arr[0:lmax] = realpart(arr_hyperdual)
	arr
end

function Pl_derivatives_allmodes!(arr::AbstractMatrix{<:Real},x::Real,
	lmax::Integer=maximum(axes(arr,1)),
	deriv::Integer=maximum(axes(arr,2)))

	arr_hyperdual = view(Pl_hyperdual_allmodes(x,lmax),0:lmax)

	@inbounds for (d_order,f) in zip(0:deriv,(realpart,eps1,eps1eps2))
		@. arr[0:lmax,d_order] = f(arr_hyperdual)
	end
	arr
end

function Pl_derivatives_allmodes(x::Real,lmax::Integer,deriv::Integer=0)
	arr = zeros(0:lmax,0:deriv)
	Pl_derivatives_allmodes!(arr,x,lmax,deriv)
	arr
end

# Convenience functions

"""
	Pl(x::Real;lmax::Integer)

Computes the Legendre Polynomials Pl(x) for the argument x and l=0:lmax.
Returns an OffsetArray P[0:lmax], where P[l] = Pl(x)
"""
Pl(x::Real;lmax::Integer) = Pl_derivatives_allmodes(x,lmax,0)[:,0]

"""
	Pl(x::Real,l::Integer)

Computes the Legendre Polynomials Pl(x) for the argument x and the degree l
"""
function Pl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	realpart(Ph)
end

"""
	dPl(x::Real;lmax::Integer)

Computes the first derivatives of Legendre Polynomials dₓPl(x) for the
argument x and l=0:lmax.
Returns an OffsetArray dP[0:lmax], where `dP[l]` = dₓPl(x)
"""
dPl(x::Real;lmax::Integer) = Pl_derivatives_allmodes(x,lmax,1)[:,1]

"""
	dPl(x::Real,l::Integer)

Computes the first derivatives of Legendre Polynomials dₓPl(x) for the
argument x and the degree l
"""
function dPl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	eps1(Ph)
end

"""
	d2Pl(x::Real;lmax::Integer)

Computes the second derivatives of Legendre Polynomials dₓ²Pl(x) for the
argument x and l=0:lmax
Returns an OffsetArray d2P[0:lmax], where d2P[l] = dₓ²Pl(x)
"""
d2Pl(x::Real;lmax::Integer) = Pl_derivatives_allmodes(x,lmax,2)[:,2]

"""
	d2Pl(x::Real,l::Integer)

Computes the second derivatives of Legendre Polynomials dₓ²Pl(x) for the
argument x and the degree l
"""
function d2Pl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	eps1eps2(Ph)
end

"""
	Pl_dPl(x::Real;lmax::Integer)

Computes the the Legendre Polynomials Pl(x) and their derivatives dₓPl(x)
for the argument x and l=0:lmax
Returns OffsetArrays P[0:lmax] and dP[0:lmax], where P[l] = Pl(x) and
dP[l] = dₓPl(x)
"""
function Pl_dPl(x::Real;lmax::Integer)
	P = Pl_derivatives_allmodes(x,lmax,1)
	P[:,0],P[:,1]
end

"""
	Pl_dPl(x::Real,l::Integer)

Computes the the Legendre Polynomials Pl(x) and their first derivatives
dPl(x) for the argument x and the degree l
"""
function Pl_dPl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	realpart(Ph),eps1(Ph)
end

"""
	Pl_dPl_d2Pl(x::Real;lmax::Integer)

Computes the the Legendre Polynomials Pl(x) and their first derivatives
dPl(x) and second derivatives dₓ²Pl(x) for the argument x and l=0:lmax
Returns OffsetArrays P[0:lmax], dP[0:lmax] and d2P[0:lmax],
where P[l] = Pl(x), dP[l] = dₓPl(x) and d2P[l] = dₓ²Pl(x)
"""
function Pl_dPl_d2Pl(x::Real;lmax::Integer)
	P = Pl_derivatives_allmodes(x,lmax,2)
	P[:,0],P[:,1],P[:,2]
end

"""
	Pl_dPl_d2Pl(x::Real,l::Integer)

Computes the the Legendre Polynomials Pl(x) and their first derivatives
dPl(x) and second derivatives dₓ²Pl(x) for the argument x and the degree l
"""
function Pl_dPl_d2Pl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	realpart(Ph),eps1(Ph),eps1eps2(Ph)
end

"""
	dPl_d2Pl(x::Real;lmax::Integer)

Computes the the first and second derivatives of Legendre Polynomials dPl(x)
and dₓ²Pl(x) for the argument x and l=0:lmax

Returns OffsetArrays dP[0:lmax] and d2P[0:lmax], where dP[l] = dₓPl(x) and
d2P[l] = dₓ²Pl(x)
"""
function dPl_d2Pl(x::Real;lmax::Integer)
	P = Pl_derivatives_allmodes(x,lmax,2)
	P[:,1],P[:,2]
end

"""
	dPl_d2Pl(x::Real,l::Integer)

Computes the the first and second derivatives of Legendre Polynomials dPl(x)
and dₓ²Pl(x) for the argument x and the degree l
"""
function dPl_d2Pl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	eps1(Ph),eps1eps2(Ph)
end

"""
	Pl_d2Pl(x::Real;lmax::Integer)

Computes the the Legendre Polynomials Pl(x) and their second derivatives
dₓ²Pl(x) for the argument x and l=0:lmax

Returns OffsetArrays P[0:lmax] and d2P[0:lmax], where P[l] = Pl(x) and
d2P[l] = dₓ²Pl(x)
"""
function Pl_d2Pl(x::Real;lmax::Integer)
	P = Pl_derivatives_allmodes(x,lmax,2)
	P[:,0],P[:,2]
end

"""
	Pl_d2Pl(x::Real,l::Integer)

Computes the the Legendre Polynomials Pl(x) and their second derivatives
dₓ²Pl(x) for the argument x and the degree l
"""
function Pl_d2Pl(x::Real,l::Integer)
	Ph = Pl_hyperdual(x,l)
	realpart(Ph),eps1eps2(Ph)
end

"""
	Pl!(arr,x;[lmax])

Computes the Legendre Polynomials Pl(x) for the argument x and l=0:lmax,
and saves it in the AbstractArray arr. Assumes that arr has 0-based indexing.

All dimensions of arr should have indices starting from 0
The axes of arr should be (0:l,0:n) where l ⫺ lmax and n ⫺ 0

At output, arr[l,0] = Pl(x)

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
Pl!(arr::AbstractVecOrMat{<:Real},x::Real;lmax::Integer=maximum(axes(arr,1))) =
	Pl_derivatives_allmodes!(arr,x,lmax,0)

"""
	Pl_dPl!(arr,x;[lmax])

Computes the Legendre Polynomials Pl(x) and their derivatives dₓPl(x)
for the argument x and l=0:lmax, and saves them in the AbstractArray arr.
Assumes that arr has 0-based indexing.

All dimensions of arr should have indices starting from 0
The axes of arr should be (0:l,0:n) where l ⫺ lmax and n ⫺ 1

At output, arr[l,0] = Pl(x) and arr[l,1] = dₓPl(x)

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
Pl_dPl!(arr::AbstractMatrix{<:Real},x::Real;lmax::Integer=maximum(axes(arr,1))) =
	Pl_derivatives_allmodes!(arr,x,lmax,min(1,maximum(axes(arr,2))))

"""
	dPl!(arr,x;[lmax])

Computes the first derivatives of Legendre Polynomials dₓPl(x)
for the argument x and l=0:lmax, and saves them in the AbstractArray arr.
Assumes that arr has 0-based indexing.

The first dimension of arr should be 0:l, where l ⫺ lmax
The second dimension -- if arr is 2D -- should contain the index 1

At output, arr[l] = dₓPl(x) if arr is a Vector
or arr[l,1] = dₓPl(x) if arr is a Matrix

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
function dPl!(arr::AbstractVector{<:Real},x::Real;lmax::Integer=maximum(axes(arr,1)))
	P = zeros(axes(arr,1),0:1)
	Pl_dPl_d2Pl!(P,x;lmax=lmax)
	@views @. arr = P[:,1]
	arr
end

function dPl!(arr::AbstractMatrix{<:Real},x::Real;
	lmax::Integer=maximum(axes(arr,1)))
	P = zeros(axes(arr,1),0:1)
	Pl_dPl_d2Pl!(P,x;lmax=lmax)
	@views @. arr[:,1] = P[:,1]
	arr
end

"""
	Pl_dPl_d2Pl!(arr,x;[lmax])

Computes the Legendre Polynomials Pl(x) and their first and second
derivatives dₓPl(x) and dₓ²Pl(x), for the argument x and l=0:lmax,
and saves them in the AbstractArray arr. Assumes that arr has 0-based indexing.

All dimensions of arr should have indices starting from 0
The axes of arr should be (0:l,0:n) where l ⫺ lmax and n ⫺ 2

At output, arr[l,0] = Pl(x), arr[l,1] = dₓPl(x) and arr[l,2] = dₓ²Pl(x)

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
Pl_dPl_d2Pl!(arr::AbstractMatrix{<:Real},x::Real;lmax::Integer=maximum(axes(arr,1))) =
	Pl_derivatives_allmodes!(arr,x,lmax,min(2,maximum(axes(arr,2))))

"""
	d2Pl!(arr,x;[lmax])

Computes the second derivatives of Legendre Polynomials dₓ²Pl(x)
for the argument x and l=0:lmax, and saves them in the AbstractArray arr.
Assumes that arr has 0-based indexing.

The first dimension of arr should be 0:l, where l ⫺ lmax
The second dimension -- if arr is 2D -- should contain the index 2

At output, arr[l] = dₓ²Pl(x) if arr is a Vector
or arr[l,2] = dₓ²Pl(x) if arr is a Matrix

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
function d2Pl!(arr::AbstractVector{<:Real},x::Real;lmax::Integer=maximum(axes(arr,1)))
	P = zeros(axes(arr,1),0:2)
	Pl_dPl_d2Pl!(P,x;lmax=lmax)
	@views @. arr = P[:,2]
	arr
end

function d2Pl!(arr::AbstractMatrix{<:Real},
	x::Real;lmax::Integer=maximum(axes(arr,1)))
	P = zeros(axes(arr,1),0:2)
	Pl_dPl_d2Pl!(P,x;lmax=lmax)
	@views @. arr[:,2] = P[:,2]
	arr
end

"""
	dPl_d2Pl!(arr,x;[lmax])

Computes the first and second derivatives dₓPl(x) and dₓ²Pl(x),
for the argument x and l=0:lmax,
and saves them in the AbstractArray arr. Assumes that arr has 0-based indexing.

The first dimension of arr should be 0:l, where l ⫺ lmax
The second dimension should be (1:n) where n ⫺ 2

At output, arr[l,1] = dₓPl(x) and arr[l,2] = dₓ²Pl(x)

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
function dPl_d2Pl!(arr::AbstractMatrix{<:Real},x::Real;
	lmax::Integer=maximum(axes(arr,1)))

	P = zeros(axes(arr,1),0:2)
	Pl_dPl_d2Pl!(P,x;lmax=lmax)

	@views @. arr[:,1:2] = P[:,1:2]

	arr
end

"""
	Pl_d2Pl!(arr,x;[lmax])

Computes the first and second derivatives dₓPl(x) and dₓ²Pl(x),
for the argument x and l=0:lmax,
and saves them in the AbstractArray arr. Assumes that arr has 0-based indexing.

The first dimension of arr should be 0:l, where l ⫺ lmax
The second dimension should be (0:n) where n ⫺ 2

At output, arr[l,0] = dₓPl(x) and arr[l,2] = dₓ²Pl(x)

The optional keyword argument lmax can specify the range of l's to compute.
It defaults to lmax = maximum(axes(arr,1))
"""
function Pl_d2Pl!(arr::AbstractMatrix{<:Real},x::Real;
	lmax::Integer=maximum(axes(arr,1)))

	P = zeros(axes(arr,1),0:2)
	Pl_dPl_d2Pl!(P,x;lmax=lmax)

	@views @. arr[:,0] = P[:,0]
	@views @. arr[:,2] = P[:,2]

	arr
end

end


text:
include("euler/euler.jl")
using .Numbers: is_palindrome
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(n::Int)::Int
    result = 0
    for i ∈ 10 ^ n ÷ 11 * 11:-11:10 ^ (n - 1)
        for j ∈ 10 ^ n - 1:-1:10 ^ (n - 1)
            candidate = i * j
            if candidate ≤ result
                break
            elseif is_palindrome(candidate)
                result = candidate
            end
        end
    end
    return result
end

compute(2)

compute(3)

@benchmark compute(3)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
#=
If we take 47, reverse and add, 47 + 74 = 121, which is palindromic.

Not all numbers produce palindromes so quickly. For example,

349 + 943 = 1292,
1292 + 2921 = 4213
4213 + 3124 = 7337

That is, 349 took three iterations to arrive at a palindrome.

Although no one has proved it yet, it is thought that some numbers, like 196, never produce a palindrome. A number that never forms a palindrome through the reverse and add process is called a Lychrel number. Due to the theoretical nature of these numbers, and for the purpose of this problem, we shall assume that a number is Lychrel until proven otherwise. In addition you are given that for every number below ten-thousand, it will either (i) become a palindrome in less than fifty iterations, or, (ii) no one, with all the computing power that exists, has managed so far to map it to a palindrome. In fact, 10677 is the first number to be shown to require over fifty iterations before producing a palindrome: 4668731596684224866951378664 (53 iterations, 28-digits).

Surprisingly, there are palindromic numbers that are themselves Lychrel numbers; the first example is 4994.

How many Lychrel numbers are there below ten-thousand?

NOTE: Wording was modified slightly on 24 April 2007 to emphasise the theoretical nature of Lychrel numbers.
=#
function islychrelnumber(n)
  for i in 1:50
    r = reverse(string(n))
    l = n + parse(BigInt, r)
    if string(l) == reverse(string(l))
      return false
    end
    n = l
  end
  true
end

function calc()
  filter((x)->x, [islychrelnumber(x) for x in 1:10_000]) |> length
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test

using Distributions, PDMats, StatsBase

struct test_dist <: Distribution{Univariate, Continuous} end
Distributions.sampler(d::test_dist) = Distributions.sampler(Distributions.Normal(0,1))

struct test_batsampler{T} <: BATSampler{T, Continuous} end

function Base.rand!(rng::AbstractRNG, s::test_batsampler, x::Integer)
    return 0.5
end
Base.length(s::test_batsampler{Multivariate}) = 2
function Base.rand!(rng::AbstractRNG, s::test_batsampler, x::AbstractArray{T, 1} where T)
    for i in indices(x)[1]
        x[i] = 1
    end
    return x
 end

@testset "random number generation" begin
    @testset "bat_sampler" begin
        @test bat_sampler(@inferred test_dist()) == Distributions.Normal(0,1)
    end

    @testset "_check_rand_compat" begin
        @test BAT._check_rand_compat(MvNormal(ones(2)), ones(2,10)) == nothing
        @test_throws DimensionMismatch BAT._check_rand_compat(MvNormal(ones(3)), ones(2,10))
    end

    @testset "rand" begin
        bsguv = BATGammaMTSampler(Gamma(4f0, 2f0))
        @test size(rand(bsguv, 5)) == (5,)
        @test typeof(rand(bsguv, 5)) == Vector{Float32}
        @test typeof(rand(bsguv)) == Float32

        res = rand!(bsguv, zeros(3))
        @test size(res) == (3,)
        @test typeof(res) == Array{Float64, 1}
        res = rand!(bsguv, zeros(3,4))
        @test size(res) == (3,4,)
        @test typeof(res) == Array{Float64, 2}

        @test rand!(test_batsampler{Univariate}(), 1) == 0.5
        x = zeros(2,3)
        rand!(test_batsampler{Multivariate}(), x)
        @test x == ones(2,3)

        res = rand(Base.GLOBAL_RNG, bsguv, Dims((2,3)))
        @test typeof(res) == Array{Float32,2}
        @test size(res) == (2,3)

        bstmv = BATMvTDistSampler(MvTDist(1.5, PDMat(diagm(ones(3)))))
        res = rand(bstmv)
        @test typeof(res) == Array{Float64, 1}
        @test size(res) == (3,)
        res = rand(bstmv, Int32(2))
        @test typeof(res) == Array{Float64, 2}
        @test size(res) == (3, 2)
    end

    @testset "issymmetric_around_origin" begin
        @test issymmetric_around_origin(Normal(0.0, 3.2)) == true
        @test issymmetric_around_origin(Normal(1.7, 3.2)) == false
        @test issymmetric_around_origin(Gamma(4.2, 2.2)) == false
        @test issymmetric_around_origin(Chisq(20.3)) == false
        @test issymmetric_around_origin(MvNormal(zeros(2), ones(2))) == true
        @test issymmetric_around_origin(MvNormal(ones(2), ones(2))) == false
        @test issymmetric_around_origin(MvTDist(1.5, zeros(2), PDMat(diagm(ones(2)))))
        @test issymmetric_around_origin(MvTDist(1.5, ones(2), PDMat(diagm(ones(2))))) == false
    end

    @testset "BATGammaMTSampler" begin
        @test typeof(@inferred bat_sampler(Gamma(0.5, 20.3))) <: BATGammaMTSampler

        @test rand(MersenneTwister(7002), BATGammaMTSampler(Gamma(0.5, 20.3))) ≈ 0.872195443488103
        @test rand(MersenneTwister(8456), BATGammaMTSampler(Gamma(0.5, 20.3))) ≈ 0.2723580272080763
        @test rand(MersenneTwister(3296), BATGammaMTSampler(Gamma(0.5, 20.3))) ≈ 6.984079300633383
        @test rand(MersenneTwister(5065), BATGammaMTSampler(Gamma(0.5, 20.3))) ≈ 9.387736826614171
        @test rand(MersenneTwister(7002), BATGammaMTSampler(Gamma(0.25, 2))) ≈ 0.0031794482969992074
        @test rand(MersenneTwister(7002), BATGammaMTSampler(Gamma(10, 2))) ≈ 4.923284832865044

        @test typeof(rand(BATGammaMTSampler(Gamma(4.2f0, 2.2f0)))) == Float32
        @test eltype(BATGammaMTSampler(Gamma(4.2f0, 2.2f0))) == Float32
        @test eltype(BATGammaMTSampler(Gamma(4.2, 2.2))) == Float64
        @test typeof(rand(BATGammaMTSampler(Gamma(4.2, 2.2)))) == Float64
    end

    @testset "BATChisqSampler" begin
        @test typeof(@inferred bat_sampler(Chisq(0.5))) <: BATChisqSampler

        @test rand(MersenneTwister(7002), BATChisqSampler(Chisq(0.5))) ≈ 0.0031794482969992074
        @test rand(MersenneTwister(7002), BATChisqSampler(Chisq(20))) ≈ 4.923284832865044

        @test eltype(BATChisqSampler(Chisq(20f0))) == Float32
        @test typeof(rand(BATChisqSampler(Chisq(20f0)))) == Float32
    end

    @testset "BATMvTDistSampler" begin
        d = MvTDist(1.5, PDMat([2.0 1.0; 1.0 3.0]))

        @test typeof(@inferred bat_sampler(d)) <: BATMvTDistSampler
        @test size(@inferred rand(bat_sampler(d), 5)) == (2, 5)

        @test size(rand!(bat_sampler(d), ones(2))) == (2, )
        @test size(rand!(bat_sampler(d), ones(2, 3))) == (2, 3)

        cmat = [3.76748 0.446731 0.625418; 0.446731 3.9317 0.237361; 0.625418 0.237361 3.43867]
        tmean = [1., 2, 3]
        tmv = MvTDist(3, tmean, PDMat(eye(3)))

        tmv2 = BAT.set_cov(tmv, cmat)
        @test full(BAT.get_cov(tmv2)) ≈ cmat

        bstmv = BATMvTDistSampler(tmv2)

        n = 1000
        res = rand(MersenneTwister(7002), bstmv, n)

        @test isapprox(mean(res,2), tmean; atol = 0.5)
        @test isapprox(cov(res, 2)/3, cmat; atol = 1.5)
    end
end


text:
module JLL

using GAP
using GAP_pkg_io_jll

function __init__()
    # ensure GAP finds kernel extensions or other binaries
    sopath = joinpath(GAP_pkg_io_jll.find_artifact_dir(), "lib", "gap")
    @debug "GAP package 'io' sopath = " * sopath
    setproperty!(GAP.Globals.DirectoriesPackageProgramsOverrides, :io, GapObj(sopath))
end

end # module


text:
import SIUnits
import SIUnits.ShortUnits: s
typealias SecondT{T} SIUnits.SIQuantity{T,0,0,1,0,0,0,0,0,0}

# A signal has a common timebase (an abstract vector of type T) in seconds,
# and an array of data that it iterates and indexes over (each element is an
# abstract vector of type S). The data are stored in a collection of type D,
# which is an array of rank N+1.
abstract AbstractSignal{T<:AbstractVector, S<:AbstractVector, N} <: AbstractArray{S, N}
type Signal{T, S, N, D} <: AbstractSignal{T, S, N}
    time::T
    data::D

    dims::NTuple{N, Symbol}

    meta::Dict{Symbol, Any}
    Signal(t,d,n,m) = (_checkargs(t,d); new(t,d,n,m))
end
# There are some gotchas when using typealiases in parametric dispatch, so
# the parameters are re-arranged such that the restricted parameter comes last.
# See http://github.com/JuliaLang/julia/issues/7453
typealias RegularSignal{S, N, D, T<:Range} Signal{T, S, N, D}
typealias SignalVector{T, S, D} Signal{T, S, 1, D}
typealias SignalMatrix{T, S, D} Signal{T, S, 2, D}

# Consistency checks that aren't enforced by the type itself
_checktime{T<:SecondT}(t::Range{T}) = step(t) > 0s || throw(ArgumentError("time vector must be monotonically increasing"))
_checktime{T<:SecondT}(t::AbstractVector{T}) = issorted(t,lt=(<=)) || throw(ArgumentError("time vector must be monotonically increasing"))
_checktime(t) = throw(ArgumentError("time must be a vector specified in seconds"))
function _checkargs(t, d)
    _checktime(t)
    size(d, 1) == length(t) || throw(ArgumentError("each channel must be the same length as time"))
end

# The canonical constructor. This must figure out what the element type will be.
# This is really hard with SubArrays and slices of Signals.
stagedfunction Signal{T<:SecondT, R, N}(time::AbstractVector{T}, data::AbstractArray{R, N}, dims::(Symbol...), meta::Dict{Symbol, Any})
    S = SubArray{R, 1, data, (Colon, Int), 2}
    # If the data are Signals themselves, we'll return a SignalVector
    if data <: Signal
        S = Signal{time, R, 1, S} # TODO: is S right here?
    end
    quote
        Signal{$time, $S, N-1, $data}(time, data, dims, meta)
    end
end

# A more forgiving constructor with defaults
function Signal{T,N}(time::AbstractVector, data::AbstractArray{T,N} = Array(Any, length(time), 0),
                     dims::(Symbol...) = ntuple(N-1, (i)->symbol("")),
                     meta::Dict{Symbol,Any} = Dict{Symbol, Any}())
    Signal(inseconds(time), data, dims, meta)
end
# For simple testing, allow vectorized functions. This could be improved
Signal(time::AbstractVector, fcns::Vector{Function}) = Signal(time, Float64[f(i) for i in time, f in fcns])

# Convert the time ranges and vectors to Seconds without re-allocating them
# TODO: There has got to be a better way! github.com/Keno/SIUnits.jl/issues/25
inseconds{T<:Real}(time::Range{T})  = SIUnits.SIRange{typeof(time),T,0,0,1,0,0,0,0,0,0}(time)
inseconds{T<:Real}(time::Vector{T}) = convert(Vector{SecondT{T}}, time)
inseconds{T<:Real}(time::Vector{SIUnits.SIQuantity{T}}) = convert(Vector{SecondT{T}}, time)
inseconds{T<:SecondT}(time::AbstractVector{T}) = time
inseconds{R,T}(time::SIUnits.SIRange{R,T,0,0,1,0,0,0,0,0,0}) = time
inseconds(time::AbstractVector) = throw(ArgumentError("unsupported time vector type"))

# Convert to a RegularSignal by blindly shifting time underneath the data
# TODO: perhaps I should check the variance of diff(time(sig))?
regularize(sig::RegularSignal) = sig
regularize(sig::Signal) = Signal(inseconds(linrange(float(sig.time[1]),float(sig.time[end]),length(sig.time))), sig.data, sig.dims, sig.meta)

# Test if a signal is "regular" -- that is, is it sampled at an exact interval?
isregular(::Signal) = false
isregular(::RegularSignal) = true

## Index and iterate over the first dimension of the data
Base.size(sig::Signal) = size(sig.data)[2:end]
Base.size(sig::Signal, I::Int) = size(sig.data, I+1)
Base.elsize(sig::Signal) = size(sig.time)
# TODO: Perhaps allow custom names/indexes like DataFrames?
Base.getindex(sig::AbstractSignal, idxs::Union(Colon,Int,Array{Int,1},Range{Int})...) = sub(sig, idxs...)

# TODO: We need to propogate dimension names, maybe metadata, too?
Base.sub(sig::Signal, idx::Int) = (checkbounds(sig, idx); sub(sig.data, :, idx))
Base.sub(sig::Signal, idxs::Int...) = (checkbounds(sig, idxs...); sub(sig.data, :, idxs...))
Base.sub(sig::Signal, idxs::Union(Colon,Int,Array{Int,1},Range{Int})...) = (checkbounds(sig, idxs...); Signal(sig.time, sub(sig.data, :, idxs...)))

# Use the default fast linear indexing for iteration
Base.linearindexing(::AbstractSignal) = Base.LinearFast()

# TODO: remove this once https://github.com/JuliaLang/julia/pull/10133 is merged
# (and propagates through the nightlies and Travis)
import Base: checkbounds, trailingsize
checkbounds(sz::Int, ::Colon) = nothing
function checkbounds(A::AbstractArray, I::Union(Real,Colon,AbstractArray)...)
    n = length(I)
    if n > 0
        for dim = 1:(n-1)
            checkbounds(size(A,dim), I[dim])
        end
        checkbounds(trailingsize(A,n), I[n])
    end
end

Base.reshape(sig::Signal, idxs::(Int64...,)) = Signal(sig.time, reshape(sig.data, tuple(length(sig.time), idxs...)))

# Information specific to regular signals:
# Hack to get around poor typing in SIUnits math. I know that these are s⁻¹ & s.
# TODO: should I give up trying to store time as an SI array? It can be a pain.
samplingfreq(sig::RegularSignal) = (v = float(1/step(sig.time)); SIUnits.SIQuantity{typeof(v),0,0,-1,0,0,0,0,0,0}(v))
samplingrate(sig::RegularSignal) = (v = float(step(sig.time));   SIUnits.SIQuantity{typeof(v),0,0,1,0,0,0,0,0,0}(v))


text:
export EntropyConeLift, equalonsubsetsof!, equalvariable!

mutable struct EntropyConeLift{T<:Real} <: AbstractEntropyCone{T}
    n::Vector{Int}
    poly::Polyhedron{T}

    function EntropyConeLift{T}(n::Vector{Int}, poly::Polyhedron{T}) where {T}
        new{T}(n, poly)
    end

    function EntropyConeLift{T}(n::Vector{Int}, A::AbstractMatrix{T}, equalities::BitSet) where {T}
        if sum(ntodim(n)) != size(A, 2)
            error("The dimensions in n does not agree with the number of columns of A")
        end
        if !isempty(equalities) && last(equalities) > size(A, 1)
            error("Equalities should range from 1 to the number of rows of A")
        end
        ine = hrep(A, spzeros(T, size(A, 1)), equalities)
        new{T}(n, polyhedron(ine))
    end

end

Polyhedra.fulldim(h::EntropyConeLift) = Int(sum(ntodim(h.n)))

EntropyConeLift(n::Vector{Int}, A::AbstractMatrix{T}) where {T<:Real} = EntropyConeLift(n, A, BitSet())

Base.convert(::Type{EntropyConeLift{S}}, H::EntropyConeLift{T}) where {S<:Real,T<:Real} = EntropyConeLift{S}(H.n, Polyhedron{S}(H.poly))

Base.convert(::Type{EntropyConeLift{T}}, h::EntropyCone{T}) where {T<:Real} = EntropyConeLift([h.n], h.poly)

#Base.getindex{T<:Real}(H::EntropyConeLift{T}, i) = DualEntropyLift(H.n, H.A[i,:], i in H.equalities)

function offsetfor(h::EntropyConeLift, id::Integer)
    id == 1 ? 0 : sum(map(ntodim, h.n[1:(id-1)]))
end
function indset(h::AbstractEntropyCone, id::Integer)
    indset(h.n[id])
end
function rangefor(h::EntropyConeLift, id::Integer)
    offset = offsetfor(h, id)
    offset .+ indset(h, id)
end

promote_rule(::Type{EntropyConeLift{T}}, ::Type{EntropyCone{T}}) where {T<:Real} = EntropyConeLift{T}

function (*)(x::AbstractEntropyCone{T}, y::AbstractEntropyCone{T}) where {T<:Real}
    # A = [x.A zeros(T, size(x.A, 1), N2); zeros(T, size(y.A, 1), N1) y.A]
    # equalities = copy(x.equalities)
    # for eq in y.equalities
    #   push!(equalities, size(x.A, 1) + eq)
    # end
    EntropyConeLift{T}([x.n; y.n], x.poly * y.poly)
end

function equalonsubsetsof!(H::EntropyConeLift{T}, id1, id2, S::EntropyIndex, I::EntropyIndex=emptyset(), σ=collect(1:H.n[id1])) where {T}
    if S == emptyset()
        return
    end
    nrows = (1<<(card(S)))-1
    A = spzeros(T, nrows, Polyhedra.fulldim(H))
    cur = 1
    offset1 = offsetfor(H, id1)
    offset2 = offsetfor(H, id2)
    for K in setsto(S)
        if K ⊆ S && !(K ⊆ I)
            A[cur, offset1+K] = 1
            A[cur, offset2+mymap(σ, K, H.n[id2])] = -1
            cur += 1
        end
    end
    ine = MixedMatHRep(A, spzeros(T, nrows), BitSet(1:nrows))
    intersect!(H, ine)
end
equalonsubsetsof!(H::EntropyConeLift, id1, id2, s::Signed) = equalonsubsetsof!(H, id1, id2, set(s))

function equalvariable!(h::EntropyConeLift{T}, id::Integer, i::Signed, j::Signed) where {T}
    if id < 1 || id > length(h.n) || min(i,j) < 1 || max(i,j) > h.n[id]
        error("invalid")
    end
    if i == j
        @warn "useless"
        return
    end
    nrows = 1 << (h.n[id]-1)
    A = spzeros(T, nrows, Polyhedra.fulldim(h))
    offset = offsetfor(h, id)
    cur = 1
    for S in indset(h, id)
        if myin(i, S)
            A[cur, offset+S] = 1
            Q = union(setdiff(S, set(i)), set(j))
            A[cur, offset+Q] = -1
            cur += 1
        end
    end
    intersect!(h, MixedMatHRep(A, spzeros(T, nrows), BitSet(1:nrows)))
end

ninneradh(n, J::EntropyIndex, K::EntropyIndex) = n
nselfadh(n, J::EntropyIndex, I::EntropyIndex) = n + card(setdiff(J, I))
nadh(n, J::EntropyIndex, K::EntropyIndex, adh::Type{Val{:Inner}}) = ninneradh(n, J, K)
nadh(n, J::EntropyIndex, K::EntropyIndex, adh::Type{Val{:Self}}) = nselfadh(n, J, K)
nadh(n, J::EntropyIndex, K::EntropyIndex, adh::Type{Val{:NoAdh}}) = n
nadh(n, J::EntropyIndex, K::EntropyIndex, adh::Symbol) = nadh(n, J, K, Val{adh})

function inneradhesivelift(h::EntropyCone{T}, J::EntropyIndex, K::EntropyIndex) where {T}
    cur = polymatroidcone(T, ninneradh(h.n, J, K))
    intersect!(cur, submodulareq(cur.n, J, K))
    lift = h * cur
    I = J ∩ K
    equalonsubsetsof!(lift, 1, 2, J)
    equalonsubsetsof!(lift, 1, 2, K, I)
    lift
end
function selfadhesivelift(h::EntropyCone{T}, J::EntropyIndex, I::EntropyIndex) where {T}
    newn = nselfadh(h.n, J, I)
    K = setdiff(fullset(newn), fullset(h.n)) ∪ I
    cur = polymatroidcone(T, newn)
    intersect!(cur, submodulareq(cur.n, fullset(h.n), K, I))
    lift = h * cur
    equalonsubsetsof!(lift, 1, 2, fullset(h.n))
    themap = Vector{Int}(undef, h.n)
    cur = h.n
    for i in 1:h.n
        if myin(i, I)
            themap[i] = i
        elseif myin(i, J)
            cur += 1
            themap[i] = cur
        else
            themap[i] = -1
        end
    end
    @assert cur == newn
    equalonsubsetsof!(lift, 1, 2, J, I, themap)
    lift
end
adhesivelift(h::EntropyCone, J::EntropyIndex, K::EntropyIndex, adh::Type{Val{:Inner}}) = inneradhesivelift(h, J, K)
adhesivelift(h::EntropyCone, J::EntropyIndex, K::EntropyIndex, adh::Type{Val{:Self}}) = selfadhesivelift(h, J, K)
adhesivelift(h::EntropyCone, J::EntropyIndex, K::EntropyIndex, adh::Symbol) = adhesivelift(h, J, K, Val{adh})


text:
# a simple benchmark

using Clustering

m = 100
n = 10000
k = 50

x = rand(m, n)

# warming
kmeans(x[:,1:1000], 2; display=:none, max_iter=200)

println("kmeans on $n samples (of dimension $m) with K = $k ...")
t_new = @elapsed r_new = kmeans(x, k; display=:none, max_iter=200)
println("\telapsed = $t_new, per_iteration = $(t_new / r_new.iterations)")
println()


text:
# Collection of utility functions for dealing with arrays

"""
Convert a matrix of real values into complex values by treating alternating
columns as real and imaginary parts
"""
function real_cols_to_complex(mat::Matrix{T}) where {T}
    ncols = size(mat,2)
    out = Matrix{Complex{T}}(undef, size(mat,1), div(ncols,2))
    for col in 1:div(ncols,2)
        out[:,col] .= mat[:,2col-1] .+ mat[:,2col]*im
    end
    return out
end

"Turns vector of vectors into vector by concatenating them together"
function flatten(arr::Vector{Vector{T}}) where {T}
    flattened_arr = Vector{T}(undef, 0)
    for i = 1:length(arr)
        append!(flattened_arr, arr[i])
    end
    return flattened_arr
end

"Turns n-dim array into vector"
flatten(a::AbstractArray) = reshape(a, :)

function rel_abs_error(exact::Complex{T}, vec::Complex{T})::T where {T}
    abs(exact - vec)/abs(exact)
end

function max_rel_abs_error(exact::Vector{Complex{T}}, vec::Vector{Complex{T}})::T where {T}
    current_max = zero(T)
    for i in eachindex(exact)
        current_max = max(current_max, rel_abs_error(exact[i], vec[i]))
    end
    return current_max
end

"""
    offsetmap(fnc::Function, inds::StepRange)::OffsetArray

Maps `fnc` over `inds`, returning an OffsetArray with indices of `inds`.
"""
offsetmap(fnc, inds) = OffsetArray(map(fnc, inds), inds)


text:
function process_pokes(ongoing::AbstractDict)
    curr_data = adjust_matfile(ongoing)
    dropmissing!(curr_data, :Trial)
    process_pokes(curr_data)
end

function process_pokes(curr_data::AbstractDataFrame)
    curr_data[!,:Poke_within_Trial] = Vector{Union{Float64,Missing}}(undef,size(curr_data,1))
    curr_data[!,:Pre_Interpoke] = Vector{Union{Float64,Missing}}(undef,size(curr_data,1))
    curr_data[!,:Post_Interpoke] = Vector{Union{Float64,Missing}}(undef,size(curr_data,1))
    curr_data[!,:LastPoke] .= false
    combine(groupby(curr_data,:Trial)) do dd
        dd[:,:Poke_within_Trial] = collect(1:nrow(dd))
        dd[:,:Pre_Interpoke] =  dd[!,:PokeIn] .-lag(dd[!,:PokeOut],default = missing)
        dd[:,:Post_Interpoke] = lead(dd[!,:PokeIn],default = missing).- dd[!,:PokeOut]
        dd[end,:LastPoke] = true
    end
    curr_data.Poke_Hierarchy = Vector{Union{Float64,Missing}}(undef,nrow(curr_data))
    combine(groupby(curr_data,:Trial)) do dd
         dd[:,:Poke_Hierarchy] = ismissing(dd[1,:Reward]) ? [missing] : Switching.get_hierarchy(dd.Reward)
    end
    curr_data[!,:ReverseTrial] .= reverse(curr_data[:,:Trial])
    curr_data[!,:Correct] = [ismissing(x) ? false : true for x in curr_data[!,:Poke]]
    curr_data[!,:PokeDuration] = curr_data.PokeOut - curr_data.PokeIn
    curr_data[!,:Stim_Day] .= length(union(curr_data.Stim)) > 1
    # filter!(r -> r.PokeDuration > 0.10,curr_data)
    return curr_data
end

function adjust_matfile(ongoing::AbstractDict)
    # switching_type = match(r"switching\d{1}",filename).match
    # extract outcomes and protocols of both sides in 2 simpler variable to be later pointed
    outcomesR = ongoing["saved_history"]["switching2_RrewList"]
    outcomesL = ongoing["saved_history"]["switching2_LrewList"]
    protocolsR = ongoing["saved_history"]["switching2_current_rightMax"]
    protocolsL = ongoing["saved_history"]["switching2_current_leftMax"]
    Stimulation = ongoing["saved_history"]["LaserSection_currentStim"]
    #this is what contains the series of pokes
    trials_vec = ongoing["saved_history"]["ProtocolsSection_parsed_events"]
    poke_count = 0
    falsetrial_count = 0
    #create an empty dataframe wher to push single poke rows info
    df = DataFrame(Poke = Union{Int64,Missing}[], Trial = Union{Int64,Missing}[],
        PokeIn = Union{Float64,Missing}[], PokeOut = Union{Float64,Missing}[],
        Side = Union{String,Missing}[], Reward = Union{Bool,Missing}[],
        ROI_In = Float64[], ROI_Out = Float64[], Protocol = Union{String,Missing}[],
        Stim = Union{Bool,Missing}[])

    #loops over the trials and index them
    for (idx,t) in enumerate(trials_vec)
        # retrieve safe alternative for starting and ending time of pokes when value is missing
        default_in = t["states"]["state_0"][1,2]
        default_out = t["states"]["state_0"][2,1]
        # retrieve the entry and exit time to the ROI
        ROI_activity = t["pokes"]["C"]
        if size(ROI_activity,1)  == 0
            ROI_activity_in = default_in
            ROI_activity_out = default_out
        else
            ROI_activity_in = isnan(ROI_activity[1]) ? default_in : ROI_activity[1]
            ROI_activity_out = isnan(ROI_activity[end]) ? default_out : ROI_activity[end]
        end
        # define in which side was the trial, if no poke was made returns nothing
        if length(t["pokes"]["R"]) > 0
            p = t["pokes"]["R"]
            side = "R"
            outcomes = outcomesR[idx]
            protocol = string(protocolsR[idx])
            stim = Bool(Stimulation[idx])
        elseif length(t["pokes"]["L"]) > 0
            p = t["pokes"]["L"]
            side = "L"
            outcomes = outcomesL[idx]
            protocol = string(protocolsL[idx])
            stim = Bool(Stimulation[idx])
        else
            p = nothing
        end
        # if poke isn't nothing loops over the pokes to push all the info in 1 row
        if !isnothing(p)
            for i  = 1:size(p,1)
                poke_count += 1
                rew = i>20 ? false : Bool(outcomes[i])
                p[1] = isnan(p[1]) ? default_in : p[1]
                p[end] = isnan(p[end]) ? default_out : p[end]
                push!(df,[poke_count,idx-falsetrial_count,p[i,1],p[i,2],side,rew,ROI_activity_in,ROI_activity_out,protocol,stim])
            end
        else
            falsetrial_count += 1
            push!(df,[missing,missing,missing,missing,missing,missing,ROI_activity_in,ROI_activity_out,missing,missing])
        end

    end
    #test = size(df,1) == length(trials_vec)
    return df
end


text:
using Franklin, Test, Markdown, Dates, Random, Literate
const F = Franklin
const R = @__DIR__
const D = joinpath(dirname(dirname(pathof(Franklin))), "test", "_dummies")

F.FD_ENV[:SILENT_MODE] = true
# F.FD_ENV[:DEBUG_MODE] = true
F.FD_ENV[:STRUCTURE] = v"0.1" # legacy, it's switched up in the tests.

# UTILS
println("UTILS-1")
include("utils/folder_structure.jl")
include("utils/paths_vars.jl"); include("test_utils.jl")
# ---
println("UTILS-2")
include("utils/misc.jl")
include("utils/errors.jl")
include("utils/html.jl")
include("regexes.jl")
println("🍺")

# MANAGER folder
println("MANAGER")
include("manager/utils.jl")
include("manager/utils_fs2.jl")
include("manager/rss.jl")
include("manager/config.jl")
include("manager/config_fs2.jl")
include("manager/dir_utils.jl")
include("manager/page_vars_html.jl")
println("🍺")

# PARSER folder
println("PARSER/MD+LX")
include("parser/1-tokenize.jl")
include("parser/2-blocks.jl")
include("parser/markdown+latex.jl")
include("parser/markdown-extra.jl")
include("parser/footnotes+links.jl")
include("parser/latex++.jl")
include("parser/indentation++.jl")
include("parser/md-dbb.jl")
println("🍺")

# EVAL
println("EVAL")
include("eval/module.jl")
include("eval/run.jl")
include("eval/io.jl")
include("eval/io_fs2.jl")
include("eval/codeblock.jl")
include("eval/eval.jl")
include("eval/eval_fs2.jl")
include("eval/integration.jl")
include("eval/extras.jl")

# CONVERTER folder
println("CONVERTER/MD")
include("converter/md/markdown.jl")
include("converter/md/markdown2.jl")
include("converter/md/markdown3.jl")
include("converter/md/markdown4.jl")
include("converter/md/hyperref.jl")
include("converter/md/md_defs.jl")
include("converter/md/md_defs2.jl")
include("converter/md/tags.jl")
println("🍺")
println("CONVERTER/HTML")
include("converter/html/html.jl")
include("converter/html/html2.jl")
include("converter/html/html_for.jl")
println("🍺")
println("CONVERTER/LX")
include("converter/lx/input.jl")
include("converter/lx/simple.jl")
include("converter/lx/simple_fs2.jl")
println("🍺")

fs2()

println("GLOBAL")
include("global/cases1.jl")
include("global/cases2.jl")
include("global/ordering.jl")
include("global/html_esc.jl")

begin
    # create temp dir to do complete integration testing (has to be here in
    # order to locally play nice with node variables etc, otherwise it's a big
    # headache)
    p = joinpath(D, "..", "__tmp");
    # after errors, this may not have been deleted properly
    isdir(p) && rm(p; recursive=true, force=true)
    # make dir, go in it, do the tests, then get completely out (otherwise
    # windows can't delete the folder)
    mkdir(p); cd(p);
    include("global/postprocess.jl");
    include("global/rss.jl")
    cd(p)
    include("global/eval.jl")
    cd(joinpath(D, ".."))
    # clean up
    rm(p; recursive=true, force=true)
end
cd(dirname(dirname(pathof(Franklin))))

println("TEMPLATING")
include("templating/for.jl")
include("templating/fill.jl")

println("UTILS FILE")
include("utils_file/basic.jl")

println("INTEGRATION")
include("integration/literate.jl")
include("integration/literate_fs2.jl")
include("integration/literate_extras.jl")
include("integration/hfuns.jl")

flush_td()
cd(joinpath(dirname(dirname(pathof(Franklin)))))

println("HTML validation")
include("html/closep.jl")
include("html/closep_lx.jl")
include("html/pdiv.jl")

println("COVERAGE")
include("coverage/extras1.jl")
include("coverage/paths.jl")

println("😅 😅 😅 😅")

# check quickly if the IPs in IP_CHECK are still ok
println("Verifying ip addresses, if online these should succeed.")
for (addr, name) in F.IP_CHECK
    println(rpad("Ping $name:", 13), ifelse(F.check_ping(addr), "✓", "✗"), ".")
end


text:
struct SimulatedAnnealing{Tn, Ttemp} <: RandomSearch
    neighbor::Tn
    temperature::Ttemp
    n_per_temperature::Int
end
Base.summary(::SimulatedAnnealing) = "Simulated Annealing"

"""
# SimulatedAnnealing
## Constructor
```julia
SimulatedAnnealing(; neighbor = default_neighbor, temperature = log_temperature)
```

The constructor takes two keywords:
* `neighbor = a(x_now, [x_candidate])`, a function that generates a new iterate based on the current
* `temperature = b(iteration)`, a function of the current iteration that returns a temperature

## Description
Simulated Annealing is a derivative free method for optimization. It is based on the
Metropolis-Hastings algorithm that was originally used to generate samples from a
thermodynamics system and is often used to generate draws from a posterior when doing
Bayesian inference. As such, it is a probabilistic method for finding the minimum of a
function, often over a quite large domains. For the historical reasons given above, the
algorithm uses terms such as cooling, temperature, and acceptance probabilities.
"""
SimulatedAnnealing(;
  neighbor = default_neighbor,
  temperature = log_temperature,
  n_per_temperature = 1) =
    SimulatedAnnealing(neighbor, temperature, n_per_temperature)

log_temperature(t) = 1 / log(t)^2

function default_neighbor(x_best)
  T = eltype(x_best)
  n = length(x_best)
  return x_best .+ T.(randn(n))
end

function solve(prob::OptimizationProblem, x0, method::SimulatedAnnealing, options::OptimizationOptions)
  T = eltype(x0)
  t0 = time()

  x_best = copy(x0)
  f_best = value(prob, x_best)
  x_now = copy(x0)
  f_now = f_best
  f0 = f_now
  temperature = f_best
  iter = 0
  is_converged = converged(method, f_now, options)
  while iter ≤ options.maxiter && !(is_converged)
    iter += 1
    # Determine the temperature for current iteration
    temperature = method.temperature(iter)

    for n_temp = 1:method.n_per_temperature
      # Randomly generate a neighbor of our current state
      x_candidate = method.neighbor(x_best)

      # Evaluate the cost function at the proposed state
      f_candidate = value(prob, x_candidate)

      if f_candidate <= f_now # this handles non-finite values as well
        # If proposal is superior, we always move to it
        x_now = copy(x_candidate)
        f_now = f_candidate

        # If the new state is the best state yet, keep a record of it
        if f_candidate < f_best
          x_best = copy(x_now)
          f_best = f_now
        end
      else
        # If proposal is inferior, we move to it with probability p
        p = exp(-(f_candidate - f_now) / temperature)
        if rand() <= p
          x_now = copy(x_candidate)
          f_now = f_candidate
        end
      end
    end # n_per_temperature
    is_converged = converged(method, f_now, options)
  end

  ConvergenceInfo(method, (solution=x_best, minimum=f_best, f_now=f_now, x_now=x_now, temperature=temperature, f0=f0, iter=iter, time=time()-t0), options)
end
function converged(method::SimulatedAnnealing, fz, options)
  f_converged = false
  f_converged = f_converged || fz ≤ options.f_limit
  f_converged = f_converged || isnan(fz)
end


text:
export ExactPermutationTest, ApproximatePermutationTest

function ptstats(x,y,f)
    xy = vcat(x,y)
    rx = 1:length(x)
    ry = (length(xy)-length(y)+1):length(xy)
    (xy, rx, ry)
end

struct PermutationTest{T<:Real} <: HypothesisTest
    observation::T
    samples::Vector{T}
end

"""
    ExactPermutationTest(x::Vector, y::Vector, f::Function)

Perform a permutation test (a.k.a. randomization test) of the null hypothesis
that `f(x)` is equal to `f(y)`.  All possible permutations are sampled.
"""
function ExactPermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
                              f::Function) where {R<:Real,S<:Real}
    xy, rx, ry = ptstats(x,y,f)
    P = permutations(xy)
    samples = [f(view(p,rx)) - f(view(p,ry)) for p in P]
    PermutationTest(f(x) - f(y), samples)
end

"""
    ApproximatePermutationTest(x::Vector, y::Vector, f::Function, n::Int)

Perform a permutation test (a.k.a. randomization test) of the null hypothesis
that `f(x)` is equal to `f(y)`.  `n` of the `factorial(length(x)+length(y))`
permutations are sampled at random.
"""
function ApproximatePermutationTest(x::AbstractVector{R}, y::AbstractVector{S},
                                    f::Function, n::Int) where {R<:Real,S<:Real}
    xy, rx, ry = ptstats(x,y,f)
    samples = [(shuffle!(xy); f(view(xy,rx)) - f(view(xy,ry))) for i = 1:n]
    PermutationTest(f(x) - f(y), samples)
end

function pvalue(apt::PermutationTest; tail=:both)
    if tail == :both
        count = sum(abs(apt.observation) >= abs(x) for x in apt.samples)
    elseif tail == :left
        count = sum(apt.observation >= x for x in apt.samples)
    elseif tail == :right
        count = sum(apt.observation <= x for x in apt.samples)
    end
    return count / length(apt.samples)
end

testname(::PermutationTest) = "Permutation Test"

function show_params(io::IO, apt::PermutationTest, ident)
    println(io, ident, "observation: ", apt.observation)
    println(io, ident, "samples: ", apt.samples)
end


text:
using ..JuliaRL.FeatureCreators
import ..build_algorithm_dict
import Flux: Descent

"""
    Agent for learning online in the four rooms environment.
"""
mutable struct IntFourRoomsContAgent{O, M<:Resampling.AbstractPolicy, P<:Resampling.AbstractPolicy} <: JuliaRL.AbstractAgent
    μ::M
    sample_policy::P
    WER::WeightedExperienceReplay
    opt::O
    tilecoder::TileCoder
    value_arr::Array{Resampling.SparseLayer, 1}
    lu::Resampling.LearningUpdate
    gvf::GVF
    training_gap::Int64
    buffer_size::Int64
    batch_size::Int64
    warm_up::Int64
    counter::Int64
    state_t::Tuple{Array{Float64, 1}, Bool}
    action_t::Int64
    α_arr::Array{Float64, 1}
    max_is::Float64
    function IntFourRoomsContAgent(μ::M, sample_policy::P, gvf::GVF, opt::O, num_tilings, num_tiles, α_arr::Array{Float64, 1},
                                   training_gap, buffer_size, batch_size,
                                   warm_up, env_size; max_is_ratio=1.0) where {O, M<:Resampling.AbstractPolicy, P<:Resampling.AbstractPolicy}

        # build tile coder....
        tilecoder = TileCoder(num_tilings, num_tiles, 2)
        num_features = num_tilings*(num_tiles+1)^2
        value_arr = [SparseLayer(num_features, 1) for α in α_arr]
        er_names = [:s_t, :s_tp1, :a_t, :a_tp1, :r, :γ_tp1, :terminal, :mu_t, :pi_t, :ρ]
        state_type = Array{Int64, 1}
        er_types = [state_type, state_type, Int64, Int64, Float64, Float64, Bool, Float64, Float64, Float64]
        WER = WeightedExperienceReplay(buffer_size, er_types; column_names=er_names)
        new{O, M, P}(μ, sample_policy, WER, opt, tilecoder, value_arr, Resampling.BatchTD(),
                     gvf, training_gap, buffer_size, batch_size, warm_up,
                     warm_up, ([0.0,0.0], false), 0, α_arr, max_is_ratio)
    end
end

JuliaRL.get_action(agent::IntFourRoomsContAgent, s; rng=Random.GLOBAL_RNG) =
    StatsBase.sample(rng, agent.μ, s[1])

function JuliaRL.start!(agent::IntFourRoomsContAgent, env_s_tp1; rng=Random.GLOBAL_RNG, kwargs...)
    agent.counter = agent.warm_up
    agent.state_t = env_s_tp1
    agent.action_t = JuliaRL.get_action(agent, env_s_tp1; rng=rng)
    return agent.action_t
end

function JuliaRL.step!(agent::IntFourRoomsContAgent, env_s_tp1, r, terminal; rng=Random.GLOBAL_RNG, kwargs...)

    # After taking action at timestep t
    μ_prob = get(agent.μ, agent.state_t[1], agent.action_t, nothing, nothing, nothing)
    sample_prob = get(agent.sample_policy, agent.state_t[1], agent.action_t, nothing, nothing, nothing)
    c, γ, π_prob = get(agent.gvf, agent.state_t, agent.action_t, env_s_tp1, nothing, nothing)
    action_tp1 = JuliaRL.get_action(agent, agent.state_t; rng=rng)

    experience = (create_features(agent.tilecoder, agent.state_t[1]./11.0),
                  create_features(agent.tilecoder, env_s_tp1[1]./11.0),
                  agent.action_t, action_tp1, c, γ, terminal, μ_prob, π_prob, (π_prob/sample_prob))

    add!(agent.WER, experience, sample_prob/μ_prob)

    if agent.counter == 0
        agent.counter = agent.training_gap
        train_value_functions(agent; rng=rng)
    end
    agent.counter -= 1

    agent.state_t = env_s_tp1
    agent.action_t = action_tp1 #JuliaRL.get_action(agent, agent.state_t; rng=rng)
    return agent.action_t
end

function train_value_functions(agent::IntFourRoomsContAgent; rng=Random.GLOBAL_RNG)

    α_arr = agent.α_arr

    # samp_er = sample(agent.ER, agent.batch_size; rng=rng)
    samp_wer = sample(agent.WER, agent.batch_size; rng=rng)

    arg_wer = (samp_wer[:ρ], samp_wer[:s_t], samp_wer[:s_tp1],
               samp_wer[:r], samp_wer[:γ_tp1], samp_wer[:terminal],
               samp_wer[:a_t], samp_wer[:a_tp1], agent.gvf.policy)

    for (α_idx, α) in enumerate(α_arr)
        agent.opt.eta = α
        update!(agent.value_arr[α_idx], agent.opt, agent.lu, arg_wer...; corr_term=1.0)
    end
end

# function predict(agent::IntFourRoomsContAgent, key, α_idx, eval_states::Array{Array{Float64, 1}, 1})
#     states = [create_features(agent.tilecoder, eval_state./11.0) for eval_state in eval_states]
#     return agent.value_dict[key][α_idx].(states)
# end

function predict!(agent::IntFourRoomsContAgent, eval_states::Array{Array{Float64, 1}, 1}, predict_arr::Array{Array{Float64, 1}})
    states = [create_features(agent.tilecoder, eval_state./11.0) for eval_state in eval_states]
    if length(predict_arr) == 0
        for i in 1:length(agent.α_arr)
            push!(predict_arr, zeros(length(eval_states)))
        end
    end
    for α_idx in 1:length(agent.α_arr)
        predict_arr[α_idx] .= agent.value_arr[α_idx].(states)
    end
end


text:
# Mortgage.jl
module Mortgage

function payment(amount, rate, years)
    monthly_rate = rate / 12.0
    factor = (1 + monthly_rate) ^ (years * 12.0)
    return amount * (monthly_rate * factor / (factor - 1))
end

end # module


text:
using Documenter
using ContextManagers

makedocs(
    sitename = "ContextManagers",
    format = Documenter.HTML(),
    modules = [ContextManagers],
    doctest = false,  # tested via test/runtests.jl
    checkdocs = :exports,  # ignore complains about non-exported docstrings
    strict = lowercase(get(ENV, "CI", "false")) == "true",
)

deploydocs(
    repo = "github.com/tkf/ContextManagers.jl",
    push_preview = true,
    # See: https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
)


text:
export NoCorrection, TinyCorrection, FixedCorrection, ClampEigenValuesCorrection

abstract type AbstractCorrection end

# Correction regularization terms for matrices

"""
    correction!(strategy, matrix)
    correction!(strategy, real)

Modifies the `matrix` with a specified correction strategy. Matrix must be squared.
Also supports real values, with the same strategies.

See also: [`NoCorrection`](@ref), [`TinyCorrection`](@ref)
"""
function correction! end

"""
    NoCorrection

One of the correction strategies for `correction!` function. Does not modify matrix and returns an original one.

See also: [`correction!`](@ref), [`TinyCorrection`](@ref)
"""
struct NoCorrection <: AbstractCorrection end

correction!(::NoCorrection, value::Real)            = value
correction!(::NoCorrection, matrix::AbstractMatrix) = matrix

"""
    TinyCorrection

One of the correction strategies for `correction!` function. Adds `ReactiveMP.tiny` term to the `matrix`'s diagonal.

See also: [`correction!`](@ref), [`NoCorrection`](@ref), [`FixedCorrection`](@ref), [`ClampEigenValuesCorrection`](@ref)
"""
struct TinyCorrection <: AbstractCorrection end

correction!(::TinyCorrection, value::Real) = clamp(value, tiny, typemax(value))

function correction!(::TinyCorrection, matrix::AbstractMatrix)
    s = size(matrix)
    @assert length(s) == 2 && s[1] === s[2]
    for i in 1:s[1]
        @inbounds matrix[i, i] += tiny
    end
    return matrix
end

"""
    FixedCorrection

One of the correction strategies for `correction!` function. Adds fixed `v` term to the `matrix`'s diagonal.

# Arguments
- `v`: fixed value to add to the matrix diagonal

See also: [`correction!`](@ref), [`NoCorrection`](@ref), [`TinyCorrection`](@ref), [`ClampEigenValuesCorrection`](@ref)
"""
struct FixedCorrection{T} <: AbstractCorrection
    v :: T
end

correction!(correction::FixedCorrection, value::Real) = clamp(value, correction.v, Inf)

function correction!(correction::FixedCorrection, matrix::AbstractMatrix)
    s = size(matrix)
    @assert length(s) == 2 && s[1] === s[2]
    for i in 1:s[1]
        @inbounds matrix[i, i] += correction.v
    end
    return matrix
end

"""
    ClampEigenValuesCorrection

One of the correction strategies for `correction!` function. Clamps eigen values of matrix to be equal or greater than fixed `v` term.

# Arguments
- `v`: fixed value used to clamp eigen values

See also: [`correction!`](@ref), [`NoCorrection`](@ref), [`FixedCorrection`](@ref), [`TinyCorrection`](@ref)
"""
struct ClampEigenValuesCorrection{T} <: AbstractCorrection
    v :: T
end

correction!(correction::ClampEigenValuesCorrection, value::Real) = clamp(value, correction.v, Inf)

function correction!(correction::ClampEigenValuesCorrection, matrix::AbstractMatrix)
    s = size(matrix)
    @assert length(s) == 2 && s[1] === s[2]

    F = svd(matrix)
    clamp!(F.S, correction.v, Inf)
    R = lmul!(Diagonal(F.S), F.Vt)
    M = mul!(matrix, F.U, R)

    return M
end


text:
## ----------  Reductions

"""
	$(SIGNATURES)

Apply `reduceFct` to each variable in each `GroupStats` object for which this is defined.
"""
function reduce_stats(
        scV :: Vector{StatsCollection{T}},
        reduceFct
        ) where T

    nSim = length(scV);
	meanStats = deepcopy(scV[1]);

	grps = get_groups(scV[1]);
	for grp in grps
		# Make vector of this field across all simulations
		# E.g., gpaYpS
        mStatsV = [get_mstats(scV[j], grp)  for j = 1 : nSim];
        newField = reduce_stats(mStatsV, reduceFct);
        replace_mstats!(meanStats, newField);
    end
    return meanStats
end

"""
	$(SIGNATURES)

Apply a scalar reduction function to all numeric fields in a `Vector` of `ModelStats`.

# Example
```
reduce
```
"""
function reduce_stats(
        statV :: AbstractVector{ModelStats{T}},
        reduceFct
        ) where T

    statsOut = deepcopy(statV[1]);
    varNameV = var_names(statV[1]);
    for varName in varNameV
        # Value from each `ModelStats` object.
        varMean = reduce_one_field(statV, varName, reduceFct);
        # Cannot reduce all types of fields.
        if !isnothing(varMean)
            set_values!(statsOut, varName, varMean);
        end
    end
    return statsOut
end

"""
    $(SIGNATURES)

Apply `reduceFct` to one field in a vector of objects.
Returns `nothing` if reduction not possible. E.g., because `reduceFct` is not defined for a given field or it returns the wrong type.
"""
function reduce_one_field(
        oVecV :: AbstractVector{ModelStats{T1}},
        pn,
        reduceFct :: Function
        ) where T1

    n = length(oVecV);
    fieldV = [get_values(oVecV[iObj], pn)  for iObj = 1 : n];
    if isa(fieldV[1], AbstractArray)
        xM = reduce_array_vector(fieldV, reduceFct);
    else
        xM = reduce_scalar_vector(fieldV, reduceFct);
    end
    return xM
end

# ------------


text:
# TensorDesc

# descriptor

mutable struct TensorDesc
    ptr::cudnnTensorDescriptor_t
end

unsafe_free!(td::TensorDesc) = cudnnDestroyTensorDescriptor(td.ptr)

Base.unsafe_convert(::Type{cudnnTensorDescriptor_t}, td::TensorDesc) = 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
    td = Ref{cudnnTensorDescriptor_t}()
    cudnnCreateTensorDescriptor(td)
    cudnnSetTensorNdDescriptor(td[], cudnnDataType(T), length(sz), sz, st)
    this = TensorDesc(td[])
    finalizer(unsafe_free!, this)
    return this
end

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

# wrappers

function cudnnAddTensor(C::DenseCuArray{T,N}, A::DenseCuArray{T,N};
                        alpha=1, beta=1) where {T,N}
    cudnnAddTensor(handle(),
                   scalingParameter(T, alpha), TensorDesc(A), A,
                   scalingParameter(T, beta ), TensorDesc(C), C)
    return C
end

# OpTensorDesc

# descriptor

mutable struct OpTensorDesc
    ptr::cudnnOpTensorDescriptor_t
end

unsafe_free!(otd::OpTensorDesc) = cudnnDestroyOpTensorDescriptor(otd.ptr)

Base.unsafe_convert(::Type{cudnnOpTensorDescriptor_t}, otd::OpTensorDesc) = otd.ptr

function OpTensorDesc(op::cudnnOpTensorOp_t, T::Type;
                      opTensorNanOpt=CUDNN_NOT_PROPAGATE_NAN)
    otd = Ref{cudnnOpTensorDescriptor_t}()
    cudnnCreateOpTensorDescriptor(otd)
    cudnnSetOpTensorDescriptor(otd[], op, cudnnDataType(T), opTensorNanOpt)
    this = OpTensorDesc(otd[])
    finalizer(unsafe_free!, this)
    return this
end

OpTensorDesc(op::cudnnOpTensorOp_t, a::DenseCuArray) = OpTensorDesc(op, eltype(a))

# wrappers

function cudnnOpTensor(op::cudnnOpTensorOp_t,
                       A::DenseCuArray{T,N}, B::DenseCuArray{T,N}, C::DenseCuArray{T,N};
                       alpha1=true, alpha2=true, beta=false) where {T,N}
    cudnnOpTensor(handle(), OpTensorDesc(op, T),
                  scalingParameter(T, alpha1), TensorDesc(A), A,
                  scalingParameter(T, alpha2), TensorDesc(B), B,
                  scalingParameter(T, beta  ), TensorDesc(C), C)
    return C
end

# ReduceTensorDesc

# descriptor

mutable struct ReduceTensorDesc
    ptr::cudnnReduceTensorDescriptor_t
end

unsafe_free!(rtd::ReduceTensorDesc) = cudnnDestroyReduceTensorDescriptor(rtd.ptr)

Base.unsafe_convert(::Type{cudnnReduceTensorDescriptor_t}, rtd::ReduceTensorDesc) = rtd.ptr

function ReduceTensorDesc(op::cudnnReduceTensorOp_t, T::Type;
                          reduceTensorNanOpt=CUDNN_NOT_PROPAGATE_NAN,
                          reduceTensorIndices=CUDNN_REDUCE_TENSOR_NO_INDICES,
                          reduceTensorIndicesType=CUDNN_32BIT_INDICES)
    rtd = Ref{cudnnReduceTensorDescriptor_t}()
    cudnnCreateReduceTensorDescriptor(rtd)
    cudnnSetReduceTensorDescriptor(rtd[], op, cudnnDataType(T), reduceTensorNanOpt,
                                   reduceTensorIndices, reduceTensorIndicesType)
    this = ReduceTensorDesc(rtd[])
    finalizer(unsafe_free!, this)
    return this
end

ReduceTensorDesc(op::cudnnReduceTensorOp_t, a::DenseCuArray) = ReduceTensorDesc(op, eltype(a))

# wrappers

function cudnnGetReductionIndicesSize(op::cudnnReduceTensorOp_t,
                                      A::DenseCuArray{T,N}, C::DenseCuArray{T,N}) where {T,N}
    size=@argout(
        cudnnGetReductionIndicesSize(
            handle(), ReduceTensorDesc(op, A),
            TensorDesc(A), TensorDesc(C),
            out(Ref{Csize_t}()))
        )[]
    return size
end

function cudnnReduceTensor(op::cudnnReduceTensorOp_t,
                           A::DenseCuArray{T,N}, C::DenseCuArray{T,N};
                           alpha=true, beta=false) where {T,N}
    # indices = Array{UInt64, 1}(undef, N)
    indicesSizeInBytes = cudnnGetReductionIndicesSize(op, A, C)
    @workspace size=@argout(
        cudnnGetReductionWorkspaceSize(
            handle(), ReduceTensorDesc(op, A),
            TensorDesc(A), TensorDesc(C),
            out(Ref{Csize_t}()))
        )[] workspace->begin
            cudnnReduceTensor(handle(), ReduceTensorDesc(op, A),
                              C_NULL, indicesSizeInBytes,
                              workspace, sizeof(workspace),
                              scalingParameter(T, alpha), TensorDesc(A), A,
                              scalingParameter(T, beta ), TensorDesc(C), C)
        end
    return C
end


text:
# Abbreviation to help keep code short!
const AMSA = AbstractMultiScaleArray

Base.map!(f::F, m::AMSA, A0::AbstractArray, As::AbstractArray...) where {F} =
        broadcast!(f, m, A0, As...)
Base.map!(f::F, m::AMSA, A0, As...) where {F} =
        broadcast!(f, m, A0, As...)

Base.BroadcastStyle(::Type{<:AMSA}) = Broadcast.ArrayStyle{AMSA}()
Base.BroadcastStyle(::Type{<:AbstractMultiScaleArrayLeaf}) = Broadcast.ArrayStyle{AbstractMultiScaleArrayLeaf}()

@inline function Base.copy(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AMSA}})
    first_amsa = find_amsa(bc)
    out = similar(first_amsa)
    copyto!(out,bc)
    out
end

@inline function Base.copy(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AbstractMultiScaleArrayLeaf}})
    first_amsa = find_amsa(bc)
    out = similar(first_amsa)
    copyto!(out,bc)
    out
end

@inline function Base.copyto!(dest::AMSA, bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AMSA}})
    if !any_non_amsa(bc)
        N = length(dest.nodes)
        for i in 1:N
            copyto!(dest.nodes[i], unpack(bc, i))
        end
        copyto!(dest.values,unpack(bc, nothing))
    else
        copyto!(dest,convert(Base.Broadcast.Broadcasted{Broadcast.DefaultArrayStyle{length(axes(bc))}}, bc))
    end
    dest
end

@inline function Base.copyto!(dest::AbstractMultiScaleArrayLeaf, bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{AbstractMultiScaleArrayLeaf}})
    if !any_non_amsa(bc)
        copyto!(dest.values,unpack(bc,nothing))
    else
        copyto!(dest,convert(Base.Broadcast.Broadcasted{Broadcast.DefaultArrayStyle{length(axes(bc))}}, bc))
    end
    dest
end

# drop axes because it is easier to recompute
@inline unpack(bc::Broadcast.Broadcasted, i) = Broadcast.Broadcasted(bc.f, unpack_args(i, bc.args))
unpack(x,::Any) = x
unpack(x::AMSA, i) = x.nodes[i]
unpack(x::AMSA, ::Nothing) = x.values

@inline unpack_args(i, args::Tuple) = (unpack(args[1], i), unpack_args(i, Base.tail(args))...)
unpack_args(i, args::Tuple{Any}) = (unpack(args[1], i),)
unpack_args(::Any, args::Tuple{}) = ()

nnodes(A) = 0
nnodes(A::AMSA) = length(A.nodes)
nnodes(bc::Broadcast.Broadcasted) = _nnodes(bc.args)
nnodes(A, Bs...) = common_number(nnodes(A), _nnodes(Bs))

@inline _nnodes(args::Tuple) = common_number(nnodes(args[1]), _nnodes(Base.tail(args)))
_nnodes(args::Tuple{Any}) = nnodes(args[1])
_nnodes(args::Tuple{}) = 0

"`A = find_amsa(As)` returns the first AMSA among the arguments."
find_amsa(bc::Base.Broadcast.Broadcasted) = find_amsa(bc.args)
find_amsa(args::Tuple) = !isempty(args) && find_amsa(find_amsa(args[1]), Base.tail(args))
find_amsa(x) = x
find_amsa(a::AMSA, rest) = a
find_amsa(::Any, rest) = find_amsa(rest)

any_non_amsa(bc::Base.Broadcast.Broadcasted) = any_non_amsa(bc.args)
any_non_amsa(args::Tuple) = any_non_amsa(any_non_amsa(args[1]), Base.tail(args))
any_non_amsa(x::AMSA) = false
any_non_amsa(x::Number) = false
any_non_amsa(x::Any) = true
any_non_amsa(x::AbstractArray) = true
any_non_amsa(x::Bool, rest) = isempty(rest) ? x : x || any_non_amsa(rest)

## utils
common_number(a, b) =
    a == 0 ? b :
    (b == 0 ? a :
     (a == b ? a :
      throw(DimensionMismatch("number of nodes must be equal"))))


text:
using MathProgComplex, DataStructures

function main()

    # problem = buildPOP_WB2(v2max=1.022, setnetworkphase=false)
    # relax_ctx = set_relaxation(problem; hierarchykind=:Real,
    #                                     # symmetries=[PhaseInvariance],
    #                                     d = 1)

    # problem = buildPOP_WB2(setnetworkphase=false)
    problem = buildPOP_WB5()
    relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                        # symmetries=[PhaseInvariance],
                                        d = 2,
                                        params = Dict(:opt_outlev=>1,
                                                      :opt_outmode=>0,
                                                      :opt_outcsv=>0))

    # println("\n--------------------------------------------------------")
    # println("problem = \n$problem")

    # println("\n--------------------------------------------------------")
    # println("relax_ctx = \n$relax_ctx")

    ########################################
    # Construction du sparsity pattern, extension chordale, cliques maximales.
    max_cliques = get_maxcliques(relax_ctx, problem)
    relax_ctx.relaxparams[:opt_nb_cliques] = length(max_cliques)

    # println("\n--------------------------------------------------------")
    # println("max cliques =")
    # println(max_cliques)

    ########################################
    # Compute moment and localizing matrices parameters: order et variables
    momentmat_param, localizingmat_param = build_sparsity(relax_ctx, problem, max_cliques)

    # println("\n--------------------------------------------------------")
    # println("moment params =")
    # for (cliquename, dcl) in momentmat_param
    #     println("Moment matrix, $cliquename \t -> dcl = $dcl")
    # end
    # for (key, (val1, val2)) in localizingmat_param
    #     @printf("%15s \t -> di-ki = %i, \tcliques = ", key, val2)
    #     for clique in val1 print("$clique, ") end
    #     @printf("\b\b \n")
    # end

    ########################################
    # Build the moment relaxation problem
    mmtrel_pb = build_momentrelaxation(relax_ctx, problem, momentmat_param, localizingmat_param, max_cliques)

    # println("\n--------------------------------------------------------")
    # println("mmtrel_pb = $mmtrel_pb")

    ########################################
    # Convert to a primal SDP problem
    sdpinstance = build_SOSrelaxation(relax_ctx, mmtrel_pb)

    # println("\n--------------------------------------------------------")
    # println("sdpinstance = \n$sdpinstance")

    path = joinpath(pwd(), "Mosek_runs", "worksdp")
    mkpath(path)
    export_SDP(sdpinstance, path)
    sdp_instance = read_SDPPrimal(path)

    # println("VAR_TYPES size:     $(size(sdp_instance.VAR_TYPES))")
    # println("BLOCKS size:        $(size(sdp_instance.BLOCKS))")
    # println("LINEAR size:        $(size(sdp_instance.LINEAR))")
    # println("CONST size:         $(size(sdp_instance.CONST))")

    sdp = SDP_Problem()

    set_constraints!(sdp, sdp_instance)
    set_vartypes!(sdp, sdp_instance)
    set_blocks!(sdp, sdp_instance)
    set_linvars!(sdp, sdp_instance)

    set_matrices!(sdp, sdp_instance)
    set_linear!(sdp, sdp_instance)
    set_const!(sdp, sdp_instance)

    # println(sdp)

    primal = SortedDict{Tuple{String,String,String}, Float64}()
    dual = SortedDict{Tuple{String, String, String}, Float64}()

    primobj, dualobj = solve_mosek(sdp::SDP_Problem, primal, dual, sol_info=relax_ctx.relaxparams)

    final_output(relax_ctx)
    # # println("Primal solution")
    # # for ((blockname, var1, var2), val) in primal
    # # @printf("%15s %5s %5s %f\n", blockname, var1, var2, val)
    # # end

    # # println("\nDual solution NEGATED")
    # # for var in problem.variables
    # #     ctrname = get_momentcstrname()
    # #     var1 = var[1]
    # #     var2 = "1"
    # #     val = dual[(ctrname, var1, var2)]
    # #     println("($(ctrname), $(var1), $(var2)) = $(-val)")
    # # end

    # println("Objectives : $primobj, $dualobj")
end

main()


text:
@testset "524.longest-word-in-dictionary-through-deleting.jl" begin
    @test find_longest_word("abpcplea", ["ale","apple","monkey","plea"]) == "apple"
    @test find_longest_word("abpcplea", ["a","b","c"]) == "a"
end