text:
@testset "PWF to Dict" begin
    @testset "Intermediary functions" begin
        file = open(joinpath(@__DIR__,"data/pwf/test_system.pwf"))

        sections = PWF._split_sections(file)
        @test isa(sections, Vector{Vector{String}})
        @test length(sections) == 5
        @test sections[1][1] == "TITU"

        data = Dict{String, Any}()
        PWF._parse_section!(data, sections[1])
        @test haskey(data, "TITU")
        PWF._parse_section!(data, sections[2])
        @test haskey(data, "DOPC IMPR")
        PWF._parse_section!(data, sections[3])
        @test haskey(data, "DCTE")
        PWF._parse_section!(data, sections[4])
        @test haskey(data, "DBAR")
        PWF._parse_section!(data, sections[5])
        @test haskey(data, "DLIN")
    end

    @testset "Resulting Dict" begin
        file = open(joinpath(@__DIR__,"data//pwf/test_system.pwf"))
        dict = PWF._parse_pwf_data(file)

        @testset "Keys" begin
            @test haskey(dict, "TITU")
            @test haskey(dict, "DOPC IMPR")
            @test haskey(dict, "DCTE")
            @test haskey(dict, "DBAR")
            @test haskey(dict, "DLIN")
        end

        @testset "Types" begin
            @test isa(dict, Dict)
            @test isa(dict["TITU"], String)
            @test isa(dict["DOPC IMPR"], Dict)
            @test isa(dict["DCTE"], Dict)
            @test isa(dict["DBAR"], Dict)
            @test isa(dict["DLIN"], Dict)
        end

        @testset "Lengths" begin
            @test length(dict["DOPC IMPR"]) == 13
            @test length(dict["DCTE"]) == 67
            @test length(dict["DBAR"]) == 9
            @test length(dict["DLIN"]) == 7

            for (idx,item) in dict["DBAR"]
                @test length(item) == 30
            end
            for (idx,item) in dict["DLIN"]
                @test length(item) == 31
            end
        end

        @testset "DBAR" begin
            for (idx,item) in dict["DBAR"]
                @test isa(item["NUMBER"], Int)
                @test isa(item["OPERATION"], Char)
                @test isa(item["STATUS"], Char)
                @test isa(item["TYPE"], Int)
                @test isa(item["BASE VOLTAGE GROUP"], String)
                @test isa(item["NAME"], String)
                @test isa(item["VOLTAGE LIMIT GROUP"], String)
                @test isa(item["VOLTAGE"], Float64)
                @test isa(item["ANGLE"], Float64)
                @test isa(item["ACTIVE GENERATION"], Float64)
                @test isa(item["REACTIVE GENERATION"], Float64)
                @test isa(item["MINIMUM REACTIVE GENERATION"], Float64)
                @test isa(item["MAXIMUM REACTIVE GENERATION"], Float64)
                @test isa(item["CONTROLLED BUS"], Int)
                @test isa(item["ACTIVE CHARGE"], Float64)
                @test isa(item["REACTIVE CHARGE"], Float64)
                @test isa(item["TOTAL REACTIVE POWER"], Float64)
                @test isa(item["AREA"], Int)
                @test isa(item["CHARGE DEFINITION VOLTAGE"], Float64)
                @test isa(item["VISUALIZATION"], Int)
                @test isa(item["REACTIVE CHARGE"], Float64)
            end
        end

        @testset "DLIN" begin
            for (idx,item) in dict["DLIN"]
                @test isa(item["FROM BUS"], Int)
                @test isa(item["OPENING FROM BUS"], Char)
                @test isa(item["OPERATION"], Char)
                @test isa(item["OPENING TO BUS"], Char)
                @test isa(item["TO BUS"], Int)
                @test isa(item["CIRCUIT"], Int)
                @test isa(item["STATUS"], Char)
                @test isa(item["OWNER"], Char)
                @test isa(item["RESISTANCE"], Float64)
                @test isa(item["REACTANCE"], Float64)
                @test isa(item["SHUNT SUSCEPTANCE"], Float64)
                @test isa(item["TAP"], Nothing) || isa(item["TAP"], Float64)
                @test isa(item["MINIMUM TAP"], Nothing)
                @test isa(item["MAXIMUM TAP"], Nothing)
                @test isa(item["LAG"], Float64)
                @test isa(item["CONTROLLED BUS"], Int)
                @test isa(item["NORMAL CAPACITY"], Float64)
                @test isa(item["EMERGENCY CAPACITY"], Float64)
                @test isa(item["NUMBER OF TAPS"], Int)
                @test isa(item["EQUIPAMENT CAPACITY"], Float64)
                @test isa(item["TRANSFORMER"], Bool)
            end
        end

        @testset "DCTE" begin
            for (key, value) in dict["DCTE"]
                @test isa(key, String)
                @test isa(value, Float64)
                @test length(key) == 4
            end
        end

        @testset "DOPC" begin
            for (key, value) in dict["DOPC IMPR"]
                @test isa(key, String)
                @test isa(value, Char)
                @test length(key) == 4
            end
        end

        @testset "TITU" begin
            @test occursin("Caso do Anderson - P gina 38", dict["TITU"])
        end
    end

    @testset "Default values" begin
        pwf = Dict{String,Any}("TITU" => "Default values test", "name" => "test_defaults",
                               "DCTE" => Dict("TEPA" => 0.1, "TEPR" => 0.1, "TLPR" => 0.1, "TLVC" => .5,
                                            "TLTC" => 0.01, "TETP" => 5.0, "TBPA" => 5.0, "TSFR" => 0.01, "TUDC" => 0.001,
                                            "TADC" => 0.01, "BASE" => 100.0, "DASE" => 100.0, "ZMAX" => 500.0, "ACIT" => 30,
                                            "LPIT" => 50, "LFLP" => 10, "LFIT" => 10, "DCIT" => 10, "VSIT" => 10, "LCRT" => 23,
                                            "LPRT" => 60, "LFCV" => 1, "TPST" => 0.2, "QLST" => 0.2, "EXST" => 0.4, "TLPP" => 1.0,
                                            "TSBZ" => 0.01, "TSBA" => 5.0, "PGER" => 30.0, "VDVN" => 40.0, "VDVM" => 200.0,
                                            "ASTP" => 0.05, "VSTP" => 5.0, "CSTP" => 5.0, "VFLD" => 70, "HIST" => 0, "ZMIN" => 0.001,
                                            "PDIT" => 10, "ICIT" => 50, "FDIV" => 2.0, "DMAX" => 5, "ICMN" => 0.05, "VART" => 5.0,
                                            "TSTP" => 33, "TSDC" => 0.02, "ASDC" => 1, "ICMV" => 0.5, "APAS" => 90, "CPAR" => 70,
                                            "VAVT" => 2.0, "VAVF" => 5.0, "VMVF" => 15.0, "VPVT" => 2.0, "VPVF" => 5.0,
                                            "VPMF" => 10.0, "VSVF" => 20.0, "VINF" => 1.0, "VSUP" => 1.0, "TLSI" => 0.0),
                                "DBAR" => Dict("1" => Dict("NUMBER" => 1, "OPERATION" => 'A', "STATUS" => 'L', "TYPE" => 0, "BASE VOLTAGE GROUP" => " 0", "NAME" => "       Bus 1","VOLTAGE LIMIT GROUP" => " 0", "VOLTAGE" => 1.0, "ANGLE" => 0.0, "ACTIVE GENERATION" => 0.0, "REACTIVE GENERATION" => 0.0, "MINIMUM REACTIVE GENERATION" => 0.0, "MAXIMUM REACTIVE GENERATION" => 0.0, "CONTROLLED BUS" => 1, "ACTIVE CHARGE" => 0.0, "REACTIVE CHARGE" => 0.0, "TOTAL REACTIVE POWER" => 0.0, "AREA" => 1, "CHARGE DEFINITION VOLTAGE" => 1.0, "VISUALIZATION" => 0, "AGGREGATOR 1" => nothing, "AGGREGATOR 2" => nothing, "AGGREGATOR 3" => nothing, "AGGREGATOR 4" => nothing, "AGGREGATOR 5" => nothing, "AGGREGATOR 6" => nothing, "AGGREGATOR 7" => nothing, "AGGREGATOR 8" => nothing, "AGGREGATOR 9" => nothing, "AGGREGATOR 10" => nothing),
                                               "2" => Dict("NUMBER" => 2, "OPERATION" => 'A', "STATUS" => 'L', "TYPE" => 2, "BASE VOLTAGE GROUP" => " 0", "NAME" => "       Bus 2","VOLTAGE LIMIT GROUP" => " 0", "VOLTAGE" => 1.0, "ANGLE" => 0.0, "ACTIVE GENERATION" => 0.0, "REACTIVE GENERATION" => 0.0, "MINIMUM REACTIVE GENERATION" => -9999.0, "MAXIMUM REACTIVE GENERATION" => 99999.0, "CONTROLLED BUS" => 2, "ACTIVE CHARGE" => 0.0, "REACTIVE CHARGE" => 0.0, "TOTAL REACTIVE POWER" => 0.0, "AREA" => 1, "CHARGE DEFINITION VOLTAGE" => 1.0, "VISUALIZATION" => 0, "AGGREGATOR 1" => nothing, "AGGREGATOR 2" => nothing, "AGGREGATOR 3" => nothing, "AGGREGATOR 4" => nothing, "AGGREGATOR 5" => nothing, "AGGREGATOR 6" => nothing, "AGGREGATOR 7" => nothing, "AGGREGATOR 8" => nothing, "AGGREGATOR 9" => nothing, "AGGREGATOR 10" => nothing)),
                                "DLIN" => Dict("1" => Dict("FROM BUS" => 1, "OPENING FROM BUS" => 'L', "OPERATION" => 'A', "OPENING TO BUS" => 'L', "TO BUS" => 2, "CIRCUIT" => 0, "STATUS" => 'L', "OWNER" => 'F', "RESISTANCE" => 0.0, "REACTANCE" => nothing, "SHUNT SUSCEPTANCE" => 0.0, "TAP" => 1.0, "MINIMUM TAP" => 0.95, "MAXIMUM TAP" => 1.05, "LAG" => 0.0, "CONTROLLED BUS" => 1, "NORMAL CAPACITY" => Inf, "EMERGENCY CAPACITY" => Inf, "NUMBER OF TAPS" => 33, "EQUIPAMENT CAPACITY" => Inf, "AGGREGATOR 1" => nothing, "AGGREGATOR 2" => nothing, "AGGREGATOR 3" => nothing, "AGGREGATOR 4" => nothing, "AGGREGATOR 5" => nothing, "AGGREGATOR 6" => nothing, "AGGREGATOR 7" => nothing, "AGGREGATOR 8" => nothing, "AGGREGATOR 9" => nothing, "AGGREGATOR 10" => nothing, "TRANSFORMER" => true),
                                               "2" => Dict("FROM BUS" => 2, "OPENING FROM BUS" => 'L', "OPERATION" => 'A', "OPENING TO BUS" => 'L', "TO BUS" => 1, "CIRCUIT" => 0, "STATUS" => 'L', "OWNER" => 'F', "RESISTANCE" => 0.0, "REACTANCE" => nothing, "SHUNT SUSCEPTANCE" => 0.0, "TAP" => 1.0, "MINIMUM TAP" => nothing, "MAXIMUM TAP" => nothing, "LAG" => 0.0, "CONTROLLED BUS" => 2, "NORMAL CAPACITY" => 0.0, "EMERGENCY CAPACITY" => 0.0, "NUMBER OF TAPS" => 33, "EQUIPAMENT CAPACITY" => 0.0, "AGGREGATOR 1" => nothing, "AGGREGATOR 2" => nothing, "AGGREGATOR 3" => nothing, "AGGREGATOR 4" => nothing, "AGGREGATOR 5" => nothing, "AGGREGATOR 6" => nothing, "AGGREGATOR 7" => nothing, "AGGREGATOR 8" => nothing, "AGGREGATOR 9" => nothing, "AGGREGATOR 10" => nothing, "TRANSFORMER" => false)),
                                "DGER" => Dict("1" => Dict("NUMBER" => 1, "OPERATION" => 'A', "MINIMUM ACTIVE GENERATION" => 0.0, "MAXIMUM ACTIVE GENERATION" => 99999.0, "PARTICIPATION FACTOR" => 0.0, "REMOTE CONTROL PARTICIPATION FACTOR" => 100., "NOMINAL POWER FACTOR" => nothing, "ARMATURE SERVICE FACTOR" => nothing,"ROTOR SERVICE FACTOR" => nothing, "CHARGE ANGLE" => nothing, "MACHINE REACTANCE" => nothing, "NOMINAL APPARENT POWER" => nothing)),
                                "DGBT" => Dict("1" => Dict("GROUP" => " 0", "VOLTAGE" => 1.0)),
                                "DGLT" => Dict("1" => Dict("GROUP" => " 0", "LOWER BOUND" => 0.8, "UPPER BOUND" => 1.2, "LOWER EMERGENCY BOUND" => 0.8, "UPPER EMERGENCY BOUND" => 1.2), "2" => Dict("GROUP" => " 1", "LOWER BOUND" => 0.9, "UPPER BOUND" => 1.1, "LOWER EMERGENCY BOUND" => 0.9, "UPPER EMERGENCY BOUND" => 1.1)),
                                "DCBA" => Dict("1" => Dict("NUMBER" => 10, "OPERATION" => 'A', "TYPE" => 1, "POLARITY" => '+', "NAME" => "RET         ", "VOLTAGE LIMIT GROUP" => " 0", "VOLTAGE" => 1.0, "GROUND ELECTRODE" => 0.0, "DC LINK" => 1),
                                               "2" => Dict("NUMBER" => 20, "OPERATION" => 'A', "TYPE" => 0, "POLARITY" => '-', "NAME" => "INV         ", "VOLTAGE LIMIT GROUP" => " 0", "VOLTAGE" => 1.0, "GROUND ELECTRODE" => 0.0, "DC LINK" => 1),
                                               "3" => Dict("NUMBER" => 30, "OPERATION" => 'A', "TYPE" => 0, "POLARITY" => '0', "NAME" => "NEUR        ", "VOLTAGE LIMIT GROUP" => " 0", "VOLTAGE" => 0.0, "GROUND ELECTRODE" => 0.0, "DC LINK" => 1),
                                               "4" => Dict("NUMBER" => 40, "OPERATION" => 'A', "TYPE" => 0, "POLARITY" => '0', "NAME" => "NEUI        ", "VOLTAGE LIMIT GROUP" => " 0", "VOLTAGE" => 0.0, "GROUND ELECTRODE" => 0.0, "DC LINK" => 1)),
                                "DCLI" => Dict("1" => Dict("FROM BUS" => 10, "OPERATION" => 'A', "TO BUS" => 20, "CIRCUIT" => 0, "OWNER" => nothing, "RESISTANCE" => nothing, "INDUCTANCE" => 0.0, "CAPACITY" => Inf)),
                                "DELO" => Dict("1" => Dict("NUMBER" => 1, "OPERATION" => 'A', "VOLTAGE" => nothing, "BASE" => 100.0, "NAME" => nothing, "HI MVAR MODE" => 'N', "STATUS" => 'L')),
                                "DBSH" => Dict("1" => Dict("FROM BUS" => 1, "OPERATION" => 'A', "TO BUS" => nothing, "CIRCUIT" => 1, "CONTROL MODE" => 'C', "MINIMUM VOLTAGE" => 0.8, "MAXIMUM VOLTAGE" => 1.2, "CONTROLLED BUS" => 1, "INITIAL REACTIVE INJECTION" => 0.0, "CONTROL TYPE" => 'C', "ERASE DBAR" => 'N', "EXTREMITY" => nothing, "REACTANCE GROUPS" => Dict("1" => Dict("GROUP" => 10, "OPERATION" => 'A', "STATUS" => 'L', "UNITIES" => 1, "OPERATING UNITIES" => 1, "REACTANCE" => nothing), "2" => Dict("GROUP" => 20, "OPERATION" => 'A', "STATUS" => 'L', "UNITIES" => 2, "OPERATING UNITIES" => 2, "REACTANCE" => nothing)))),
                                "DCSC" => Dict("1" => Dict("FROM BUS" => 1, "OPERATION" => 'A', "TO BUS" => 2, "CIRCUIT" => 0, "STATUS" => 'L', "OWNER" => 'F', "BYPASS" => 'D', "MINIMUM VALUE" => -9999.0, "MAXIMUM VALUE" => 0.0, "INITIAL VALUE" => 0.0, "CONTROL MODE" => 'X', "SPECIFIED VALUE" => nothing, "MEASUREMENT EXTREMITY" => 1, "NUMBER OF STAGES" => nothing, "NORMAL CAPACITY" => Inf, "EMERGENCY CAPACITY" => Inf, "EQUIPAMENT CAPACITY" => Inf, "AGGREGATOR 1" => nothing, "AGGREGATOR 2" => nothing, "AGGREGATOR 3" => nothing, "AGGREGATOR 4" => nothing, "AGGREGATOR 5" => nothing, "AGGREGATOR 6" => nothing, "AGGREGATOR 7" => nothing, "AGGREGATOR 8" => nothing, "AGGREGATOR 9" => nothing, "AGGREGATOR 10" => nothing)),
                                "DCAI" => Dict("1" => Dict("BUS" => 1, "OPERATION" => 'A', "GROUP" => nothing, "STATUS" => 'L', "UNITIES" => 2, "OPERATING UNITIES" => 2, "ACTIVE CHARGE" => 0.0, "REACTIVE CHARGE" => 0.0, "PARAMETER A" => nothing, "PARAMETER B" => nothing, "PARAMETER C" => nothing, "PARAMETER D" => nothing, "VOLTAGE" => 70.0, "CHARGE DEFINITION VOLTAGE" => 1.0)),
                                "DGEI" => Dict("1" => Dict("BUS" => 1, "OPERATION" => 'A', "AUTOMATIC MODE" => 'N', "GROUP" => nothing, "STATUS" => 'L', "UNITIES" => 2, "OPERATING UNITIES" => 2, "MINIMUM OPERATING UNITIES" => 1, "ACTIVE GENERATION" => 0.0, "REACTIVE GENERATION" => 0.0, "MINIMUM REACTIVE GENERATION" => -9999.0, "MAXIMUM REACTIVE GENERATION" => 99999.0, "ELEVATOR TRANSFORMER REACTANCE" => nothing, "XD" => 0.0, "XQ" => 0.0, "XL" => 0.0, "POWER FACTOR" => 1.0, "APPARENT POWER" => 99999.0, "MECHANICAL LIMIT" => 99999.0)))
        parsed_pwf = PWF.parse_file(joinpath(@__DIR__,"data/pwf/test_defaults.pwf"), pm = false)

        @test parsed_pwf == pwf
    end
end


text:
module MagMechFEM_Matlab2Julia

include(joinpath("IGA","nurbs_toolbox","NURBStoolbox.jl"));

greet() = print("Hello World!")

end # module

ext: jl


text:
using LinearAlgebra, ForwardDiff, SparseArrays
using Convex, SCS

# problem setup
n = 2

_c = 1.0 * ones(n)
_b = 1.0

"Convex.jl"
x = Variable(n)
prob = minimize(_c' * x)
prob.constraints += norm(x) <= _b

optimizer = SCS.Optimizer()
@time solve!(prob, optimizer)

optimizer.cone
optimizer.sol

# cone program data
A = Matrix(SparseMatrixCSC{Float64, Int64}(optimizer.data.A.m,
    optimizer.data.A.n,
    optimizer.data.A.colptr,
    optimizer.data.A.rowval,
    optimizer.data.A.nzval))

b = optimizer.data.b
c = optimizer.data.c

m = optimizer.data.A.m
n = optimizer.data.A.n

# cone program solution
x = optimizer.data.primal
y = optimizer.data.dual
s = optimizer.data.slack

optimizer.data.dimension

optimizer.data.num_rows

prob.model.model.optimizer.data.primal

optimizer.data.primal
optimizer.data.dual

c
s' * y
# # cones
# function κ_free(z)
#     z
# end
#
# function κ_no(z)
#     return max.(0.0, z)
# end
#
# function κ_soc(z)
#     z1 = z[1:end-1]
#     z2 = z[end]
#
#     z_proj = zero(z)
#
#     if norm(z1) <= z2
#         z_proj = copy(z)
#     elseif norm(z1) <= -z2
#         z_proj = zero(z)
#     else
#         a = 0.5 * (1.0 + z2 / norm(z1))
#         z_proj[1:end-1] = a * z1
#         z_proj[end] = a * norm(z1)
#     end
#     return z_proj
# end
#
# function Pκ(z)
#     z_proj = zero(z)
#
#     z[1] = κ_free(z[1])
#     z[2] = κ_no(z[2])
#     z[3:5] = κ_soc(z[3:5])
#
#     z_proj
# end
#
# Pκ(s)
# s

b - A * x

# Q = Array([zeros(n, n) transpose(A) c;
#            -A zeros(m, m) b;
#            -c' -b'      0.0])


text:
struct ParamsDump <: Output
  path::String
end
ParamsDump(path::AbstractString, ::Integer) = ParamsDump(path)

function beforetrajectories(d::ParamsDump, params::Simulation.SimState)
  @info "saving full parameters to $(d.path)"
  file = jldopen(d.path, "w", compress=true)
  try
    Simulation.saveparams(file, params)
  finally
    close(file)
  end
end


text:
module TestAsNumber

using Base.Test
using ArbFloats

af0 = zero(ArbFloat)
af1 = one(ArbFloat)
af2 = af1 + af1
afhalf = af1 / 2
afNaN = ArbFloat(NaN)
afInf = ArbFloat(Inf)

@test af1 <  af2
@test af1 <= af2
@test af2 >  af1
@test af2 >= af1

@test  asin(log(acosh(cosh(exp(sin(afhalf)))))) > afhalf
@test  upperbound(asin(log(atanh(tanh(exp(sin(afhalf))))))) >= afhalf
@test  lowerbound(asin(log(atanh(tanh(exp(sin(afhalf))))))) <= afhalf

end # module


text:
s = BenchmarkGroup()
suite["shortestpaths"] = s

fname, fio = mktemp()
close(fio)
Random.seed!(17)
for G in GLIST
    for n in [10_000]
        g = G(n, 5n, seed=1)
        s["dijkstra","$g"] = @benchmarkable dijkstra_shortest_paths($g, 1)
    end
end


text:
## Based on "Multi-Threading and One-Sided Communication in Parallel LU Factorization"
## http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.138.4361&rank=7

function hpl_seq(A::Matrix, b::Vector)

    blocksize = 5

    n = size(A,1)
    A = [A b]

    B_rows = linspace(0, n, div(n,blocksize)+1)
    B_rows[end] = n
    B_cols = [B_rows, [n+1]]
    nB = length(B_rows)
    depend = zeros(Bool, nB, nB) # In parallel, depend needs to be able to hold futures

    ## Small matrix case
    if nB <= 1
        x = A[1:n, 1:n] \ A[:,n+1]
        return x
    end

    ## Add a ghost row of dependencies to boostrap the computation
    for j=1:nB; depend[1,j] = true; end

    for i=1:(nB-1)
        ## Threads for panel factorizations
        I = (B_rows[i]+1):B_rows[i+1]
        #(depend[i+1,i], panel_p) = spawn(panel_factor_seq, I, depend[i,i])
        (depend[i+1,i], panel_p) = panel_factor_seq(A, I, depend[i,i])

        ## Threads for trailing updates
        for j=(i+1):nB
            J = (B_cols[j]+1):B_cols[j+1]
            #depend[i+1,j] = spawn(trailing_update_seq, I, J, panel_p, depend[i+1,i],depend[i,j])
            depend[i+1,j] = trailing_update_seq(A, I, J, panel_p, depend[i+1,i],depend[i,j])
        end
    end

    ## Completion of the last diagonal block signals termination
    #wait(depend[nB, nB])

    ## Solve the triangular system
    x = triu(A[1:n,1:n]) \ A[:,n+1]

    return x

end ## hpl()

### Panel factorization ###

function panel_factor_seq(A, I, col_dep)
    n = size (A, 1)

    ## Enforce dependencies
    #wait(col_dep)

    ## Factorize a panel
    K = I[1]:n
    panel_p = lufact!(sub(A, K, I))[:p] # Economy mode

    ## Panel permutation
    panel_p = K[panel_p]

    return (true, panel_p)

end ## panel_factor_seq()

### Trailing update ###

function trailing_update_seq(A, I, J, panel_p, row_dep, col_dep)

    n = size (A, 1)

    ## Enforce dependencies
    #wait(row_dep, col_dep)

    ## Apply permutation from pivoting
    K = (I[end]+1):n
    A[I[1]:n, J] = A[panel_p, J]

    ## Compute blocks of U
    L = tril(A[I,I],-1) + eye(length(I))
    A[I, J] = L \ A[I, J]

    ## Trailing submatrix update
    if !isempty(K)
        A[K,J] = A[K,J] - A[K,I]*A[I,J]
    end

    return true

end ## trailing_update_seq()

# This version is written for a shared memory implementation.
# The matrix A is local to the first Worker, which allocates work to other Workers
# All updates to A are carried out by the first Worker. Thus A is not distributed

hpl_par(A::Matrix, b::Vector) = hpl_par(A, b, max(1, div(max(size(A)),4)), true)

hpl_par(A::Matrix, b::Vector, bsize::Integer) = hpl_par(A, b, bsize, true)

function hpl_par(A::Matrix, b::Vector, blocksize::Integer, run_parallel::Bool)

    n = size(A,1)
    A = [A b]

    if blocksize < 1
       throw(ArgumentError("hpl_par: invalid blocksize: $blocksize < 1"))
    end

    B_rows = linspace(0, n, div(n,blocksize)+1)
    B_rows[end] = n
    B_cols = [B_rows, [n+1]]
    nB = length(B_rows)
    depend = cell(nB, nB)

    ## Small matrix case
    if nB <= 1
        x = A[1:n, 1:n] \ A[:,n+1]
        return x
    end

    ## Add a ghost row of dependencies to boostrap the computation
    for j=1:nB; depend[1,j] = true; end
    for i=2:nB, j=1:nB; depend[i,j] = false; end

    for i=1:(nB-1)
        #println("A=$A") #####
        ## Threads for panel factorizations
        I = (B_rows[i]+1):B_rows[i+1]
        K = I[1]:n
        (A_KI, panel_p) = panel_factor_par(A[K,I], depend[i,i])

        ## Write the factorized panel back to A
        A[K,I] = A_KI

        ## Panel permutation
        panel_p = K[panel_p]
        depend[i+1,i] = true

        ## Apply permutation from pivoting
        J = (B_cols[i+1]+1):B_cols[nB+1]
        A[K, J] = A[panel_p, J]
        ## Threads for trailing updates
        #L_II = tril(A[I,I], -1) + eye(length(I))
        L_II = tril(sub(A,I,I), -1) + eye(length(I))
        K = (I[length(I)]+1):n
        A_KI = A[K,I]

        for j=(i+1):nB
            J = (B_cols[j]+1):B_cols[j+1]

            ## Do the trailing update (Compute U, and DGEMM - all flops are here)
            if run_parallel
                A_IJ = A[I,J]
                #A_KI = A[K,I]
                A_KJ = A[K,J]
                depend[i+1,j] = @spawn trailing_update_par(L_II, A_IJ, A_KI, A_KJ, depend[i+1,i], depend[i,j])
            else
                depend[i+1,j] = trailing_update_par(L_II, A[I,J], A[K,I], A[K,J], depend[i+1,i], depend[i,j])
            end
        end

        # Wait for all trailing updates to complete, and write back to A
        for j=(i+1):nB
            J = (B_cols[j]+1):B_cols[j+1]
            if run_parallel
                (A_IJ, A_KJ) = fetch(depend[i+1,j])
            else
                (A_IJ, A_KJ) = depend[i+1,j]
            end
            A[I,J] = A_IJ
            A[K,J] = A_KJ
            depend[i+1,j] = true
        end

    end

    ## Completion of the last diagonal block signals termination
    @assert depend[nB, nB]

    ## Solve the triangular system
    x = triu(A[1:n,1:n]) \ A[:,n+1]

    return x

end ## hpl()

### Panel factorization ###

function panel_factor_par(A_KI, col_dep)

    @assert col_dep

    ## Factorize a panel
    panel_p = lufact!(A_KI)[:p] # Economy mode

    return (A_KI, panel_p)

end ## panel_factor_par()

### Trailing update ###

function trailing_update_par(L_II, A_IJ, A_KI, A_KJ, row_dep, col_dep)

    @assert row_dep
    @assert col_dep

    ## Compute blocks of U
    A_IJ = L_II \ A_IJ

    ## Trailing submatrix update - All flops are here
    if !isempty(A_KJ)
        m, k = size(A_KI)
        n = size(A_IJ,2)
        blas_gemm('N','N',m,n,k,-1.0,A_KI,m,A_IJ,k,1.0,A_KJ,m)
        #A_KJ = A_KJ - A_KI*A_IJ
    end

    return (A_IJ, A_KJ)

end ## trailing_update_par()

### using DArrays ###

function hpl_par2(A::Matrix, b::Vector)
    n = size(A,1)
    A = [A b]

    C = distribute(A, 2)
    nB = length(C.pmap)

    ## case if only one processor
    if nB <= 1
        x = A[1:n, 1:n] \ A[:,n+1]
        return x
    end

    depend = Array(RemoteRef, nB, nB)

    #pmap[i] is where block i's stuff is
    #block i is dist[i] to dist[i+1]-1
    for i = 1:nB
        #println("C=$(convert(Array, C))") #####
        ##panel factorization
        panel_p = remotecall_fetch(C.pmap[i], panel_factor_par2, C, i, n)

        ## Apply permutation from pivoting
        for j = (i+1):nB
           depend[i,j] = remotecall(C.pmap[j], permute, C, i, j, panel_p, n, false)
        end
        ## Special case for last column
        if i == nB
           depend[nB,nB] = remotecall(C.pmap[nB], permute, C, i, nB+1, panel_p, n, true)
        end

        ##Trailing updates
        (i == nB) ? (I = (C.dist[i]):n) :
                    (I = (C.dist[i]):(C.dist[i+1]-1))
        C_II = C[I,I]
        L_II = tril(C_II, -1) + eye(length(I))
        K = (I[length(I)]+1):n
        if length(K) > 0
            C_KI = C[K,I]
        else
            C_KI = zeros(0)
        end

        for j=(i+1):nB
            dep = depend[i,j]
            depend[j,i] = remotecall(C.pmap[j], trailing_update_par2, C, L_II, C_KI, i, j, n, false, dep)
        end

        ## Special case for last column
        if i == nB
            dep = depend[nB,nB]
            remotecall_fetch(C.pmap[nB], trailing_update_par2, C, L_II, C_KI, i, nB+1, n, true, dep)
        else
            #enforce dependencies for nonspecial case
            for j=(i+1):nB
                wait(depend[j,i])
            end
        end
    end

    A = convert(Array, C)
    x = triu(A[1:n,1:n]) \ A[:,n+1]
end ## hpl_par2()

function panel_factor_par2(C, i, n)
    (C.dist[i+1] == n+2) ? (I = (C.dist[i]):n) :
                           (I = (C.dist[i]):(C.dist[i+1]-1))
    K = I[1]:n
    C_KI = C[K,I]
    #(C_KI, panel_p) = lu!(C_KI) #economy mode
    panel_p = lu!(C_KI)[2]
    C[K,I] = C_KI

    return panel_p
end ##panel_factor_par2()

function permute(C, i, j, panel_p, n, flag)
    if flag
        K = (C.dist[i]):n
        J = (n+1):(n+1)
        C_KJ = C[K,J]

        C_KJ = C_KJ[panel_p,:]
        C[K,J] = C_KJ
    else
        K = (C.dist[i]):n
        J = (C.dist[j]):(C.dist[j+1]-1)
        C_KJ = C[K,J]

        C_KJ = C_KJ[panel_p,:]
        C[K,J] = C_KJ
    end
end ##permute()

function trailing_update_par2(C, L_II, C_KI, i, j, n, flag, dep)
    if isa(dep, RemoteRef); wait(dep); end
    if flag
        #(C.dist[i+1] == n+2) ? (I = (C.dist[i]):n) :
        #                       (I = (C.dist[i]):(C.dist[i+1]-1))
        I = C.dist[i]:n
        J = (n+1):(n+1)
        K = (I[length(I)]+1):n
        C_IJ = C[I,J]
        if length(K) > 0
            C_KJ = C[K,J]
        else
            C_KJ = zeros(0)
        end
        ## Compute blocks of U
        C_IJ = L_II \ C_IJ
        C[I,J] = C_IJ
    else
        #(C.dist[i+1] == n+2) ? (I = (C.dist[i]):n) :
        #                       (I = (C.dist[i]):(C.dist[i+1]-1))

        I = (C.dist[i]):(C.dist[i+1]-1)
        J = (C.dist[j]):(C.dist[j+1]-1)
        K = (I[length(I)]+1):n
        C_IJ = C[I,J]
        if length(K) > 0
            C_KJ = C[K,J]
        else
            C_KJ = zeros(0)
        end

        ## Compute blocks of U
        C_IJ = L_II \ C_IJ
        C[I,J] = C_IJ
        ## Trailing submatrix update - All flops are here
        if !isempty(C_KJ)
            cm, ck = size(C_KI)
            cn = size(C_IJ,2)
            blas_gemm('N','N',cm,cn,ck,-1.0,C_KI,cm,C_IJ,ck,1.0,C_KJ,cm)
            #C_KJ = C_KJ - C_KI*C_IJ
            C[K,J] = C_KJ
        end
    end
end ## trailing_update_par2()

## Test n*n matrix on np processors
## Prints 5 numbers that should be close to zero
function test(n, np)
    A = rand(n,n); b = rand(n);
    X = (@elapsed x = A \ b);
    Y = (@elapsed y = hpl_par(A,b, max(1,div(n,np))));
    Z = (@elapsed z = hpl_par2(A,b));
    for i=1:(min(5,n))
        print(z[i]-y[i], " ")
    end
    println()
    return (X,Y,Z)
end

## test k times and collect average
function test(n,np,k)
    sum1 = 0; sum2 = 0; sum3 = 0;
    for i = 1:k
        (X,Y,Z) = test(n,np)
        sum1 += X
        sum2 += Y
        sum3 += Z
    end
    return (sum1/k, sum2/k, sum3/k)
end

max_stars_repo_name: alanedelman/julia

max_issues_repo_name: Aerlinger/julia

max_forks_repo_name: Aerlinger/julia


text:
using PyPlot, PyCall, Revise
using Kwant
kwant = Kwant; pyplot = PyPlot

include("superconductor.jl")

syst = make_system(W=10)

# Check that the system looks as intended.
kwant.plot(syst)

# Finalize the system.
system = syst.finalized()

# Check particle-hole symmetry of the scattering matrix
check_PHS(system)

# Compute and plot the conductance
plot_conductance(system, energies=[0.002 * i for i in range(-10, stop=100)])


text:
"""
    InitStateBC

Set the value at the boundary to match the `init_state_prognostic!` function. This is
mainly useful for cases where the problem has an explicit solution.

# TODO: This should be fixed later once BCs are figured out (likely want
# different things here?)
"""
struct InitStateBC end
function atmos_boundary_state!(
    ::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    bc::InitStateBC,
    m::AtmosModel,
    state⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    aux⁻::Vars,
    bctype,
    t,
    _...,
)
    init_state_prognostic!(m, state⁺, aux⁺, aux⁺.coord, t)
end

function atmos_normal_boundary_flux_second_order!(
    nf,
    bc::InitStateBC,
    atmos,
    fluxᵀn,
    n⁻,
    state⁻,
    diff⁻,
    hyperdiff⁻,
    aux⁻,
    state⁺,
    diff⁺,
    hyperdiff⁺,
    aux⁺,
    bctype,
    t,
    args...,
)

    normal_boundary_flux_second_order!(
        nf,
        atmos,
        fluxᵀn,
        n⁻,
        state⁻,
        diff⁻,
        hyperdiff⁻,
        aux⁻,
        state⁺,
        diff⁺,
        hyperdiff⁺,
        aux⁺,
        bc,
        t,
        args...,
    )

end

function boundary_state!(
    ::NumericalFluxSecondOrder,
    m::AtmosModel,
    state⁺::Vars,
    diff⁺::Vars,
    aux⁺::Vars,
    n⁻,
    state⁻::Vars,
    diff⁻::Vars,
    aux⁻::Vars,
    bc::InitStateBC,
    t,
    args...,
)
    init_state_prognostic!(m, state⁺, aux⁺, aux⁺.coord, t)
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
include("solver_preamble.jl")
import CDCS
# Iterations:
# dsos_concave_then_convex_cubic : > 2000, < 3000
# chebyshev : > 12000, < 12500
factory = optimizer_with_attributes(CDCS.Optimizer, "verbose" => 0, "maxIter" => 12500)
# chebyshev : > 2e-3, < 3e-3
config = MOI.Test.TestConfig(atol=3e-3, rtol=3e-3, query=false)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config, [
        "sdsos_options_pricing"
    ])
end
@testset "SDP" begin
    Tests.sd_test(factory, config, [
        "sos_options_pricing"
    ])
end


text:
using CUDA
using GemmKernels
using GPUArrays

M = parse(Int, ARGS[1])
N = parse(Int, ARGS[2])
K = parse(Int, ARGS[3])

transpose_a = (ARGS[4] == "n" ? false : ARGS[4] == "t" ? true : error("Invalid memory layout for A: $(ARGS[4])"))
transpose_b = (ARGS[5] == "n" ? false : ARGS[5] == "t" ? true : error("Invalid memory layout for B: $(ARGS[5])"))

# Large matrices take way too long
if N > 4096
    exit()
end

function benchmark_matmul(a, b, c, d)
    CUDA.@sync begin
        GPUArrays.generic_matmatmul!(c, b, a, 1, 1)
    end
end

a_h = rand(Complex{Float16}, (M, K)) / sqrt(Float16(K))
b_h = rand(Complex{Float16}, (K, N)) / sqrt(Float16(K))
c_h = rand(Complex{Float32}, (M, N))

# Transpose input if necessary
a_h = transpose_a ? transpose(a_h) : a_h
b_h = transpose_b ? transpose(b_h) : b_h

a   = CuArray(a_h)
b   = CuArray(b_h)
c   = CuArray(c_h)
d   = similar(c)

# warmup
benchmark_matmul(a, b, c, d)

# profile
for i = 1 : 10
    CUDA.@profile benchmark_matmul(a, b, c, d)
end


text:
# We illustrate here the use of `ListofStates` in dealing with a warm start procedure.

# `ListofStates` is designed to store the of the iteration process.
# In this tutorial, we compare the resolution of a convex unconstrained problem with 3 variants:
#  - a steepest descent method
#  - an inverse-BFGS method
#  - a mix of 5 steps of steepest descent and then switching to a BFGS initialized with the 5 previous steps.

using Stopping, NLPModels, LinearAlgebra, Test, Printf

# First, we introduce our two implementations that both uses an backtracking Armijo linesearch.

import Stopping.armijo
function armijo(xk, dk, fk, slope, f)
  t = 1.0
  fk_new = f(xk + dk)
  while f(xk + t * dk) > fk + 1.0e-4 * t * slope
    t /= 1.5
    fk_new = f(xk + t * dk)
  end
  return t, fk_new
end

function steepest_descent(stp :: NLPStopping)

  xk = stp.current_state.x
  fk, gk = objgrad(stp.pb, xk)

  OK = update_and_start!(stp, fx = fk, gx = gk)

  @printf "%2s %9s %7s %7s %7s\n" "k" "fk" "||∇f(x)||" "t" "λ"
  @printf "%2d %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score)
  while !OK
    dk = - gk
    slope = dot(dk, gk)
    t, fk = armijo(xk, dk, fk, slope, x->obj(stp.pb, x))
    xk += t * dk
    fk, gk = objgrad(stp.pb, xk)

    OK = update_and_stop!(stp, x = xk, fx = fk, gx = gk)

    @printf "%2d %9.2e %7.1e %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score) t slope
  end
  return stp
end

function bfgs_quasi_newton_armijo(stp :: NLPStopping; Hk = nothing)

  xk = stp.current_state.x
  fk, gk = objgrad(stp.pb, xk)
  gm = gk

  dk, t = similar(gk), 1.
  if isnothing(Hk)
    Hk = I #start from identity matrix
  end

  OK = update_and_start!(stp, fx = fk, gx = gk)

  @printf "%2s %7s %7s %7s %7s\n" "k" "fk" "||∇f(x)||" "t" "cos"
  @printf "%2d %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score)

  while !OK
    if stp.meta.nb_of_stop != 0
      sk = t * dk
      yk = gk - gm
      ρk = 1/dot(yk, sk)
      #we need yk'*sk > 0 for instance yk'*sk ≥ 1.0e-2 * sk' * Hk * sk
      Hk = ρk ≤ 0.0 ? Hk : (I - ρk * sk * yk') * Hk * (I - ρk * yk * sk') + ρk * sk * sk'
      if norm(sk) ≤ 1e-14
        break
      end
      #H2 = Hk + sk * sk' * (dot(sk,yk) + yk'*Hk*yk )*ρk^2 - ρk*(Hk * yk * sk' + sk * yk'*Hk)
    end
    dk = - Hk * gk
    slope = dot(dk, gk) # ≤ -1.0e-4 * norm(dk) * gnorm
    t, fk = armijo(xk, dk, fk, slope, x->obj(stp.pb, x))

    xk = xk + t * dk
    gm = copy(gk)
    gk = grad(stp.pb, xk)

    OK = update_and_stop!(stp, x = xk, fx = fk, gx = gk)
    @printf "%2d %7.1e %7.1e %7.1e %7.1e\n" stp.meta.nb_of_stop fk norm(stp.current_state.current_score) t slope
  end
  stp.stopping_user_struct = Dict( :Hk => Hk)
  return stp
end

# We consider the following convex unconstrained problem model using `ADNLPModels.jl` and defines a related `NLPStopping`.

fH(x) = (x[2]+x[1].^2-11).^2+(x[1]+x[2].^2-7).^2
nlp = ADNLPModel(fH, [10., 20.])
stp = NLPStopping(nlp, optimality_check = unconstrained_check,
                 atol = 1e-6, rtol = 0.0, max_iter = 100)

reinit!(stp, rstate = true, x = nlp.meta.x0)
steepest_descent(stp)

@test status(stp) == :Optimal
@test stp.listofstates == VoidListofStates()

elapsed_time(stp)
nlp.counters

reinit!(stp, rstate = true, x = nlp.meta.x0, rcounters = true)
bfgs_quasi_newton_armijo(stp)

@test status(stp) == :Optimal
@test stp.listofstates == VoidListofStates()

elapsed_time(stp)
nlp.counters

NLPModels.reset!(nlp)
stp_warm = NLPStopping(nlp, optimality_check = unconstrained_check,
                      atol = 1e-6, rtol = 0.0, max_iter = 5,
                      n_listofstates = 5) #shortcut for list = ListofStates(5, Val{NLPAtX{Float64,Array{Float64,1},Array{Float64,2}}}()))
steepest_descent(stp_warm)
@test status(stp_warm) == :IterationLimit
@test length(stp_warm.listofstates) == 5

Hwarm = I
for i=2:5
  sk = stp_warm.listofstates.list[i][1].x - stp_warm.listofstates.list[i-1][1].x
  yk = stp_warm.listofstates.list[i][1].gx - stp_warm.listofstates.list[i-1][1].gx
  ρk = 1/dot(yk, sk)
  if ρk > 0.0
    global Hwarm = (I - ρk * sk * yk') * Hwarm * (I - ρk * yk * sk') + ρk * sk * sk'
  end
end

reinit!(stp_warm)
stp_warm.meta.max_iter = 100
bfgs_quasi_newton_armijo(stp_warm, Hk = Hwarm)
status(stp_warm)

elapsed_time(stp_warm)
nlp.counters


text:
"""
    xfit(object::Union{Pca, Plsr}, X; nlv = nothing)
Matrix fitting from a PCA or PLS model
* `object` : The fitted model.
* `X` : X-data to be approximatred from the model.
* `nlv` : Nb. components (PCs or LVs) to consider. If nothing,
    it is the maximum nb. of components.

Compute an approximate of matrix `X` (X_fit) from a PCA
or PLS fitted on `X`.
"""
function xfit(object::Union{Pca, Plsr}, X; nlv = nothing)
    xfit!(object, copy(X); nlv = nlv)
end

function xfit!(object::Union{Pca, Plsr}, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    if(nlv == 0)
        m = size(X, 1)
        @inbounds for i = 1:m
            X[i, :] .= object.xmeans
        end
    else
        P = @view(object.P[:, 1:nlv])
        mul!(X, transform(object, X; nlv = nlv), P')
        center!(X, -object.xmeans)
    end
    X
end

"""
    xresid(object::Union{Pca, Plsr}, X; nlv = nothing)
Residual matrix after fitting by a PCA or PLS model
* `object` : The fitted model.
* `X` : X-data for which the residuals have to be computed.
* `nlv` : Nb. components (PCs or LVs) to consider. If nothing,
    it is the maximum nb. of components.

Compute the residual matrix E = X - X_fit.
"""
function xresid(object::Union{Pca, Plsr}, X; nlv = nothing)
    xresid!(object, copy(X); nlv = nlv)
end

function xresid!(object::Union{Pca, Plsr}, X; nlv = nothing)
    a = size(object.T, 2)
    isnothing(nlv) ? nlv = a : nlv = min(nlv, a)
    if(nlv == 0)
        center!(X, object.xmeans)
    else
        X .= X .- xfit(object, X; nlv = nlv)
    end
    X
end


text:
### THIS IS WIP, runner under construction

### Serial tempering Monte Carlo runner

export SerialTempMC

immutable SerialTemperingMCRunner <: SerialMCRunner
  burnin::Int
  thinning::Int
  nsteps::Int
  r::Range{Int}
  swapperiod::Int

  function SerialTemperingMCRunner(r::Range{Int}, p::Int)
    burnin = first(r)-1
    thinning = r.step
    nsteps = last(r)
    @assert burnin >= 0 "Number of burn-in iterations should be non-negative."
    @assert thinning >= 1 "Thinning should be >= 1."
    @assert nsteps > burnin "Total number of MCMC iterations should be greater than number of burn-in iterations."
    @assert p > 0 "Swap period should be positive."
    new(burnin, thinning, nsteps, r, p)
  end
end

SerialTemperingMCRunner(r::UnitRange{Int}, p::Int) = SerialTemperingMCRunner(first(r):1:last(r), p)

SerialTemperingMCRunner(; burnin::Int=0, thinning::Int=1, nsteps::Int=100, swapperiod::Int=5) =
  SerialTemperingMCRunner((burnin+1):thinning:nsteps, swapperiod)

typealias SerialTempMC SerialTemperingMCRunner

function run(systems::Vector{MCSystem})
  tic()

  local nchains = length(systems)
  local modelsize = systems[1].model.size
  local runner = systems[1].runner
  local jobtype = typeof(systems[1].job)

  @assert all(map(s->s.model.size, systems) .== modelsize) "Models must have the same parameter vector size."
  @assert isa(runner, SerialTempMC) "Runners must be of SerialTempMC type."
  @assert all(map(s->s.runner, systems) .== runner) "The runners in the vector of MC systems must be identical."
  @assert all(map(s->s.job, systems) .== jobtype) "Jobs must be of the same type."

	local logweights = zeros(nchains)  # log of task weights that will be adapted
  local at1 = 1  # Pick starting task

  # Pre-allocation for storing results
  mcchain::MCChain = MCChain(modelsize, length(runner.r))

  local mcstate1 = systems[at1].job.receive()
  local ppars, logtarget, pars
  ppars, logtarget, pars = mcstate1.ppars, mcstate1.logtarget, mcstate1.pars

  # Sampling loop
  i::Int = 1
  for j in 1:runner.nsteps
    if j % swapperiod == 0  # Attempt a task switch
      at2 = rand(1:(nchains-1)) # Pick another task at1 random
      at2 = at2 >= at1 ? at2+1 : at2

      mcstate2 = systems[at2].job.reset(pars)
      if rand() < exp(logtarget - mcstate2.logtarget + logweights[at2] - logweights[at1]) # accept swap ?
        at1, mcstate1 = at2, mcstate2
      end
    else
      mcstate1 = consume(systems[at1].task)
    end

    # TODO : add logweights adaptation
    ppars, logtarget, pars = mcstate1.ppars, mcstate1.logtarget, mcstate1.pars

    if in(j, runner.r)
      mcchain.samples[i, :] = ppars
      mcchain.logtargets[i] = logtarget

      i += 1
      #mcchain.misc[:mod][pos] = at1
    end
  end

  mcchain.runtime = toq()
  mcchain
end

# TODO: check that all elements of array contain MCMCTasks of the same type
function resume(systems::Array{MCMCTask}; nsteps::Int=100)
  run(MCMCTask[systems[j].model * systems[j].sampler * SerialTempMC(nsteps=nsteps, swapperiod=systems[j].runner.swapperiod)
    for j in 1:length(systems)])
end


text:
pf = API(joinpath(@__DIR__, "resources", "defs.jl"), sym -> @eval(vk, $(Meta.parse("$sym"))))

# println(pf.structs)
# println(pf.funcs)
# println(pf.consts)
# println(pf.enums)

s1 = SDefinition("MyMutableStruct", true)
s2 = SDefinition("MyStruct", false, fields=("a" => "Int", "b" => "Ptr{Cvoid}", "c" => "NTuple{16,Cfloat}"))

f1 = FDefinition("f", Signature("f", PositionalArgument.(["x", "y", "z"], ["Int", nothing, nothing]), [KeywordArgument("p", 5)]), false, [Statement("println(50)"), Statement("println(\"finished\")")])
f2 = FDefinition("g", Signature("g", [PositionalArgument("x", nothing)], []), true, [Statement("50 * x")])

c1 = CDefinition("myconst", 50)

e1 = EDefinition(name="my_enum", fields=["a", "b", "c", "d", "e", "f"], with_begin_block=false)
e2 = EDefinition(name="my_other_enum", fields=["a = 3", "b = 2", "c=200"], with_begin_block=true)

@testset "Parsing" begin
    @testset "Source parsing" begin
        @test generate(pf.structs["MyMutableStruct"]) == generate(s1)
        @test generate(pf.structs["MyStruct"]) == generate(s2)
        @test generate(pf.funcs["f"]) == generate(f1)
        @test generate(pf.funcs["g"]) == generate(f2)
        @test generate(pf.consts["myconst"]) == generate(c1)
        @test generate(pf.enums["my_enum"]) == generate(e1)
        @test generate(pf.enums["my_other_enum"]) == generate(e2)
    end
    # @testset "Declaration generation interoperability" begin

    # end
end

"""
mutable struct MyMutableStruct end
struct MyStruct
    a::Int
    b::Ptr{Cvoid}
    c::NTuple{16,Cfloat}
end

function f(x::Int, y, z; p=5)
    println(50)
    println("finished")
end

g(x) = 50 * x

@enum my_enum a b c d e f

const myconst = 50
"""


text:
function processInitPopAges( ages::Vector{Float64}, ageRes::Float64 )

    result = DataFrame()

    # Summary statistics.
    minAge = minimum( ages )
    maxAge = maximum( ages )
    ageSummary = [mean( ages ), median( ages ),
        length( ages ) == 1 ? missing : std( ages ), minAge, maxAge]

    # Age binning.
    minAge = floor( minAge / ageRes ) * ageRes
    maxAge = floor( maxAge / ageRes ) * ageRes
    ageGrid = collect( minAge:ageRes:maxAge )
    counts = map( age -> count( ages .>= age ), ageGrid )
    counts -= vcat( counts[2:end], 0 )

    # Collating.
    result[:, :age] = vcat( ageGrid, ["mean", "median", "stdev", "min", "max"] )
    result[:, :counts] = vcat( counts, ageSummary )
    return result

end  # processInitPopAges( ages::Vector{Float64}, ageRes::Float64 )


text:
rg(M) = M.*G./C./C

function dudθ(θ, tgα, R)
    cθ = cos(θ)
    sθ = sin(θ)
    (cθ+sθ.*tgα)./R./(sθ-cθ.*tgα)
end

function dduddθ(θ, tgα, R)
    cθ = cos(θ)
    sθ = sin(θ)

#        V                 W
#    ------------      -----------
#    (cθ+sθ.*tgα)./R./(sθ-cθ.*tgα)

    v = cθ+sθ.*tgα
    w = sθ-cθ.*tgα

    dv = -sθ+cθ.*tgα
    dw =  cθ-sθ.*tgα

#   1/R*{  V*W^(-1) }'
    1/R*( dv/w-v/w/w*dw)
end

function ray_angle(θ, r, dudθ)
    cθ = cos(θ)
    sθ = sin(θ)
    -atan((-r.*dudθ.*sθ+cθ)./(-r.*dudθ.*cθ-sθ))
end

function mylinspace(i,f,N)
    l = linspace(0,1,div(N,2)).^0.001
    ll = π/2-i
    lr = f-π/2
    [(1-reverse(l)).*ll+i; l.*lr+π/2]
end

function extrapolate_to_meeting_same_side(x1,y1,x2,y2)
    s1 = Spline1D(x1,y1, k=1, bc="error")
    s2 = Spline1D(x2,y2, k=1, bc="error")
    X = fzero(x->evaluate(s1,x)-evaluate(s2,x), 0.0, min(x1[end],x2[end]))
    Y = evaluate(s1,X)

    xx1 = x1[x1.<X]
    yy1 = y1[x1.<X]
    xx2 = x2[x2.<X]
    yy2 = y2[x2.<X]
    push!(xx1, X)
    push!(xx2, X)

    push!(yy1, Y)
    push!(yy2, Y)

    xx1,yy1, xx2,yy2
end

function extrapolate_to_x2!(x1,y1,x2)
    ixs = find(x1.<x2)
    x1 = x1[ixs]
    y1 = y1[ixs]

    dx1 = x1[end]-x1[end-10]
    dy1 = y1[end]-y1[end-10]
    tg1 = dy1/dx1

    DX1 = x2-x1[end]
    DY1 = DX1*tg1

    X = x2
    Y = y1[end] + DY1

    push!(x1, X)
    push!(y1, Y)

    x1, y1
end

function to_small(v)
    Float64[float(vi) for vi in v]
end

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
export AltLadder, altladder

type AltLadder <: Model
  S2 :: Int
  L :: Int
  lz :: Tuple{Float64, Float64}
  lx :: Tuple{Complex128, Complex128}
  rz :: Float64
  rx :: Complex128
  twistedmask :: Tuple{Bool, Bool, Bool}
  coeff :: BondCoeffs
  function AltLadder(S2::Integer, L::Integer, Jlz::Real, Jlx::Real, delta::Real, Jrz::Real, Jrx::Real, twistedbond::Symbol)
    c = bond_coeffs(S2)
    if twistedbond == :rung
      twistedmask = (true, false, false)
    elseif twistedbond == :plus_leg
      twistedmask = (false, true, false)
    elseif twistedbond == :minus_leg
      twistedmask = (false, false, true)
    else
      error("twistedbond should be :rung, :plus_leg, or :minus_leg")
    end
    lz = (Jlz*(1.0+delta), Jlz*(1.0-delta))
    lx = (Jlx*complex(1.0+delta), Jlx*complex(1.0-delta))
    rz = Jrz
    rx = complex(Jrx)
    new(S2, L, lz, lx, rz, rx, twistedmask, c)
  end
end
altladder(S2::Integer, L::Integer, Jlz::Real, Jlx::Real, delta::Real, Jrz::Real, Jrx::Real, twist::Symbol) = AltLadder(S2, L, Jlz, Jlx, delta, Jrz, Jrx, twist)
num_sites(ladder::AltLadder) = 2ladder.L
num_bonds(ladder::AltLadder) = 3ladder.L

function interall(ladder::AltLadder, t::Real=0.0)
  const io = open("interall.def", "w")
  header_interall(io, ladder)

  L = ladder.L

  twist = cis(t)
  cone = one(Complex128)
  bond(io, 0, L, ladder.rz, ladder.rx * ifelse(ladder.twistedmask[1], twist, cone), ladder.coeff)
  bond(io, 0, 1, ladder.lz[1], ladder.lx[1] * ifelse(ladder.twistedmask[2], twist, cone), ladder.coeff)
  bond(io, L, L+1, ladder.lz[2], ladder.lx[2] * ifelse(ladder.twistedmask[3], twist, cone), ladder.coeff)

  @inbounds for i in 1:(L-1)
    j = (i+1)%L

    parity = i%2
    # first leg
    bond(io, i, j, ladder.lz[parity+1], ladder.lx[parity+1], ladder.coeff)
    # second leg
    bond(io, i+L, j+L, ladder.lz[2-parity], ladder.lx[2-parity], ladder.coeff)
  end

  ## rung
  @inbounds for i in 1:(L-1)
    bond(io, i, i+L, ladder.rz, ladder.rx, ladder.coeff)
  end

  close(io)
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
"""
    MC_NewtonGS(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                         YdH_mc::Array{SMCg{N,T},2},YH_mc::Array{SMCg{N,T},2},
                         nx::Int64,np::Int64,optc::Vector{Any})

Performs a Newton Gauss-Seidel bounding iteration using McCormick relaxations
where `x_mc` is the initial state variable, `z_mc` are relaxations of the affine
function in `x_mc`, `YdH_mc` is the preconditioned jacobian, `YH_mc` is the
equality constraints jacobian, `nx` size of the state variables, `np` is the
size of the decision space, and `optc` is a vector with stored values.
"""
function MC_NewtonGS!(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                     YdH_mc::Union{Vector{SMCg{N,T}},Array{SMCg{N,T},2}},
                     YH_mc::Vector{SMCg{N,T}},nx::Int64) where {N,T<:AbstractFloat}

  S1::SMCg{N,T} = zero(SMCg{N,T})
  x_mc_int::Vector{SMCg{N,T}} = copy(x_mc)
  for i=1:nx
    S1 = zero(SMCg{N,T})
    for j=1:nx
      if (i<j)
        S1 = S1 + YdH_mc[i,j]*(x_mc[j]-z_mc[j])
      elseif (j>i)
        S1 = S1 + YdH_mc[i,j]*(x_mc[j]-z_mc[j])
      end
    end
    x_mc[i] = z_mc[i] - (YH_mc[i]+S1)/YdH_mc[i,i]
    x_mc[i] = Final_Cut(x_mc[i],x_mc_int[i])
  end
end

"""
    MC_KrawczykCW(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                           YdH_mc::Array{SMCg{N,T},2},YH_mc::Array{SMCg{N,T},2},
                           nx::Int64,np::Int64,optc::Vector{Any})

Performs a componentwise Krawczyk bounding iteration using McCormick relaxations
where `x_mc` is the initial state variable, `z_mc` are relaxations of the affine
function in `x_mc`, `YdH_mc` is the preconditioned jacobian, `YH_mc` is the
equality constraints jacobian, `nx` size of the state variables, `np` is the
size of the decision space, and `optc` is a vector with stored values.
"""
function MC_KrawczykCW!(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                        YdH_mc::Union{Vector{SMCg{N,T}},Array{SMCg{N,T},2}},
                        YH_mc::Vector{SMCg{N,T}},nx::Int64) where{N,T<:AbstractFloat}
  #println("z_mc: $z_mc")
  #println("x_mc: $x_mc")
  #println("YdH_mc: $YdH_mc")
  #println("YH_mc: $YH_mc")
  #println("nx: $nx")
  S1::SMCg{N,T} = zero(SMCg{N,T})
  x_mc_int::Vector{SMCg{N,T}} = copy(x_mc)
  #println("YH_mc: $S1")
  #println("x_mc_int: $x_mc_int")
  for i=1:nx
    #println("i: $i")
    S1 = zero(SMCg{N,T})
    for j=1:nx
      #println("j: $j")
      if (i<j)
        S1 = S1 - (YdH_mc[i,j])*(x_mc[j]-z_mc[j])
      elseif (j>i)
        S1 = S1 - (YdH_mc[i,j])*(x_mc[j]-z_mc[j])
      else
        S1 = S1 + (one(T)-YdH_mc[i,j])*(x_mc[j]-z_mc[j])
      end
    end
    x_mc[i] =  z_mc[i] - YH_mc[i] + S1
    #println("x_mc[$i]: $(x_mc[i])")
    x_mc[i] = Final_Cut(x_mc[i],x_mc_int[i])
    #println("x_mc[$i]: $(x_mc[i])")
  end
end

"""
    MCn_NewtonGS(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                 YdH_mc::Union{Vector{SMCg{N,T}},SparseMatrixCSC{SMCg{N,T},Int64}},
                 YH_mc::Array{SMCg{N,T},2},nx::Int64,np::Int64,optc::Vector{Any})

Performs a Newton Gauss-Seidel bounding iteration using McCormick relaxations
where `x_mc` is the initial state variable, `z_mc` are relaxations of the affine
function in `x_mc`, `YdH_mc` is the preconditioned jacobian, `YH_mc` is the
equality constraints jacobian, `nx` size of the state variables, `np` is the
size of the decision space, and `optc` is a vector with stored values.
"""
function MCn_NewtonGS!(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                     YdH_mc::Union{Vector{SMCg{N,T}},SparseMatrixCSC{SMCg{N,T},Int64}},
                     YH_mc::Vector{SMCg{N,T}},nx::Int64) where {N,T<:AbstractFloat}

  S1::SMCg{N,T} = zero(SMCg{N,T})
  S2::SMCg{N,T} = zero(SMCg{N,T})
  S3::SMCg{N,T} = zero(SMCg{N,T})
  x_mc_int::Vector{SMCg{N,T}} = copy(x_mc)
  for i=1:nx
    S1 = zero(SMCg{N,T})
    S2 = zero(SMCg{N,T})
    S3 = zero(SMCg{N,T})
    for q=(YdH_mc.colptr[i]):(YdH_mc.colptr[i+1]-1)
      if (i < YdH_mc.rowval[q])
        S1 += YdH_mc.nzval[q]*(x_mc[YdH_mc.rowval[q]]-z_mc[YdH_mc.rowval[q]])
      elseif (i > YdH_mc.rowval[q])
        S2 += YdH_mc.nzval[q]*(x_mc[YdH_mc.rowval[q]]-z_mc[YdH_mc.rowval[q]])
      else
        S3 = YdH_mc.nzval[q]
      end
    end
    x_mc[i] = z_mc[i] - (YH_mc[i]+S1+S2)/S3
    x_mc[i] = Final_Cut(x_mc[i],x_mc_int[i])
  end
end

"""
    MCn_KrawczykCW(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                           YdH_mc::Array{SMCg{N,T},2},YH_mc::Array{SMCg{N,T},2},
                           nx::Int64,np::Int64,optc::Vector{Any})

Performs a componentwise Krawczyk bounding iteration using McCormick relaxations
where `x_mc` is the initial state variable, `z_mc` are relaxations of the affine
function in `x_mc`, `YdH_mc` is the preconditioned jacobian, `YH_mc` is the
equality constraints jacobian, `nx` size of the state variables, `np` is the
size of the decision space, and `optc` is a vector with stored values.
"""
function MCn_KrawczykCW!(z_mc::Vector{SMCg{N,T}},x_mc::Vector{SMCg{N,T}},
                        YdH_mc::Union{Vector{SMCg{N,T}},SparseMatrixCSC{SMCg{N,T},Int64}},
                        YH_mc::Vector{SMCg{N,T}},nx::Int64) where{N,T<:AbstractFloat}
  S1::SMCg{N,T} = zero(SMCg{N,T})
  x_mc_int::Vector{SMCg{N,T}} = copy(x_mc)
  for i=1:nx
    S1 = zero(SMCg{N,T})
    for q=(YdH_mc.colptr[i]):(YdH_mc.colptr[i+1]-1)
      if (i == YdH_mc.rowval[q])
        S1 += (one(T)-YdH_mc.nzval[q])*(x_mc[YdH_mc.rowval[q]]-z_mc[YdH_mc.rowval[q]])
      else
        S1 += -YdH_mc.nzval[q]*(x_mc[YdH_mc.rowval[q]]-z_mc[YdH_mc.rowval[q]])
      end
    end
    x_mc[i] =  z_mc[i] - YH_mc[i] + S1
    x_mc[i] = Final_Cut(x_mc[i],x_mc_int[i])
  end
end


text:
println("Immutable:")
a = 10
b = a
b = 20
@show a

println("\nNo copy:")
a = [10]
b = a
b[1] = 20
@show a

println("\nCopy:")
a = [10]
b = copy(a)
b[1] = 20
@show a

println("\nShallow copy:")
a = [[10]]
b = copy(a)
b[1][1] = 20
@show a

println("\nDeep copy:")
a = [[10]]
b = deepcopy(a)
b[1][1] = 20
@show a;

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
#=
    Tiled LU factorization in pure Julia. The serial performance essentially
    comes from the `TriangularSolve` and `LoopVectorization`, and
    `RecursiveFactorization` packages. The parallelization is handled by
    `DataFlowTask`s.
=#

lu(A::Matrix,args...) = lu!(deepcopy(A),args...)

function lu!(A::Matrix,s=TILESIZE[],tturbo::Val{T}=Val(false)) where {T}
    _lu!(PseudoTiledMatrix(A,s),tturbo)
end

function _lu!(A::PseudoTiledMatrix,tturbo::Val{T}=Val(false)) where {T}
    m,n = size(A)
    for i in 1:m
        Aii = A[i,i]
        # TODO: for simplicity, no pivot is allowed. Pivoting the diagonal
        # blocks requires permuting the corresponding row/columns before continuining
        @dspawn RecursiveFactorization.lu!(Aii,NoPivot(),tturbo) (Aii,) (RW,)
        # @dspawn LinearAlgebra.lu!(Aii) (Aii,) (RW,)
        for j in i+1:n
            Aij = A[i,j]
            Aji = A[j,i]
            @dspawn begin
                TriangularSolve.ldiv!(UnitLowerTriangular(Aii),Aij,tturbo)
                TriangularSolve.rdiv!(Aji,UpperTriangular(Aii),tturbo)
            end (Aii,Aij,Aji) (R,RW,RW)
            # TriangularSolve.ldiv!(UnitLowerTriangular(Aii),Aij,tturbo)
            # TriangularSolve.rdiv!(Aji,UpperTriangular(Aii),tturbo)
        end
        for j in i+1:m
            for k in i+1:n
                Ajk = A[j,k]
                Aji = A[j,i]
                Aik = A[i,k]
                @dspawn schur_complement!(Ajk,Aji,Aik,tturbo) (Ajk,Aji,Aik) (RW,R,R)
                # schur_complement!(Ajk,Aji,Aik,tturbo)
            end
        end
    end
    # wait for all computations before returning
    DataFlowTasks.sync()
    return LU(A.data,BlasInt[],zero(BlasInt))
end

# a fork-join approach for comparison with the data-flow parallelism
function _lu_forkjoin!(A::PseudoTiledMatrix,tturbo::Val{T}=Val(false)) where {T}
    m,n = size(A)
    for i in 1:m
        Aii = A[i,i]
        # FIXME: for simplicity, no pivot is allowed. Pivoting the diagonal
        # blocks requires permuting the corresponding row/columns before continuining
        RecursiveFactorization.lu!(Aii,NoPivot(),tturbo)
        Threads.@threads for j in i+1:n
            Aij = A[i,j]
            Aji = A[j,i]
            TriangularSolve.ldiv!(UnitLowerTriangular(Aii),Aij,tturbo)
            TriangularSolve.rdiv!(Aji,UpperTriangular(Aii))
        end
        @sync for j in i+1:m
            for k in i+1:n
                Ajk = A[j,k]
                Aji = A[j,i]
                Aik = A[i,k]
                Threads.@spawn schur_complement!(Ajk,Aji,Aik,tturbo)
            end
        end
    end
    # wait for all computations before returning
    DataFlowTasks.sync()
    return LU(A.data,BlasInt[],zero(BlasInt))
end


text:
function install_mfem()
    PWD = pwd()
    change_directory()
    http_file("https://bit.ly/mfem-4-1", "mfem-4-1.tgz")
    uncompress("mfem-4-1.tgz", "mfem-4.1")
    str = String(read("mfem-4.1/CMakeLists.txt"))
    str = replace(str, "add_library(mfem \${SOURCES} \${HEADERS} \${MASTER_HEADERS})"=>"""add_library(mfem SHARED \${SOURCES} \${HEADERS} \${MASTER_HEADERS})
set_property(TARGET mfem PROPERTY POSITION_INDEPENDENT_CODE ON)""")
    open("mfem-4.1/CMakeLists.txt", "w") do io
        write(io, str)
    end
    change_directory("mfem-4.1/build")
    require_file("build.ninja") do
        ADCME.cmake(CMAKE_ARGS = ["SHARED=YES", "STATIC=NO"])
    end
    require_library("mfem") do
        ADCME.make()
    end
    mv("libmfem.so.4.1.0", joinpath(ADCME.LIBDIR, "libmfem.so"), force=true)
    cd(PWD)
end


text:
using Markdown
struct YAXVarInfo
  project::String
  longname::String
  units::String
  url::String
  comment::String
  reference::String
end
Base.isless(a::YAXVarInfo, b::YAXVarInfo) = isless(string(a.project, a.longname),string(b.project, b.longname))

import Base.show
function show(io::IO,::MIME"text/markdown",v::YAXVarInfo)
    un=v.units
    url=v.url
    re=v.reference
    pr = v.project
    ln = v.longname
    co = v.comment
    mdt=md"""
### $ln
*$(co)*

* **Project** $(pr)
* **units** $(un)
* **Link** $(url)
* **Reference** $(re)
"""
    mdt[3].items[1][1].content[3]=[" $pr"]
    mdt[3].items[2][1].content[3]=[" $un"]
    mdt[3].items[3][1].content[3]=[" $url"]
    mdt[3].items[4][1].content[3]=[" $re"]
    show(io,MIME"text/markdown"(),mdt)
end
show(io::IO,::MIME"text/markdown",v::Vector{YAXVarInfo})=foreach(x->show(io,MIME"text/markdown"(),x),v)

"""
    cubeinfo(cube)

Shows the metadata and citation information on variables contained in a cube.
"""
function cubeinfo(ds::YAXArray, variable="unknown")
    p = ds.properties
    vi=YAXVarInfo(
      get(p,"project_name", "unknown"),
      get(p,"long_name",variable),
      get(p,"units","unknown"),
      get(p,"url","no link"),
      get(p,"comment",variable),
      get(p,"references","no reference")
    )
end

export cubeinfo


text:
"""
Arguments for configuring loss surface computation.

# Parameters
- `xmin::Float32`, `xmax::Float32`, `xnum::Int32`: define span of the surface
and amount of point in the `x` direction.
- `ymin::Float32`, `ymax::Float32`, `ynum::Int32`: define span of the surface
and amount of point in the `y` direction.
- `use_gpu::Bool`: Whether to use gpu. If `true` then it is up to you to
transfer mini-batch in `evaluation_step` function to the gpu.
- `x_directions_file::Union{Nothing, String}`: If provided, directions
for `x` axis will be loaded from it. Otherwise, random initialized.
Should match in shape with model's weights.
- `y_directions_file::Union{Nothing, String}`: If provided, directions
for `y` axis will be loaded from it. Otherwise, random initialized.
Should match in shape with model's weights.

!!! note
    If `use_gpu = true` then it is up to you to transfer
    mini-batch in `evaluation_step` function to the gpu.

!!! note
    Both x- and y- directions files in should contain directions
    under `checkpoint_weights` key.

    Saving weights like this in code would look like
    `BSON.@save file checkpoint_weights`.
"""
@with_kw struct SurfaceArgs
    xmin::Float32 = -1
    xmax::Float32 = 1
    xnum::Int32 = 20

    ymin::Float32 = -1
    ymax::Float32 = 1
    ynum::Int32 = 20

    batch_size::Int32 = 64
    use_gpu::Bool = true

    x_directions_file::Union{Nothing, String} = nothing
    y_directions_file::Union{Nothing, String} = nothing
end

function evaluation_loop(
    model, dataloader::DataLoader, evaluation_step::Function, args::SurfaceArgs,
)::Float32
    total = Float32(length(dataloader) * dataloader.batchsize)
    total_loss = 0.0f0
    for batch in dataloader
        batch_loss = evaluation_step(model, batch, args) * dataloader.batchsize
        total_loss += batch_loss |> cpu
    end
    total_loss / total
end

function get_directions(
    model_parameters::Params, args::SurfaceArgs,
)::Tuple{Params, Params}
    directions = Array{Params}(undef, 2)

    directions_files = (args.x_directions_file, args.y_directions_file)
    is_targets = [isa(df, String) && isfile(df) for df in directions_files]

    @inbounds for (i, (is_target, target_file)) in enumerate(zip(is_targets, directions_files))
        if is_target
            @load target_file checkpoint_weights
            i_directions = create_directions(model_parameters, params(checkpoint_weights))
        else
            i_directions = create_directions(model_parameters)
        end
        directions[i] = i_directions
    end
    tuple(directions...)
end

"""
```julia
function create_surface(
    model, dataloader::DataLoader, evaluation_step::Function, args::SurfaceArgs,
)
```

Create loss surface.

# Parameters
- `model`: Model to use in loss function.
- `dataloader::DataLoader`: Dataset on which to evaluate loss function.
- `evaluation_step`:
    Custom-defined function which given model, mini-batch and args,
    computes loss on that mini-batch.
- `args::SurfaceArgs`: Parameters used when computing surface.

!!! note
    If you specified `use_gpu` in `args`, then it is up to you,
    to transfer mini-batch in `evaluation_step` function to gpu.
"""
function create_surface(
    model, dataloader::DataLoader, evaluation_step::Function, args::SurfaceArgs,
)
    model_parameters = params(model)
    directions = get_directions(model_parameters, args)

    if args.use_gpu
        model = model |> gpu
    end
    testmode!(model)

    x_coordinates = range(args.xmin, stop=args.xmax, length=args.xnum) |> collect
    y_coordinates = range(args.ymin, stop=args.ymax, length=args.ynum) |> collect

    loss_surface = fill(-1.0f0, args.xnum, args.ynum)
    it = product(enumerate(y_coordinates), enumerate(x_coordinates))
    @showprogress 2 "Creating surface " for ((j, y), (i, x)) in it
        shifted = shift_weights(model_parameters, directions, (x, y))
        loadparams!(model, shifted)

        eval_loss = evaluation_loop(model, dataloader, evaluation_step, args)
        loss_surface[i, j] = eval_loss
    end

    (x_coordinates, y_coordinates), loss_surface
end


text:
#   Debug.Runtime:
# ==================
# Scope data type used at runtime

module Runtime
using Debug.AST, Debug.Meta
import Base: haskey, isequal, getindex, setindex!, keys, ==
export Scope, ModuleScope, LocalScope, getter, setter, get_eval
export Frame, parent_frame, enclosing_scope_frame, scope_frameof

# ---- Scope: runtime symbol table with getters and setters -------------------

abstract Scope

type ModuleScope <: Scope
    eval::Function
end
type LocalScope <: Scope
    parent::Scope
    syms::Dict
    env::Env
end

haskey(s::ModuleScope, sym::Symbol) = false
haskey(s::LocalScope,  sym::Symbol) = haskey(s.syms, sym) || haskey(s.parent, sym)

get_entry(scope::ModuleScope, sym::Symbol) = error("No local variable `$sym` found in scope")
function get_entry(scope::LocalScope, sym::Symbol)
    haskey(scope.syms, sym) ? scope.syms[sym] : get_entry(scope.parent, sym)
end

getter(scope::Scope, sym::Symbol) = get_entry( scope, sym)[1]
setter(scope::Scope, sym::Symbol) = get_entry( scope, sym)[2]
getindex(scope::Scope, sym::Symbol) = getter(scope,sym)()
setindex!(scope::Scope, x,  sym::Symbol) = setter(scope, sym)(x)

get_eval(scope::ModuleScope) = scope.eval
get_eval(scope::LocalScope)  = get_eval(scope.parent)

keys(scope::Scope) = getkeys!(Set{Symbol}(), scope)
getkeys!(syms::Set{Symbol}, scope::ModuleScope) = syms
function getkeys!(syms::Set{Symbol}, scope::LocalScope)
    union!(getkeys!(syms, scope.parent), keys(scope.syms))
end

# ---- Frame: Runtime node instance -------------------------------------------

type Frame
    node::Node
    scope::Scope
end
isequal(f1::Frame, f2::Frame) = (f1.node === f2.node && f1.scope === f2.scope)
==(f1::Frame, f2::Frame)      = (f1.node === f2.node && f1.scope === f2.scope)

function parent_frame(f::Frame)
    @assert !is_function(f.node)
    node = parentof(f.node)
    Frame(node, introduces_scope(node) ? f.scope.parent : f.scope)
end

enclosing_scope_frame(f::Frame) = scope_frameof(parent_frame(f))

scope_frameof(f::Frame) = is_scope_node(f.node) ? f : parent_frame(f)

end # module


text:
# The material below are explicity exported with the package EffectiveWaves.jl, so no need to include this file.

## possible materials

# note that many of these materials below are automatically loaded, see "src/EffectiveWaves.jl".

# http://www.rfcafe.com/references/general/velocity-sound-media.htm
  Brick        = Medium(ρ=1.80*1000, c = 3650.0)
  IronArmco    = Medium(ρ=7.85*1000, c = 5960.0)
  LeadAnnealed = Medium(ρ=11.4*1000, c = 2160.)
  RubberGum    = Medium(ρ=0.95*1000, c = 1550.0)
  FusedSilica  = Medium(ρ=2.2*1000,  c = 5760.0)
  GlassPyrex   = Medium(ρ=2.32*1000, c = 5640.0)
  ClayRock     = Medium(ρ=2.2*1000,  c = 3480.0)

# Fluids # http://www.rshydro.co.uk/sound-speeds/
  WaterDistilled= Medium(ρ=0.998*1000, c = 1496.0)
  Glycerol      = Medium(ρ=1.26*1000,  c = 1904.0)
  Hexadecane    = Medium(ρ=0.773*1000, c = 1338.0)
  Acetone       = Medium(ρ=0.791*1000, c = 1174.0)
  Benzene       = Medium(ρ=0.87*1000,  c = 1295.0)
  Nitrobenzene  = Medium(ρ=1.204*1000, c = 1415.0)
  OliveOil      = Medium(ρ=0.912*1000, c = 1431.0)
  SodiumNitrate = Medium(ρ=1.8*1000,   c = 1760.)
# gases
  AirDry       = Medium(ρ=1.293,     c = 331.4)

# non-porous minerals, i.e. add air-bubbles to be more realistic
# https://www.spec2000.net/05-mineralprops.htm
  ρ=1.58*1000;
  Clay         = Medium(ρ=ρ, c = p_speed(ρ,1.5*10^9,1.4*10^9))
  ρ=2.65*1000;
  SilicaQuartz = Medium(ρ=ρ, c = p_speed(ρ,36.6*10^9,45.0*10^9))
  ρ=2.71*1000;
  Calcite = Medium(ρ=ρ, c = p_speed(ρ,70.8*10^9,32.0*10^9)) # approximately limestone
  LimeStone    = Medium(ρ=2460.0,  c = 4855.0) # 4% porosity, base level http://www.sciencedirect.com/science/article/pii/S1365160915300101


text:
module SummitGem

using ..Ahorn, Maple

# Not placeable, way to hardcoded
const placements = Ahorn.PlacementDict()

function gemSprite(entity::Maple.SummitGem)
    index = get(entity.data, "gem", 0)

    return "collectables/summitgems/$index/gem00.png"
end

function Ahorn.selection(entity::Maple.SummitGem)
    x, y = Ahorn.position(entity)
    sprite = gemSprite(entity)

    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SummitGem, room::Maple.Room)
    sprite = gemSprite(entity)

    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

# Nothing to render, way to hardcoded to work with
Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Maple.SummitGemManager, room::Maple.Room) = true

end


text:
using Oceananigans.BoundaryConditions
using Oceananigans.BoundaryConditions: BCType

import Oceananigans.BoundaryConditions: bctype_str, print_condition

struct HaloCommunication <: BCType end

const HaloCommunicationBC = BoundaryCondition{<:HaloCommunication}

bctype_str(::HaloCommunicationBC) ="HaloCommunication"

HaloCommunicationBoundaryCondition(val; kwargs...) = BoundaryCondition(HaloCommunication, val; kwargs...)

struct HaloCommunicationRanks{F, T}
    from :: F
      to :: T
end

HaloCommunicationRanks(; from, to) = HaloCommunicationRanks(from, to)

print_condition(hcr::HaloCommunicationRanks) = "(from rank $(hcr.from) to rank $(hcr.to))"

function inject_halo_communication_boundary_conditions(field_bcs, local_rank, connectivity)
    rank_east = connectivity.east
    rank_west = connectivity.west
    rank_north = connectivity.north
    rank_south = connectivity.south
    rank_top = connectivity.top
    rank_bottom = connectivity.bottom

    east_comm_ranks = HaloCommunicationRanks(from=local_rank, to=rank_east)
    west_comm_ranks = HaloCommunicationRanks(from=local_rank, to=rank_west)
    north_comm_ranks = HaloCommunicationRanks(from=local_rank, to=rank_north)
    south_comm_ranks = HaloCommunicationRanks(from=local_rank, to=rank_south)
    top_comm_ranks = HaloCommunicationRanks(from=local_rank, to=rank_top)
    bottom_comm_ranks = HaloCommunicationRanks(from=local_rank, to=rank_bottom)

    east_comm_bc = HaloCommunicationBoundaryCondition(east_comm_ranks)
    west_comm_bc = HaloCommunicationBoundaryCondition(west_comm_ranks)
    north_comm_bc = HaloCommunicationBoundaryCondition(north_comm_ranks)
    south_comm_bc = HaloCommunicationBoundaryCondition(south_comm_ranks)
    top_comm_bc = HaloCommunicationBoundaryCondition(top_comm_ranks)
    bottom_comm_bc = HaloCommunicationBoundaryCondition(bottom_comm_ranks)

    x_bcs = CoordinateBoundaryConditions(isnothing(rank_west) ? field_bcs.west : west_comm_bc,
                                         isnothing(rank_east) ? field_bcs.east : east_comm_bc)

    y_bcs = CoordinateBoundaryConditions(isnothing(rank_south) ? field_bcs.south : south_comm_bc,
                                         isnothing(rank_north) ? field_bcs.north : north_comm_bc)

    z_bcs = CoordinateBoundaryConditions(isnothing(rank_bottom) ? field_bcs.bottom : bottom_comm_bc,
                                         isnothing(rank_top) ? field_bcs.top : top_comm_bc)

    return FieldBoundaryConditions(x_bcs, y_bcs, z_bcs)
end


text:
#=
Image augmentation will apply various transformations to the training images
and accompanying labels
=#

mutable struct ImageAugmentation
    rotation::Bool
    rotation_angles::Array{Float64,1}
    flip_x::Bool
    flip_y::Bool
    set_bw::Bool
end

function ImageAugmentation()

    default_rotations = [pi/6, pi/2, pi, 3*pi/2, -pi/6]
    ImageAugmentation(true,default_rotations,true,true,true)
end

function rotation_augmentation(im,ll,im_out,l_out,ind,aug::ImageAugmentation,count)
    for j in aug.rotation_angles
        for k=1:size(im,3)
            im_out[:,:,k,count] = imrotate(im[:,:,k,ind],j,ImageTransformations.Flat())[1:size(im,1),1:size(im,2)]
        end
        for k=1:size(ll,3)
            l_out[:,:,k,count] = imrotate(ll[:,:,k,ind],j,0)[1:size(ll,1),1:size(ll,2)]
        end
        count += 1
    end
    count
end

function flip_x_augmentation(im,ll,im_out,l_out,ind,aug::ImageAugmentation,count)
    im_out[:,:,:,count] = reverse(im[:,:,:,ind],dims=1)
    l_out[:,:,:,count] = reverse(ll[:,:,:,ind],dims=1)
    count+=1
end

function flip_y_augmentation(im,ll,im_out,l_out,ind,aug::ImageAugmentation,count)
    im_out[:,:,:,count] = reverse(im[:,:,:,ind],dims=2)
    l_out[:,:,:,count] = reverse(ll[:,:,:,ind],dims=2)
    count+=1
end

function random_block(xdim,ydim,max_length)

    x1 = rand(1:div(xdim,2))
    x2 = rand(x1:(x1+max_length))
    y1 = rand(1:div(ydim,2))
    y2 = rand(y1:(y1+max_length))

    (x1,x2,y1,y2)
end

function flip_bw_augmentation(im,ll,im_out,l_out,ind,aug::ImageAugmentation,count)

    #Random Black
    (x1,x2,y1,y2) = random_block(size(im,1),size(im,2),5)
    im_out[:,:,:,count] = im[:,:,:,ind]
    l_out[:,:,:,count] = ll[:,:,:,ind]
    im_out[x1:x2,y1:y2,:,count] .= 0.0
    count+=1

    #Random White
    (x1,x2,y1,y2) = random_block(size(im,1),size(im,2),5)
    im_out[:,:,:,count] = im[:,:,:,ind]
    l_out[:,:,:,count] = ll[:,:,:,ind]
    im_out[x1:x2,y1:y2,:,count] .= 1.0
    count+=1
end

function image_augmentation(im,ll,aug::ImageAugmentation)

    rot_size = 0
    flip_x_size = 0
    flip_y_size = 0
    set_bw_size = 0
    if aug.rotation
        rot_size += length(aug.rotation_angles)
    end
    if aug.flip_x
        flip_x_size = 1
    end
    if aug.flip_y
        flip_y_size = 1
    end
    if aug.set_bw
        set_bw_size = 2
    end
    array_size = rot_size + flip_x_size + flip_y_size + set_bw_size + 1

    im_out = zeros(Float32,size(im,1),size(im,2),size(im,3),size(im,4) * array_size)
    l_out = zeros(Float32,size(ll,1),size(ll,2),size(ll,3),size(ll,4) * array_size)

    count=1

    for i=1:size(im,4)
        im_out[:,:,:,count] = im[:,:,:,i]
        l_out[:,:,:,count] = ll[:,:,:,i]
        count += 1
        #rotations
        if aug.rotation
            count = rotation_augmentation(im,ll,im_out,l_out,i,aug,count)
        end

        #Flip x
        if aug.flip_x
            count = flip_x_augmentation(im,ll,im_out,l_out,i,aug,count)
        end

        #Flip Y
        if aug.flip_y
            count = flip_y_augmentation(im,ll,im_out,l_out,i,aug,count)
        end

        if aug.set_bw
            count = flip_bw_augmentation(im,ll,im_out,l_out,i,aug,count)
        end

        # 0.75 Scale (zoom in)

        #1.25 Scale (Reflect around places where image is empty)
    end

    #Shuffle Positions
    myinds=Random.shuffle(collect(1:size(im_out,4)));
    im_out=im_out[:,:,:,myinds]
    l_out=l_out[:,:,:,myinds]

    (im_out,l_out)
end

function normalize_images(ii)

    mean_img = mean(ii,dims=4)[:,:,:,1]
    std_img = std(ii,dims=4)[:,:,:,1]
    std_img[std_img .== 0.0] .= 1

    ii = (ii .- mean_img) ./ std_img

    min_ref = minimum(ii)
    ii = ii .- min_ref

    max_ref = maximum(ii)
    ii = ii ./ max_ref

    (mean_img,std_img,min_ref,max_ref)
end

function normalize_new_images(ii::KnetArray,mean_img::Array,std_img,min_ref,max_ref)
    normalize_new_images(ii,convert(KnetArray,mean_img),convert(KnetArray,std_img),min_ref,max_ref)
end

function normalize_new_images(ii,mean_img)
    ii = ii ./ 255
    ii = ii .- mean_img
end


text:
module Test_clyap

using LinearAlgebra
using MatrixEquations
using Test

@testset "Testing continuous Lyapunov equation solvers" begin

n = 10
Ty = Float64

@testset "Continuous Lyapunov equations" begin

reltol = eps(float(100))
a = -1; b = 2im; @time x = lyapc(a,b)
@test abs(a*x+x*a'+b) < reltol

reltol = eps(float(100f0))
a = 1f0-2f0im; b = 2f0; @time x = lyapc(a,b)
@test abs(a*x+x*a'+b) < reltol

try
  x = lyapc(zeros(1,1),ones(1,1))
  @test false
catch
  @test true
end

try
  x = lyapc([0 1;  -1 0 ],ones(2,2))
  @test false
catch
  @test true
end

try
  x = lyapc([1 1;  0 -1],ones(2,2))
  @test false
catch
  @test true
end

try
  x = lyapc([1 1 1 1;  -1 1 0 1; 0 0 -1 1; 0 0 -1 -1],ones(4,4))
  @test false
catch
  @test true
end

for Ty in (Float64, Float32)

ar = rand(Ty,n,n);
ac = rand(Ty,n,n)+im*rand(Ty,n,n);
c = rand(Ty,n,n)+im*rand(Ty,n,n);
qc = c'*c;
Qr = real(qc);
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time x = lyapc(ac,qc);
@test norm(ac*x+x*ac'+qc)/norm(x)/norm(ac) < reltol

α = 3+im; # α = im  # SingularException
@time x = lyapc(α*I,qc);
@test norm(α*x+x*α'+qc)/norm(x) < reltol

α = 3+im; # α = im  # SingularException
@time x = lyapc(α,Qr);
@test norm(α*x+x*α'+Qr)/norm(x) < reltol

@time x = lyapc(ac',qc);
@test norm(ac'*x+x*ac+qc)/norm(x)/norm(ac) < reltol

@time x = lyapc(ar,Qr)
@test norm(ar*x+x*ar'+Qr)/norm(x)/norm(ar)  < reltol

@time x = lyapc(ar',Qr)
@test norm(ar'*x+x*ar+Qr)/norm(x)/norm(ar) < reltol

@time x = lyapc(ac,Qr);
@test norm(ac*x+x*ac'+Qr)/norm(x)/norm(ac) < reltol

@time x = lyapc(ac',Qr);
@test norm(ac'*x+x*ac+Qr)/norm(x)/norm(ac) < reltol

@time x = lyapc(ar,qc);
@test norm(ar*x+x*ar'+qc)/norm(x)/norm(ar) < reltol

@time x = lyapc(ar',qc);
@test norm(ar'*x+x*ar+qc)/norm(x)/norm(ar)  < reltol
end
end

@testset "Continuous generalized Lyapunov equations" begin

reltol = eps(float(100))
a = -1+im; ee = 3im; b = 2; @time x = lyapc(a,ee,b)
@test abs(a*x*ee'+ee*x*a'+b) < reltol

reltol = eps(float(100f0))
a = 1f0-2f0im; ee = 3f0im; b = 2f0; @time x = lyapc(a,ee,b)
@test abs(a*x*ee'+ee*x*a'+b) < reltol

try
  x = lyapc(zeros(1,1),ones(1,1),ones(1,1))
  @test false
catch
  @test true
end

try
  x = lyapc([0 2;  -2 0 ],[2 0; 0 2],ones(2,2))
  @test false
catch
  @test true
end

try
  x = lyapc([2 2;  0 -2],[2 0; 0 2],ones(2,2))
  @test false
catch
  @test true
end

try
  x = lyapc(2*[1 1 1 1;  -1 1 0 1; 0 0 -1 1; 0 0 -1 -1],2*Matrix{Float64}(I,4,4), ones(4,4))
  @test false
catch
  @test true
end

for Ty in (Float64, Float32)

ar = rand(Ty,n,n)
ac = rand(Ty,n,n)+im*rand(Ty,n,n)
er = rand(Ty,n,n)
ec = er+im*rand(Ty,n,n)

c = rand(Ty,n,n)+im*rand(Ty,n,n)
qc = c'*c
Qr = real(qc)
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

@time x = lyapc(ac,ec,qc);
@test norm(ac*x*ec'+ec*x*ac'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

β = 3
@time x = lyapc(ac,β*I,qc);
@test norm(ac*x*β'+β*x*ac'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

β = (1+im);
@time x = lyapc(ac,β,qc);
@test norm(ac*x*β'+β*x*ac'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

α = 2+3im
@time x = lyapc(α*I,ec,qc);
@test norm(α*x*ec'+ec*x*α'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

α = 2+3im; β = (1+im); # α = im; β = 1  # SingularException
@time x = lyapc(α,β,qc);
@test norm(α*x*β'+β*x*α'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

α = 2+3im; β = (1+im); # α = im; β = 1  # SingularException
@time x = lyapc(α*I,β*I,qc);
@test norm(α*x*β'+β*x*α'+qc)/norm(x)/norm(ac)/norm(ec) < reltol

@time x = lyapc(ac',ec',qc);
@test norm(ac'*x*ec+ec'*x*ac+qc)/norm(x)/norm(ac)/norm(ec)  < reltol

β = (1+im);
@time x = lyapc(ac',β*I,qc);
@test norm(ac'*x*β+β'*x*ac+qc)/norm(x)/norm(ac)/norm(ec)  < reltol

α = 2+3im
@time x = lyapc(α*I,ec',qc);
@test norm(α'*x*ec+ec'*x*α+qc)/norm(x)/norm(ac)/norm(ec) < reltol

@time x = lyapc(ac',ec,qc);
@test norm(ac'*x*ec'+ec*x*ac+qc)/norm(x)/norm(ac)/norm(ec)  < reltol

@time x = lyapc(ac,ec',qc);
@test norm(ac*x*ec+ec'*x*ac'+qc)/norm(x)/norm(ac)/norm(ec)  < reltol

@time x = lyapc(ar,er,Qr);
@test norm(ar*x*er'+er*x*ar'+Qr)/norm(x)/norm(ar)/norm(er)  < reltol

@time x = lyapc(ar,er,qc);
@test norm(ar*x*er'+er*x*ar'+qc)/norm(x)/norm(ar)/norm(er)  < reltol

@time x = lyapc(ar',er',Qr);
@test norm(ar'*x*er+er'*x*ar+Qr)/norm(x)/norm(ar)/norm(er) < reltol

@time x = lyapc(ar',er,Qr);
@test norm(ar'*x*er'+er*x*ar+Qr)/norm(x)/norm(ar)/norm(er) < reltol

@time x = lyapc(ar,er',Qr);
@test norm(ar*x*er+er'*x*ar'+Qr)/norm(x)/norm(ar)/norm(er) < reltol

@time x = lyapc(ac,ec,Qr);
@test norm(ac*x*ec'+ec*x*ac'+Qr)/norm(x)/norm(ac)/norm(ec)   < reltol

@time x = lyapc(ac',ec',Qr);
@test norm(ac'*x*ec+ec'*x*ac+Qr)/norm(x)/norm(ac)/norm(ec)   < reltol

@time x = lyapc(ac',ec,Qr);
@test norm(ac'*x*ec'+ec*x*ac+Qr)/norm(x)/norm(ac)/norm(ec)   < reltol

@time x = lyapc(ac,ec',Qr);
@test norm(ac*x*ec+ec'*x*ac'+Qr)/norm(x)/norm(ac)/norm(ec)   < reltol

@time x = lyapc(ar,er,qc);
@test norm(ar*x*er'+er*x*ar'+qc)/norm(x)/norm(ar)/norm(er) < reltol

@time x = lyapc(ar',er',qc);
@test norm(ar'*x*er+er'*x*ar+qc)/norm(x)/norm(ar)/norm(er) < reltol

@time x = lyapc(ar',er,qc);
@test norm(ar'*x*er'+er*x*ar+qc)/norm(x)/norm(ar)/norm(er) < reltol

@time x = lyapc(ar,er',qc);
@test norm(ar*x*er+er'*x*ar'+qc)/norm(x)/norm(ar)/norm(er) < reltol
end
end

@testset "Continuous Lyapunov equations - Schur form" begin

for Ty in (Float64, Float32)

ar = rand(Ty,n,n);
ac = rand(Ty,n,n)+im*rand(Ty,n,n);
er = rand(Ty,n,n);
ec = er+im*rand(Ty,n,n);
as, es = schur(ar,er);
acs, ecs = schur(ac,ec);

c = rand(Ty,n,n)+im*rand(Ty,n,n);
qc = c'*c;
Qr = real(qc);
Ty == Float64 ? reltol = eps(float(100)) : reltol = eps(100*n*one(Ty))

x = copy(qc)
@time lyapcs!(acs,ecs,x);
@test norm(acs*x*ecs'+ecs*x*acs'+qc)/norm(x)/norm(acs)/norm(ecs) < reltol

x = copy(qc)
@time lyapcs!(acs,I,x);
@test norm(acs*x+x*acs'+qc)/norm(x)/norm(acs)/norm(ecs) < reltol

x = copy(qc)
@time lyapcs!(acs,ecs,x,adj=true);
@test norm(acs'*x*ecs+ecs'*x*acs+qc)/norm(x)/norm(acs)/norm(ecs) < reltol

# test
x = copy(Qr);
@time lyapcs!(as,es,x);
@test norm(as*x*es'+es*x*as'+Qr)/norm(x)/norm(as)/norm(es) < reltol

x = copy(Qr);
@time lyapcs!(as,es,x,adj=true);
@test norm(as'*x*es+es'*x*as+Qr)/norm(x)/norm(as)/norm(es) < reltol

x = copy(Qr)
@time lyapcs!(as,I,x);
@test norm(as*x+x*as'+Qr)/norm(x)/norm(as) < reltol

x = copy(Qr)
@time lyapcs!(as,x);
@test norm(as*x+x*as'+Qr)/norm(x)/norm(as) < reltol

if Ty == Float64
try
  x = convert(Matrix{Complex{Float32}},copy(qc))
  @time lyapcs!(as,x);
  @test false
catch
  @test true
end
end

x = copy(Qr)
@time lyapcs!(as,es,x,adj=true);
@test norm(as'*x*es+es'*x*as+Qr)/norm(x)/norm(as)/norm(es) < reltol

x = copy(qc)
@time lyapcs!(acs,x);
@test norm(acs*x+x*acs'+qc)/norm(x)/norm(acs) < reltol

x = copy(qc)
@time lyapcs!(acs,x,adj=true);
@test norm(acs'*x+x*acs+qc)/norm(x)/norm(acs) < reltol

x = copy(Qr)
@time lyapcs!(as,x);
@test norm(as*x+x*as'+Qr)/norm(x)/norm(as) < reltol

x = copy(Qr)
@time lyapcs!(as,x,adj=true);
@test norm(as'*x+x*as+Qr)/norm(x)/norm(as) < reltol
end
end

end

end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
using QuantumLab
using Base.Test

info("$(now())  SETTING UP VARIABLES FOR FURTHER TESTING")
indent="=> "
include("setup.jl")

# test DocumentationModule
info("$(now())  TESTING:   DocumentationModule")
docModTester() = 1
@doc """
some markdown documentation
""" docModTester
@test typeof(@doc(docModTester)) == Base.Markdown.MD
@add_doc GenericCitation("generic citation") docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@add_doc BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060") docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@add_doc JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016) docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@add_doc Citation(JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)) docModTester
@test typeof(@doc(docModTester)) == Vector{DocumentationModule.Documentation}
@test isa([Base.Markdown.parse("Hellau!"),BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060")],Vector{DocumentationModule.Documentation})
@test isa([Base.Markdown.parse("Hellau!"),JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)],Vector{DocumentationModule.Documentation})
@test isa([Base.Markdown.parse("Hellau!"),DocumentationModule.Documentation(BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"))],Vector{DocumentationModule.Documentation})
@test isa([GenericCitation("Hellau!"),DocumentationModule.Documentation(BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"))],Vector{DocumentationModule.Documentation})

info("$(now())  TESTING:   Basic Functionality")
# test AtomModule
@test Element("C") == Element("C")
el1 = Element("C")
el2 = Element("C")
@test isequal(el1,el2)
@test !isequal(el1,Element("Na"))
# test BaseModule
mqn = MQuantumNumber(0,0,0)
for mqn in MQuantumNumbers(LQuantumNumber("D")) end
@test mqn == MQuantumNumber(0,0,2)
@test distance(Position(0,0,1),Position(1,0,0)) == sqrt(2)
@test length(MQuantumNumbers(LQuantumNumber("P"))) == 3
@test 2*Position(1,2,3) == Position(1,2,3)*2
@test isless(LQuantumNumber("S"),LQuantumNumber("P"))
@test 0. ≈ BaseModule.evaluateFunction(Position(1.,2.,3.), x->x.x+x.y-x.z)
@test h2o.atoms[2].element.symbol == "O"
@test -1.4191843 ≈ readGeometryXYZ("h2o.xyz").atoms[3].position.x atol=1e-7
# test readBasisSetTX93
@test sto3g.definitions[Element("C")][1].lQuantumNumber.symbol == "S"
@test sto3g.definitions[Element("C")][1].primitives[1].exponent == 71.616837
# test BasisModule
@test -0.7332137 ≈ bas.contractedBFs[3].primitiveBFs[1].center.y atol=1e-7
@test 0.35175381 ≈ BasisFunctionsModule.evaluateFunction(origin, bas.contractedBFs[3]) atol=1e-8

# test IntegralsModule
info("$(now())  TESTING:   IntegralsModule")
@test IntegralsModule.GaussianIntegral1D_Valeev(1,.4) ≈ 0.
@test 0.0 ≈ matrixOverlap[6,1]
@test 0.3129324434238492 ≈ matrixOverlap[4,1] atol=1e-8
@test 29.00 ≈ matrixKinetic[2,2] atol=1e-2
@test 0.16175843 ≈ IntegralsModule.FIntegral(2,0.3) atol=1e-8
@test -π ≈ IntegralsModule.computeIntegralNuclearAttraction(PrimitiveGaussianBasisFunction(Position(0,0,0),1.,MQuantumNumber(1,0,0)),PrimitiveGaussianBasisFunction(Position(0,0,0),1.,MQuantumNumber(1,0,0)),Atom(Element("C"),Position(0,0,0)))
@test IntegralsModule.GaussianIntegral1D_Mathematica(6,1.2) ≈ IntegralsModule.GaussianIntegral1D_Valeev(6,1.2)
@test IntegralsModule.OverlapFundamental(bas.contractedBFs[1].primitiveBFs[1],bas.contractedBFs[1].primitiveBFs[2]) ≈ computeIntegralOverlap(bas.contractedBFs[1].primitiveBFs[1],bas.contractedBFs[1].primitiveBFs[2])

# test InitialGuessModule
info("$(now())  TESTING:   InitialGuessModule")
@test -0.264689 ≈ mean(matrixSADguess)[3,2] atol=1e-6
@test_throws ErrorException computeDensityGuessSAD("NotImplemented","STO-3G",h2o)

# test HartreeFockModule
info("$(now())  TESTING:   HartreeFock (bas, h2o, density)")
@test -74.96178985 ≈ (evaluateSCF(bas,h2o,density,info=false,detailedinfo=false)[2] + computeEnergyInteratomicRepulsion(h2o)) atol=1e-7 # checked against FermiONs++
#@test_approx_eq -4.473355520007 mean(HartreeFockModule.evaluateSCFStep(bas,h2o,mean(matrixSADguess),matrixOverlap)[1])
#@test_approx_eq_eps -74.96178985 (computeEnergyHartreeFock(bas,h2o,density) + computeEnergyInteratomicRepulsion(h2o)) 1e-7 # checked against FermiONs++

# test Shells: ShellModule and LibInt2Module
info("$(now())  TESTING:   Shells/LibInt2Shells/Basis equivalence")
shell_native  = Shell(Position(0.,0.,0.),LQuantumNumber("S"),[1.,2.,3.],[.1,.2,.3],renorm=false)
shell_libint2 = LibInt2Shell([0.,0.,0.],0,3,[1.,2.,3.],[.1,.2,.3])
shell_nativefromlibint2 = Shell(shell_libint2)
@test shell_nativefromlibint2.coefficients[2] ≈ 0.41030724 atol=1e-8
@test shell_native.coefficients[2] ≈ .2
@test computeMatrixBlockOverlap(shell_nativefromlibint2,shell_nativefromlibint2) ≈ computeMatrixBlockOverlap(shell_libint2,shell_libint2)
@test computeTensorBlockElectronRepulsionIntegrals(shell_nativefromlibint2,shell_nativefromlibint2,shell_nativefromlibint2,shell_nativefromlibint2) ≈ computeTensorBlockElectronRepulsionIntegrals(shell_libint2,shell_libint2,shell_libint2,shell_libint2)
@test LibInt2Module.nbf(shell_libint2) ≈ 1
import QuantumLab.libint2_available
if (libint2_available)
  @test_throws ErrorException LibInt2Shell([0.,0.,0.],100,1,[.1],[.5])
end

# test LibInt2Module
info("$(now())  TESTING:   LibInt2Module")
libInt2Finalize()
libInt2Initialize()
tmp_lib = LibInt2Shell([0.,1.,2.],1,2,[1.,2.],[0.5,0.5];renorm=false)
tmp = Shell(tmp_lib).coefficients
@test tmp[1] ≈ tmp[2]
destroy!(tmp_lib)
S = computeMatrixOverlap(shells)
@test S ≈ computeMatrixOverlap(shells_native)
@test S ≈ matrixOverlap
T = computeMatrixKinetic(shells)
@test T ≈ computeMatrixKinetic(shells_native)
@test T ≈ matrixKinetic
J = computeMatrixCoulomb(shells,density)
@test J ≈ computeMatrixCoulomb(bas,density) atol=1e-8
@test computeMatrixBlockOverlap(shells[1],shells[1]) ≈ S[1:1,1:1]
@test computeMatrixBlockKinetic(shells[1],shells[1]) ≈ T[1:1,1:1]
@test computeTensorBlockElectronRepulsionIntegrals(shells[1],shells[5],shells[4],shells[4])[1,1,2,2] ≈ computeElectronRepulsionIntegral(bas.contractedBFs[1],bas.contractedBFs[7],bas.contractedBFs[5],bas.contractedBFs[5])

# test LaplaceModule
info("$(now())  TESTING:   LaplaceModule")
if(!isfile("hackbusch_pretables"))
  downloadLaplacePointsHackbusch("hackbusch_pretables")
end
R = transformRangeToIdealLaplace(0.5,3.)[2]
lp = transformLaplacePointFromIdealLaplace( findLaplacePointsHackbuschPretableLarger(15,R,"hackbusch_pretables")[1], 0.5)
@test_throws ErrorException findLaplacePointsHackbuschPretableSmaller(15,R,"hackbusch_pretables")
@test LaplaceModule.computeInverseByLaplaceApproximation(2.3,lp) ≈ 1./2.3 atol=1e-7
@test LaplaceModule.computeRPADenominatorByDoubleLaplace(1.2,2.3,lp) ≈ 1./(1.2^2 + 2.3^2) atol=1e-3
println("Expecting warning below:")
@test findLaplacePointsHackbuschPretableLarger(2,100.,"hackbusch_pretables") == findLaplacePointsHackbuschPretableLarger(2,50.,"hackbusch_pretables")
@test transformLaplacePointFromIdealLaplace( findLaplacePointsHackbuschPretableSmaller(15,transformRangeToIdealLaplace(0.5,6.)[2],"hackbusch_pretables")[1], 0.5) == lp

# test RI-Module
info("$(now())  TESTING:   RIModule")
@test 0.3950752513027109 ≈ mean(computeMatrixExchangeRIK(bas,bas,matrixSADguess[1]))
tensorRICoulomb = computeTensorElectronRepulsionIntegralsRICoulomb(bas,bas)
@test 0.0429373705056905 ≈ mean(tensorRICoulomb)
@test 0.0012436634924295 ≈ tensorRICoulomb[1,3,4,5]
tensorRIOverlap = computeTensorElectronRepulsionIntegralsRIOverlap(bas,bas)
@test 0.0144597945326691 ≈ tensorRIOverlap[1,3,4,5]
@test 0.0537036598506078 ≈ mean(tensorRIOverlap)

# test SpecialMatricesModule
info("$(now())  SpecialMatricesModule")
@test -0.785008186026 ≈ mean(computeMatrixFock(bas,h2o,matrixSADguess[1])) atol=1e-10                  #### !!!! ####

### display functions ###
# These are simply here, so that the test coverage isn't limited by the display functions.
# They do not however test for anything really, so do not trust them. If you do add a *real*
# test for a display function, please add it to the corresponding module block above instead
# of here.
info("$(now())  DISPLAY FUNCTIONS")
display(bas)
display(bas.contractedBFs[3])
display(JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016))
display(GenericCitation("M. Mustermann personal note"))
display(BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"))
display([BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"), JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)])
display([Base.Markdown.parse("**Why, god, why??**"),BookCitation(["C. Darwin"], "On the Origin of Species", "978-0451529060"), JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)])
display([JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016),JournalCitation(["M. Mustermann"],"J. Stup. Mistakes",1,12,2016)])
display(@doc(IntegralsModule.GaussianIntegral1D_Valeev))
display(shells[1])
display(shell_native)
display(lp)
summarize(BasisModule.Basis)
summarize(GaussianBasis)

info("\n\n     CONGRATULATIONS!! ALL TESTS COMPLETED!")


text:
export EQMag

"""
    EQMag

Earthquake magnitude container object

Field     | Type    | Meaning
--------: |:------- |:--------------
 val      | Float32 | numeric magnitude value (note: Float32!)
scale     | String  | magnitude scale (freeform)
 gap      | Float64 | azimuthal gap (°)
 nst      | Int64   | number of stations used in magnitude calculation
 src      | String  | magnitude source
"""
mutable struct EQMag
  val     ::Float32
  scale   ::String
  nst     ::Int64
  gap     ::Float64
  src     ::String

  function EQMag( val     ::Float32,
                  scale   ::String,
                  nst     ::Int64,
                  gap     ::Float64,
                  src     ::String
                )
    return new(val, scale, nst, gap, src)
  end
end
EQMag(;
      val     ::Float32   =   -5.0f0,
      scale   ::String    =   "",
      nst     ::Int64     =   zero(Int64),
      gap     ::Float64   =   zero(Float64),
      src     ::String    =   ""
      ) = EQMag(val, scale, nst, gap, src)

isempty(Mag::EQMag) = min(getfield(Mag, :val) == -5.0f0,
                          getfield(Mag, :gap) == zero(Float64),
                          getfield(Mag, :nst) == zero(Int64),
                          isempty(getfield(Mag, :scale)),
                          isempty(getfield(Mag, :src)))

function hash(Mag::EQMag)
  h = hash(getfield(Mag, :val))
  for f in (:scale, :nst, :gap, :src)
    h = hash(getfield(Mag, f), h)
  end
  return h
end

function isequal(S::EQMag, U::EQMag)
  q::Bool = isequal(getfield(S, :val), getfield(U, :val))
  for f in (:scale, :nst, :gap, :src)
    q = min(q, getfield(S,f) == getfield(U,f))
  end
  return q
end
==(S::EQMag, U::EQMag) = isequal(S, U)

sizeof(Mag::EQMag) = 52 + sizeof(Mag.src) + sizeof(Mag.scale)

function write(io::IO, M::EQMag)
  write(io, getfield(M, :val))
  write(io, getfield(M, :gap))
  write(io, getfield(M, :nst))

  scale = codeunits(getfield(M, :scale))
  write(io, Int64(length(scale)))
  write(io, scale)

  src = codeunits(getfield(M, :src))
  write(io, Int64(length(src)))
  write(io, src)
  return nothing
end

function read(io::IO, ::Type{EQMag})
  M = EQMag()
  setfield!(M, :val, fastread(io, Float32))
  setfield!(M, :gap, fastread(io, Float64))
  setfield!(M, :nst, fastread(io, Int64))

  L = fastread(io, Int64)
  setfield!(M, :scale, String(fastread(io, L)))

  L = fastread(io, Int64)
  setfield!(M, :src, String(fastread(io, L)))
  return M
end

function show(io::IO, Mag::EQMag)
  if get(io, :compact, false) == false
    println(io, "EQMag with fields:")
    for f in (:val, :scale, :nst, :gap, :src)
      fn = lpad(String(f), 5, " ")
      println(io, fn, ": ", getfield(Mag, f))
    end
  else
    c = :compact => true
    print(io, getfield(Mag, :scale), " ",
              repr(getfield(Mag, :val), context=c), " ",
              "(g ", repr(getfield(Mag, :gap), context=c), "°, ",
              "n ", getfield(Mag, :nst), ")")
  end
  return nothing
end


text:
module KitMSA

using DataFrames
using DelimitedFiles
using FastaIO
using StatsBase
using Random
using GZip

include("clean_MSA.jl")
include("read_write_MSA.jl")
include("BM_utils.jl")
include("contact_prediction.jl")

cod2amino = Dict( "ATA" => Int8(8), "ATC" => Int8(8), "ATT"=> Int8(8), "ATG"=> Int8(11),
        "ACA"=>Int8(17), "ACC"=>Int8(17), "ACG"=>Int8(17), "ACT"=> Int8(17),
        "AAC"=>Int8(12), "AAT"=>Int8(12), "AAA"=>Int8(9), "AAG"=>Int8(9),
        "AGC"=>Int8(16), "AGT"=> Int8(16), "AGA"=> Int8(15), "AGG"=> Int8(15),
        "CTA"=>Int8(10), "CTC"=>Int8(10), "CTG"=>Int8(10), "CTT"=>Int8(10),
        "CCA"=>Int8(13), "CCC"=>Int8(13), "CCG"=>Int8(13), "CCT"=>Int8(13),
        "CAC"=>Int8(7), "CAT"=>Int8(7), "CAA"=>Int8(14), "CAG"=>Int8(14),
        "CGA"=>Int8(15), "CGC"=>Int8(15), "CGG"=>Int8(15), "CGT"=>Int8(15),
        "GTA"=>Int8(18), "GTC"=>Int8(18), "GTG"=>Int8(18), "GTT"=>Int8(18),
        "GCA"=>Int8(1), "GCC"=>Int8(1), "GCG"=>Int8(1), "GCT"=>Int8(1),
        "GAC"=>Int8(3), "GAT"=>Int8(3), "GAA"=>Int8(4), "GAG"=>Int8(4),
        "GGA"=>Int8(6), "GGC"=>Int8(6), "GGG"=>Int8(6), "GGT"=>Int8(6),
        "TCA"=>Int8(16), "TCC"=>Int8(16), "TCG"=>Int8(16), "TCT"=>Int8(16),
        "TTC"=>Int8(5), "TTT"=>Int8(5), "TTA"=>Int8(10), "TTG"=>Int8(10),
        "TAC"=>Int8(20), "TAT"=>Int8(20), "TAA"=> Int8(21), "TAG"=> Int8(21),
        "TGC"=> Int8(2), "TGT"=>Int8(2) , "TGA"=> Int8(21), "TGG"=> Int8(19),
        "---" => Int8(21) )

amino2cod = Dict()
for amino in 1:21
    codons = Array{String, 1}()
    for (key, val) in cod2amino
       val == amino && push!(codons, key)
    end
    amino2cod[amino] = codons
end

export cod2amino, amino2cod

end

max_stars_repo_name: matteobisardi/KitMSA

max_issues_repo_name: matteobisardi/KitMSA

max_forks_repo_name: matteobisardi/KitMSA


text:
function euler(f::Function, x::Float64, y::Vector{Float64}, steps::Int64, h::Float64)
  nsteps = steps + 1
  res = zeros(nsteps, length(y)+1)
  for i in 1:nsteps
    res[i,:] = [x; [y[i] for i in 1:length(y)]]
    k0 = h * f(x, y)
    y += k0
    x += h
  end
  res
end

function modified_euler(f::Function, x::Float64, y::Vector{Float64}, steps::Int64, h::Float64)
  nsteps = steps + 1
  res = zeros(nsteps, length(y)+1)
  for i in 1:nsteps
    res[i,:] = [x; [y[i] for i in 1:length(y)]]
    k0 = h * f(x, y)
    k1 = h * f(x+h, y+k0)
    y += (k0+k1)/2.0
    x += h
  end
  res
end

function mid_point_euler(f::Function, x::Float64, y::Vector{Float64}, steps::Int64, h::Float64)
  nsteps = steps + 1
  res = zeros(nsteps, length(y)+1)
  for i in 1:nsteps
    res[i,:] = [x; [y[i] for i in 1:length(y)]]
    k0 = h * f(x, y)
    k1 = h * f(x+h/2.0, y+k0/2.0)
    y += k1
    x += h
  end
  res
end

function runga_kutta_4(f::Function, x::Float64, y::Vector{Float64}, steps::Int64, h::Float64)
  nsteps = steps + 1
  res = zeros(nsteps, length(y)+1)
  for i in 1:nsteps
    res[i,:] = [x; [y[i] for i in 1:length(y)]]
    k0 = h * f(x, y)
    k1 = h * f(x+h/2.0, y+k0/2.0)
    k2 = h * f(x+h/2.0, y+k1/2.0)
    k3 = h * f(x+h/2.0, y+k2)
    y += (k0+2.0*k1+2.0*k2+k3)/6.0
    x += h
  end
  res
end


text:
# solve the ground state of Lieb-Liniger model
# also demonstrates the usage of preconditioners in Optim.jl

using TensorKit
using TensorOperations
using TensorKitAD
using KrylovKit
using Zygote
using KrylovKit
using LinearAlgebra
using Plots
using Optim
using ChainRules
using ChainRulesCore

using Revise
using statmech_tm_solver

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

# parameter for lieb-liniger
c = 1
μ = 1
L = 50

# bond dimensions
d, χ = 1, 8

# random initial guess
ψ_arr = rand(ComplexF64, (χ, d+1, χ))
convert_to_cmps(ψ_arr)

# cost function (energy) and its gradient
function f(ψ_arr::Array{ComplexF64, 3})
    ψ = convert_to_cmps(ψ_arr)
    return energy_lieb_liniger(ψ, c, L, μ)
end

function g!(ψg_arr::Array{ComplexF64, 3}, ψ_arr::Array{ComplexF64, 3})
    copyto!(ψg_arr, f'(ψ_arr))
end

# use LBFGS to optimize the energy for 250 steps
res = optimize(f, g!, ψ_arr, LBFGS(), Optim.Options(show_trace=true, iterations=250))
ψm = convert_to_cmps(Optim.minimizer(res));
ψm_arr = convert_to_array(ψm);

# calculate observables. total, kinetic and potential energies
E = energy_lieb_liniger(ψm, c, L, μ)
env = finite_env(K_mat(ψm, ψm), L)
env = permute(env, (2, 3), (4, 1))

op_ρ = particle_density(ψm)
op_k = kinetic(ψm)
tr(env * op_ρ)
tr(env * op_k)

#########################################################
# now we test the performance of the preconditioner

# construct the preconditioner
function _precondprep!(P::preconditioner, ψ_arr::Array{ComplexF64, 3})
    ψ = convert_to_cmps(ψ_arr)
    P.map = tangent_map(ψ, L)
    P.proj = gauge_fixing_proj(ψ, L)
end

# use preconditioned LBFGS, optimize the energy for additional 50 steps
P0 = preconditioner(ψm, L)
res1 = optimize(f, g!, ψm_arr, LBFGS(P = P0, precondprep=_precondprep!), Optim.Options(show_trace=true, iterations=50))
ψm1 = convert_to_cmps(Optim.minimizer(res1));
E1 = energy_lieb_liniger(ψm1, c, L, μ)
@show E1 - E

# use ordinary LBFGS, optimize the energy for additional 50 steps
res2 = optimize(f, g!, ψm_arr, LBFGS(), Optim.Options(show_trace=true, iterations=50))
ψm2 = convert_to_cmps(Optim.minimizer(res2));
E2 = energy_lieb_liniger(ψm2, c, L, μ)
@show E2 - E

# from original random init guess, 100 steps of ordinary LBFGS + 100 steps of preconditioned GradientDescent + 100 steps of preconditioned LBFGS
res3 = optimize(f, g!, ψ_arr, LBFGS(), Optim.Options(show_trace=true, iterations=100))
ψ3_arr = Optim.minimizer(res3)
ψm3 = convert_to_cmps(ψ3_arr)

P0 = preconditioner(ψm3, L)
res3 = optimize(f, g!, ψ3_arr, GradientDescent(P = P0, precondprep=_precondprep!), Optim.Options(show_trace=true, iterations=100))
ψ3_arr = Optim.minimizer(res3)
ψm3 = convert_to_cmps(ψ3_arr)

P0 = preconditioner(ψm3, L)
res3 = optimize(f, g!, ψ3_arr, LBFGS(P = P0, precondprep=_precondprep!), Optim.Options(show_trace=true, iterations=100))
ψm3 = convert_to_cmps(Optim.minimizer(res3));

E3 = energy_lieb_liniger(ψm3, c, L, μ)
@show E3 - E