text:
@testset "Samplers" begin

    mt_String_array = groupmod(["b", "a", "c","e","d"])
    rng = MersenneTwister()

    design1 =constructdesign(mt_String_array,3,order=true)
    test_space = getspace(design1)
    n_samples = rand(1:20)
    a_sample = sample(rng,test_space,n_samples, with_index = false)

    @testset "sample_with_index = false " begin
        @test length(a_sample) == n_samples
        @test_throws BoundsError size(a_sample)[2]
        #number of samples oke
        @test  length(a_sample) == n_samples
        # all unique
        @test  length(Set(a_sample)) == n_samples
        # collect works
    end

    a_sample_with_index = sample(rng,test_space,n_samples, with_index = true)

    @testset "sample_with_index = true " begin
        @test length(a_sample_with_index[1]) == n_samples
        @test length(a_sample_with_index) == 2
        #number of samples oke
        @test  isequal(test_space[a_sample_with_index[2][1]],a_sample_with_index[1][1]) == true
        # collect works
    end

    Unordered_con_c_e = UnOrderedConstraint([Mod("c") , Mod("e")])
    design2 = constructdesign(mt_String_array,3,Unordered_con_c_e,order=true)
    test_space2 = getspace(design2)
    n_samples2 = rand(1:10)
    a_sample2 = sample(rng,test_space2,n_samples2)

    @testset "sample_reject" begin
        # number of samples
        @test length(a_sample2[1]) == n_samples2
        # all unique elements
        @test a_sample2[1] |> unique |> length == n_samples2
        #contrain is used
        for i in a_sample2
            @test (Mod("c") in i && Mod("e") in i) == false
        end
        @test typeof(a_sample2[2]) == Vector{Int}
    end

    design3  = constructdesign(mt_String_array,[2,3],UnOrderedConstraint([Mod("c") , Mod("e")]),order=true)
    test_space3 = getspace(design3)
    n_samples3 = rand(1:10)
    a_sample3 = sample(rng,test_space3,n_samples3)

    @testset "multiframe" begin
        # number of samples
        @test length(a_sample3[1]) == n_samples3
        # all unique elements
        @test a_sample3[1] |> unique |> length == n_samples3
        #contrain is used
        for i in a_sample3
            @test (Mod("c") in i && Mod("e") in i) == false
        end
    end

end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
using Base: @deprecate_binding

# Deprecated as of v0.3.0
@deprecate_binding LZ4Compressor LZ4FrameCompressor
@deprecate_binding LZ4Decompressor LZ4FrameDecompressor
@deprecate_binding LZ4CompressorStream LZ4FrameCompressorStream
@deprecate_binding LZ4DecompressorStream LZ4FrameDecompressorStream

ext: jl


text:
using NPZ
using CSV
using DataFrames
using Statistics
using CartesianGeneticProgramming
using Cambrian
import Cambrian.mutate
using Plots
include("utils.jl")

struct Month
    topo0
    topo1
    tp
    hs
    direction
end

dataset_path = get_dataset_path()
n_months = 236
months = []

for i in 1:n_months
    topo_path = string(dataset_path, '/', i, "_topo.npy")
    topo = npzread(topo_path)
    topo0 = norm_to_range(topo[1, :], -10.0, 10.0, -1, 1)
    topo1 = norm_to_range(topo[2, :], -10.0, 10.0, -1, 1)

    waves_path = string(dataset_path, '/', i, "_waves.csv")
    waves_csv = CSV.read(waves_path, DataFrame)

    # t = min_max_norm(waves_csv.t)
    tp_month = norm_to_range(waves_csv.tp, 0, 25, -1, 1)
    hs_month = norm_to_range(waves_csv.hs, 0, 10, -1, 1)
    direction_month = norm_to_range(waves_csv.direction, 0, 360, -1, 1)

    tp_days = []
    hs_days = []
    direction_days = []

    for i in range(0, stop=27)
        tp = tp_month[(i*8) + 1 : (i*8) + 8]
        hs = hs_month[(i*8) + 1 : (i*8) + 8]
        direction = direction_month[(i*8) + 1 : (i*8) + 8]
        push!(tp_days, tp)
        push!(hs_days, hs)
        push!(direction_days, direction)
    end
    month = Month(topo0, topo1, tp_days, hs_days, direction_days)
    global months
    push!(months, month)
end

function evaluate(ind::CGPInd, months::Vector{Any})
    errors = Vector{Float64}(undef, size(months)[1])

    for m_i in 1:size(months)[1]
        month = months[m_i]
        topo0 = deepcopy(month.topo0)
        for w_i in 1:size(month.tp)[1]
            x_row = vcat(topo0, month.tp[w_i], month.hs[w_i], month.direction[w_i])
            topo0 = process(ind, x_row)
        end
        errors[m_i] = mean(abs.(month.topo1 - topo0))
    end
    [mean(errors) * -1]
end

cfg = get_config("topo.yaml")
fit(i::CGPInd) = evaluate(i, months)
mutate(i::CGPInd) = goldman_mutate(cfg, i)
e = CGPEvolution(cfg, fit)
run!(e)

max_stars_repo_name: mahmoud-al-najar/cgp-morphcast

max_issues_repo_name: mahmoud-al-najar/cgp-morphcast

max_forks_repo_name: mahmoud-al-najar/cgp-morphcast


text:
# Test of connected components

using Graphs
using Base.Test

g = simple_adjlist(8, is_directed=false)

eds = [(1, 2), (1, 3), (2, 4), (3, 4), (5, 6), (6, 7), (7, 5)]

for i = 1 : length(eds)
    e = eds[i]
    add_edge!(g, e[1], e[2])
end

ccs = connected_components(g)

@test length(ccs) == 3

@test ccs[1] == [1, 2, 3, 4]
@test ccs[2] == [5, 6, 7]
@test ccs[3] == [8]

###########################################################
#
#   Connected components of directed graph
#
##########################################################

# test 1
eds = [(1, 2), (2, 3), (3, 1), (4, 1)]
g = simple_graph(4)

for (u, v) in eds
    add_edge!(g, u, v)
end

scc = strongly_connected_components(g)
@test length(scc) == 2
@test scc[1] == [1, 2, 3]
@test scc[2] == [4]

# test 2, from Vazirani's notes
# (http://www.cs.berkeley.edu/~vazirani/s99cs170/notes/lec12.pdf)

eds = [(1, 2), (2, 3), (2, 4), (2, 5), (3, 6),
       (4, 5), (4, 7), (5, 2), (5, 6), (5, 7),
       (6, 3), (6, 8), (7, 8), (7, 10), (8, 7),
       (9, 7), (10, 9), (10, 11), (11, 12), (12, 10)]
g = simple_graph(12)

for (u, v) in eds
    add_edge!(g, u, v)
end

scc = strongly_connected_components(g)
@test length(scc) == 4
@test scc[1] == [8, 7, 10, 9, 11, 12]
@test scc[2] == [3, 6]
@test scc[3] == [2, 4, 5]
@test scc[4] == [1]


text:
using ..Utils: LDS, CLCDS

using LazySets: CachedMinkowskiSumArray,
                isdisjoint

import LazySets.Approximations: overapproximate

"""
combine_cpas(cpa1::CartesianProductArray{N, S}, cpa2::CartesianProductArray{N, S},
                       blocks1::Vector{Int}, blocks2::Vector{Int}) where {N, S<:LazySet{N}}

Compose high-dimensional cartesian product array from 2 low-dimensional CPA's

### Input

- `cpa1`   -- Cartesian Product Array
- `cpa2`   -- Cartesian Product Array
- `blocks1`   -- Block structure of cpa1 in sense of high-dimensional set
- `blocks2`   -- Block structure of cpa2 in sense of high-dimensional set

### Output

High-dimensional cartesian product array.
"""
function combine_cpas(cpa1::CartesianProductArray{N, S}, cpa2::CartesianProductArray{N, S},
                       blocks1::Vector{Int}, blocks2::Vector{Int}) where {N, S<:LazySet{N}}
    result = Vector{S}(undef, length(array(cpa1)) + length(array(cpa2)))
    result[blocks1] = array(cpa1)
    result[blocks2] = array(cpa2)
    return CartesianProductArray(result)
end

"""
    reach(problem, options)

Interface to reachability algorithms for an LTI system.

### Input

- `problem`   -- initial value problem
- `options`   -- additional options

### Output

A sequence of [`SparseReachSet`](@ref)s.

### Notes

The numeric type of the system's coefficients and the set of initial states
is inferred from the first parameter of the system (resp. lazy set), ie.
`NUM = first(typeof(problem.s).parameters)`.
"""
function reach_mixed(problem::Union{IVP{<:CLDS{NUM}, <:LazySet{NUM}},
                              IVP{<:CLCDS{NUM}, <:LazySet{NUM}}},
               options::TwoLayerOptions
              ) where {NUM <: Real}
    # optional matrix conversion
    problem = matrix_conversion(problem, options)

    # list containing the arguments passed to any reachability function
    args = []

    # coefficients matrix
    A = problem.s.A
    push!(args, A)

    # determine analysis mode (sparse/dense) for lazy_expm mode
    if A isa SparseMatrixExp
        push!(args, Val(options[:assume_sparse]))
    end

    n = statedim(problem)
    blocks = options[:blocks]
    partition = convert_partition(options[:partition])
    T = options[:T]
    N = (T == Inf) ? nothing : ceil(Int, T / options[:δ])
    # Cartesian decomposition of the initial set
    if length(partition) == 1 && length(partition[1]) == n &&
            options[:block_options_init] == LinearMap
        info("- Skipping decomposition of X0")
        Xhat0 = LazySet{NUM}[problem.x0]
    else
        info("- Decomposing X0")
        @timing begin
            Xhat0 = array(decompose(problem.x0, options[:partition],
                                    options[:block_options_init]))
        end
    end

    # compute dimensions
    dimensions = compute_dimensions(partition, blocks)

    # determine output function: linear map with the given matrix
    output_function = options[:output_function] != nothing ?
        (x -> options[:output_function] * x) :
        nothing

    # preallocate output vector
    if output_function == nothing
        res_type = SparseReachSet{CartesianProductArray{NUM, LazySet{NUM}}}
    else
        res_type = SparseReachSet{Hyperrectangle{NUM}}
    end
    res = (N == nothing) ? Vector{res_type}() : Vector{res_type}(undef, N)

    # shortcut if only the initial set is required
    if N == 1
        res[1] = res_type(
            CartesianProductArray{NUM, LazySet{NUM}}(Xhat0[blocks]),
            zero(NUM), options[:δ])
        return Flowpipe(res)
    end
    push!(args, Xhat0)

    # inputs
    if !options[:assume_homogeneous] && inputdim(problem) > 0
        U = inputset(problem)
    else
        U = nothing
    end
    push!(args, U)

    # overapproximation function for states
    block_options_iter = options[:block_options_iter]
    if block_options_iter isa AbstractVector ||
            block_options_iter isa Dict{Int, Any}
        # individual overapproximation options per block
        overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter[i])
    else
        # uniform overapproximation options for each block
        overapproximate_fun = (i, X) -> overapproximate(X, block_options_iter)
    end
    push!(args, overapproximate_fun)

    # overapproximate function for inputs
    lazy_inputs_interval = options[:lazy_inputs_interval]
    if lazy_inputs_interval == lazy_inputs_interval_always
        overapproximate_inputs_fun = (k, i, x) -> overapproximate_fun(i, x)
    else
        # first set in a series
        function _f(k, i, x::LinearMap{NUM}) where {NUM}
            @assert k == 1 "a LinearMap is only expected in the first iteration"
            return CachedMinkowskiSumArray(LazySet{NUM}[x])
        end
        # further sets of the series
        function _f(k, i, x::MinkowskiSum{NUM, <:CachedMinkowskiSumArray}) where NUM
            if has_constant_directions(block_options_iter, i)
                # forget sets if we do not use epsilon-close approximation
                forget_sets!(x.X)
            end
            push!(array(x.X), x.Y)
            if lazy_inputs_interval(k)
                # overapproximate lazy set
                y = overapproximate_fun(i, x.X)
                return CachedMinkowskiSumArray(LazySet{NUM}[y])
            end
            return x.X
        end
        function _f(k, i, x)
            # other set types
            if lazy_inputs_interval(k)
                # overapproximate lazy set
                return overapproximate_fun(i, x.X)
            end
            return x
        end
        overapproximate_inputs_fun = _f
    end
    push!(args, overapproximate_inputs_fun)

    # ambient dimension
    push!(args, n)

    # number of computed sets
    push!(args, N)

    # output function
    push!(args, output_function)

    # add mode-specific block(s) argument
    push!(args, blocks)
    push!(args, partition)
    push!(args, dimensions)

    # time step
    push!(args, options[:δ])

    # termination function
    invariant = stateset(problem.s)
    vars = [e for block in blocks for e in partition[block]]
    if invariant isa Universe
        invariant = Universe(length(vars))
    else
        invariant = LazySets.Approximations.project(invariant, vars)
    end
    termination = get_termination_function_out(N, invariant)
    if length(options[:guards_proj]) > 0
        push!(args, options[:guards_proj])
    else
        push!(args, [Universe{NUM}(length(vars))])
    end
    push!(args, options[:block_options_iter])
    push!(args, vars)
    push!(args, termination)

    info("- Computing successors")

    # progress meter
    progress_meter = (N != nothing) ?
        Progress(N, 1, "Computing successors ") :
        nothing
    push!(args, progress_meter)

    # preallocated result
    push!(args, res)

    # call the adequate function with the given arguments list
    @timing begin
        index, skip = reach_blocks!(args...)
    end

    # shrink result array
    if skip || index < N
        if N != nothing
            info("termination after only $index of $N steps")
        end
        deleteat!(res, (skip ? index : index + 1):length(res))
    end

    # return the result
    return Flowpipe(res)
end

function reach_mixed(problem::Union{IVP{<:CLCS{NUM}, <:LazySet{NUM}},
                              IVP{<:CLCCS{NUM}, <:LazySet{NUM}}},
               options::TwoLayerOptions
              ) where {NUM <: Real}
    info("Time discretization...")
    Δ = @timing discretize(problem, options[:δ],
        algorithm=options[:discretization], exp_method=options[:exp_method],
        sih_method=options[:sih_method])
    return reach_mixed(Δ, options)
end

function get_termination_function_out(N::Nothing, invariant::Universe)
    termination_function = (k, set, t0) -> termination_unconditioned(set)
    warn("no termination condition specified; the reachability analysis will " *
         "not terminate")
    return termination_function
end

function get_termination_function_out(N::Int, invariant::Universe)
    return (k, set, t0) -> termination_N_out(N, k, t0, set)
end

function get_termination_function_out(N::Nothing, invariant::LazySet)
    return (k, set, t0) -> termination_inv_out(invariant, set, t0)
end

function get_termination_function_out(N::Int, invariant::LazySet)
    return (k, set, t0) -> termination_inv_N_out(N, invariant, k, set, t0)
end

function termination_unconditioned_out(set)
    return (false, false, set)
end

function termination_N_out(N, k, t0, set)
    return (k >= N, false, set)
end

function termination_inv_out(inv, set, t0)
    l_inters = Intersection(set, inv)
    if isempty(l_inters)
        return (true, true, EmptySet(dim(inv)))
    else
        return (false, false, l_inters)
    end
end

function termination_inv_N_out(N, inv, k, set, t0)
    if k >= N
        return (true, false, EmptySet(dim(inv)))
    end
    l_inters = Intersection(set, inv)
    if isempty(l_inters)
        return (true, true, EmptySet(dim(inv)))
    else
        return (false, false, l_inters)
    end
end


text:
"""
    MetadataArray(parent::AbstractArray, metadata)

Custom `AbstractArray` object to store an `AbstractArray` `parent` as well as some `metadata`.

# Examples

```jldoctest metadataarray
julia> v = ["John", "John", "Jane", "Louise"];

julia> s = MetadataArray(v, Dict("John" => "Treatment", "Louise" => "Placebo", "Jane" => "Placebo"))
4-element MetadataArrays.MetadataArray{String,Dict{String,String},1,Array{String,1}}:
 "John"
 "John"
 "Jane"
 "Louise"

julia> metadata(s)
Dict{String,String} with 3 entries:
  "John"   => "Treatment"
  "Jane"   => "Placebo"
  "Louise" => "Placebo"
```
"""
struct MetadataArray{T, N, M, S<:AbstractArray} <: AbstractArray{T, N}
    parent::S
    metadata::M
end

MetadataArray(v::AbstractArray{T, N}, m::M = ()) where {T, N, M} =
     MetadataArray{T, N, M, typeof(v)}(v, m)

"""
    MetadataVector{T, M, S<:AbstractArray}

Shorthand for `MetadataArray{T, 1, M, S}`.
"""
const MetadataVector{T, M, S<:AbstractArray} = MetadataArray{T, 1, M, S}

MetadataVector(v::AbstractVector, n = ()) = MetadataArray(v, n)

Base.size(s::MetadataArray) = Base.size(parent(s))

Base.axes(s::MetadataArray) = Base.axes(parent(s))

Base.IndexStyle(T::Type{<:MetadataArray}) = IndexStyle(_parent_type(T))

Base.@propagate_inbounds function Base.getindex(s::MetadataArray, x::Int...)
    @boundscheck checkbounds(parent(s), x...)
    @inbounds ret = getindex(parent(s), x...)
    return ret
end

Base.@propagate_inbounds function Base.setindex!(s::MetadataArray, val, x::Int...)
    @boundscheck checkbounds(parent(s), x...)
    @inbounds parent(s)[x...] = val
    return val
end

Base.parent(s::MetadataArray) = s.parent

_parent_type(::Type{MetadataArray{T, M, N, S}}) where {T,M,N,S} = S

"""
    metadata(s::MetadataArray)

Returns metadata for `s`.
"""
metadata(s::MetadataArray) = s.metadata

metadata(s::SubArray) = metadata(parent(s))

metadata(s::T) where {T<:AbstractArray} = nothing

Base.similar(A::MetadataArray, ::Type{S}, dims::Dims) where S =
    MetadataArray(similar(parent(A), S, dims), metadata(A))

function Base.reshape(s::MetadataArray, d::Dims)
    MetadataArray(reshape(parent(s), d), metadata(s))
end


text:
module PCurlGradMonomialBasesTests

using Test
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials
using Gridap.Arrays

xi = Point(4,2)
np = 1
x = fill(xi,np)

order = 2
D = 2
T = Float64
V = VectorValue{D,T}
G = gradient_type(V,xi)
b = PCurlGradMonomialBasis{D}(T,order)

v = V[
  (1.0, 0.0), (0.0, 1.0), (4.0, 0.0), (0.0, 2.0), (16.0, 0.0), (0.0, 4.0),
  (2.0, 0.0), (0.0, 4.0), (8.0, 0.0), (0.0, 8.0), (4.0, 0.0), (0.0, 16.0),
  (64.0, 32.0), (32.0, 16.0), (16.0, 8.0)]

g = G[
  (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0), (1.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 1.0), (8.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 4.0),
  (0.0, 1.0, 0.0, 0.0), (0.0, 0.0, 1.0, 0.0), (2.0, 4.0, 0.0, 0.0),
  (0.0, 0.0, 2.0, 4.0), (0.0, 4.0, 0.0, 0.0), (0.0, 0.0, 8.0, 0.0),
  (48.0, 0.0, 16.0, 16.0), (16.0, 16.0, 4.0, 16.0), (4.0, 16.0, 0.0, 12.0)]

  vb = evaluate(b,x)

  for (vi,vbi) in zip(v,vb)
    @test vi == vbi
  end

  ∇b = Broadcasting(gradient)(b)
  gvb = evaluate(∇b,x)
  for (vi,vbi) in zip(g,gvb)
    @test vi == vbi
  end

  @test num_terms(b) == 15
  @test get_order(b) == 2

  xi = Point(2,3,5)
  np = 5
  x = fill(xi,np)

  order = 1
  D = 3
  T = Float64
  V = VectorValue{D,T}
  G = gradient_type(V,xi)
  b = PCurlGradMonomialBasis{D}(T,order)

  @test num_terms(b) == 15
  @test get_order(b) == 1

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/cholmod.jl, to deprecate
@eval SuiteSparse.CHOLMOD begin
    using Base.LinAlg: Adjoint, Transpose
    Base.Ac_ldiv_B(A::RealHermSymComplexHermF64SSL, B::StridedVecOrMat) = \(Adjoint(A), B)
    Base.Ac_ldiv_B(L::Factor, B::Dense) = \(Adjoint(L), B)
    Base.Ac_ldiv_B(L::Factor, B::VecOrMat) = \(Adjoint(L), B)
    Base.Ac_ldiv_B(L::Factor, B::Sparse) = \(Adjoint(L), B)
    Base.Ac_ldiv_B(L::Factor, B::SparseVecOrMat) = \(Adjoint(L), B)
    Base.Ac_ldiv_B(L::FactorComponent, B) = \(Adjoint(L), B)
    Base.Ac_ldiv_B(L::FactorComponent, B::RowVector) = \(Adjoint(L), B)
    Base.Ac_mul_B(A::Sparse, B::Dense) = *(Adjoint(A), B)
    Base.Ac_mul_B(A::Sparse, B::VecOrMat) =  *(Adjoint(A), B)
    Base.Ac_mul_B(A::Sparse, B::Sparse) = *(Adjoint(A), B)
    Base.A_mul_Bc(A::Sparse{Tv}, B::Sparse{Tv}) where {Tv<:VRealTypes} = *(A, Adjoint(B))
end

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/umfpack.jl, to deprecate
@eval SuiteSparse.UMFPACK begin
    using Base.LinAlg: Adjoint, Transpose
    Base.A_ldiv_B!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
        Base.LinAlg.ldiv!(X, lu, B)
    Base.At_ldiv_B!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
        Base.LinAlg.ldiv!(X, Transpose(lu), B)
    Base.Ac_ldiv_B!(X::StridedVecOrMat{T}, lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} =
        Base.LinAlg.ldiv!(X, Adjoint(lu), B)
    Base.A_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
        Base.LinAlg.ldiv!(X, lu, B)
    Base.At_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
        Base.LinAlg.ldiv!(X, Transpose(lu), B)
    Base.Ac_ldiv_B!(X::StridedVecOrMat{Tb}, lu::UmfpackLU{Float64}, B::StridedVecOrMat{Tb}) where {Tb<:Complex} =
        Base.LinAlg.ldiv!(X, Adjoint(lu), B)
    Base.A_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(lu, B)
    Base.At_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(Transpose(lu), B)
    Base.Ac_ldiv_B!(lu::UmfpackLU{T}, B::StridedVecOrMat{T}) where {T<:UMFVTypes} = Base.LinAlg.ldiv!(Adjoint(lu), B)
    Base.A_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(lu, B)
    Base.At_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(Transpose(lu), B)
    Base.Ac_ldiv_B!(lu::UmfpackLU{Float64}, B::StridedVecOrMat{<:Complex}) = Base.LinAlg.ldiv!(Adjoint(lu), B)
end

# A[ct]_(mul|ldiv|rdiv)_B[ct][!] methods from src/spqr.jl, to deprecate
@eval SuiteSparse.SPQR begin
    using Base.LinAlg: Adjoint, Transpose
    Base.A_mul_Bc!(A::StridedMatrix, Q::QRSparseQ) = Base.LinAlg.mul!(A, Adjoint(Q))
    Base.Ac_mul_B!(Q::QRSparseQ, A::StridedVecOrMat) = Base.LinAlg.mul!(Adjoint(Q), A)
    Base.A_mul_B!(A::StridedMatrix, Q::QRSparseQ) = Base.LinAlg.mul!(A, Q)
    Base.A_mul_B!(Q::QRSparseQ, A::StridedVecOrMat) = Base.LinAlg.mul!(Q, A)
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
import Base: iterate

#Simple methods
export powm, powm!, invpowm, invpowm!

mutable struct PowerMethodIterable{matT, vecT <: AbstractVector, numT <: Number, eigvalT <: Number}
    A::matT
    x::vecT
    tol::numT
    maxiter::Int
    θ::eigvalT
    r::vecT
    Ax::vecT
    residual::numT
end

##
## Iterators
##

@inline converged(p::PowerMethodIterable) = p.residual ≤ p.tol

@inline start(p::PowerMethodIterable) = 0

@inline done(p::PowerMethodIterable, iteration::Int) = iteration > p.maxiter || converged(p)

function iterate(p::PowerMethodIterable, iteration::Int=start(p))
    if done(p, iteration) return nothing end

    mul!(p.Ax, p.A, p.x)

    # Rayleigh quotient θ = x'Ax
    p.θ = dot(p.x, p.Ax)

    # (Previous) residual vector r = Ax - λx
    copyto!(p.r, p.Ax)
    BLAS.axpy!(-p.θ, p.x, p.r)

    # Normed residual
    p.residual = norm(p.r)

    # Normalize the next approximation
    copyto!(p.x, p.Ax)
    rmul!(p.x, one(eltype(p.x)) / norm(p.x))

    p.residual, iteration + 1
end

# Transforms the eigenvalue back whether shifted or inversed
@inline transform_eigenvalue(θ, inverse::Bool, σ) = σ + (inverse ? inv(θ) : θ)

function powm_iterable!(A, x; tol = eps(real(eltype(A))) * size(A, 2) ^ 3, maxiter = size(A, 1))
    T = eltype(x)
    PowerMethodIterable(A, x, tol, maxiter, zero(T), similar(x), similar(x), floatmax(real(T)))
end

"""
    powm(B; kwargs...) -> λ, x, [history]

See [`powm!`](@ref). Calls `powm!(B, x0; kwargs...)` with
`x0` initialized as a random, complex unit vector.
"""
function powm(B; kwargs...)
    x0 = rand(Complex{real(eltype(B))}, size(B, 1))
    rmul!(x0, one(eltype(B)) / norm(x0))
    powm!(B, x0; kwargs...)
end

"""
    powm!(B, x; shift = zero(eltype(B)), inverse::Bool = false, kwargs...) -> λ, x, [history]

By default finds the approximate eigenpair `(λ, x)` of `B` where `|λ|` is largest.

# Arguments
- `B`: linear map, see the note below.
- `x`: normalized initial guess. Don't forget to use complex arithmetic when necessary.

## Keywords
- `tol::Real = eps(real(eltype(B))) * size(B, 2) ^ 3`: stopping tolerance for the residual norm;
- `maxiter::Integer = size(B,2)`: maximum number of iterations;
- `log::Bool`: keep track of the residual norm in each iteration;
- `verbose::Bool`: print convergence information during the iterations.

!!! note "Shift-and-invert"
    When applying shift-and-invert to ``Ax = λx`` with `invert = true` and `shift = ...`, note
    that the role of `B * b` becomes computing `inv(A - shift I) * b`. So rather than
    passing the linear map ``A`` itself, pass a linear map `B` that has the action of
    shift-and-invert. The eigenvalue is transformed back to an eigenvalue of the actual
    matrix ``A``.

# Return values

**if `log` is `false`**
- `λ::Number` approximate eigenvalue computed as the Rayleigh quotient;
- `x::Vector` approximate eigenvector.

**if `log` is `true`**
- `λ::Number`: approximate eigenvalue computed as the Rayleigh quotient;
- `x::Vector`: approximate eigenvector;
- `history`: convergence history.

**ConvergenceHistory keys**

- `:tol` => `::Real`: stopping tolerance;
- `:resnom` => `::Vector`: residual norm at each iteration.

# Examples

```julia
using LinearMaps
σ = 1.0 + 1.3im
A = rand(ComplexF64, 50, 50)
F = lu(A - σ * I)
Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), 50, ismutating = true)
λ, x = powm(Fmap, inverse = true, shift = σ, tol = 1e-4, maxiter = 200)
```
"""
function powm!(B, x;
    tol = eps(real(eltype(B))) * size(B, 2) ^ 3,
    maxiter = size(B, 1),
    shift = zero(eltype(B)),
    inverse::Bool = false,
    log::Bool = false,
    verbose::Bool = false
)
    history = ConvergenceHistory(partial = !log)
    history[:tol] = tol
    reserve!(history, :resnorm, maxiter)
    verbose && @printf("=== powm ===\n%4s\t%7s\n", "iter", "resnorm")

    iterable = powm_iterable!(B, x, tol = tol, maxiter = maxiter)

    for (iteration, residual) = enumerate(iterable)
        nextiter!(history, mvps = 1)
        verbose && @printf("%3d\t%1.2e\n", iteration, residual)
    end

    setconv(history, converged(iterable))

    verbose && println()

    log && shrink!(history)

    λ = transform_eigenvalue(iterable.θ, inverse, shift)
    x = iterable.x

    log ? (λ, x, history) : (λ, x)
end

"""
    invpowm(B; shift = σ, kwargs...) -> λ, x, [history]

Find the approximate eigenpair `(λ, x)` of ``A`` near `shift`, where `B`
is a linear map that has the effect `B * v = inv(A - σI) * v`.

The method calls `powm!(B, x0; inverse = true, shift = σ)` with
`x0` a random, complex unit vector. See [`powm!`](@ref)

# Examples

```julia
using LinearMaps
σ = 1.0 + 1.3im
A = rand(ComplexF64, 50, 50)
F = lu(A - σ * I)
Fmap = LinearMap{ComplexF64}((y, x) -> ldiv!(y, F, x), 50, ismutating = true)
λ, x = invpowm(Fmap, shift = σ, tol = 1e-4, maxiter = 200)
```
"""
function invpowm(B; kwargs...)
    x0 = rand(Complex{real(eltype(B))}, size(B, 1))
    rmul!(x0, one(eltype(B)) / norm(x0))
    invpowm!(B, x0; kwargs...)
end

"""
    invpowm!(B, x0; shift = σ, kwargs...) -> λ, x, [history]

Find the approximate eigenpair `(λ, x)` of ``A`` near `shift`, where `B`
is a linear map that has the effect `B * v = inv(A - σI) * v`.

The method calls `powm!(B, x0; inverse = true, shift = σ)`. See [`powm!`](@ref).
"""
invpowm!(B, x0; kwargs...) = powm!(B, x0; inverse = true, kwargs...)


text:
# Logging functionality

export Logger, PrintLogger, ErrorLogger
export push_info!, push_iteration_info!
import Base.push!
using Printf

    """
    abstract type Logger ; end

The type represents a way to log information throughout the algorithms.
Error history and other properties can be stored in
the logger. The most common `Logger`s are `PrintLogger` and
`ErrorLogger`.

As a method developer you want to use [`push_info!`](@ref) and
[`push_iteration_info!`](@ref)

 See also: [`PrintLogger`](@ref) and [`ErrorLogger`](@ref).
"""
abstract type Logger ; end

   """
    function push_info!(logger, [level,] v; continues::Bool=false)

Pushes a string `v` to the logger. If continues=true, the next
`push_info!` (or `push_iteration_info!`) is connected to this,
e.g. line-feed will be omitted.
"""
function push_info!(logger::Logger,v::String; continues::Bool=false)
    push_info!(logger,1,v;continues=continues);
end
   """
    function push_iteration_info!(logger, [level,] iter; kwargs)

Pushes information about a specific iteration `iter`.
Standardized keywords are `λ`, `err` and `v`.

"""
function push_iteration_info!(logger::Logger,iter; kwargs...)
    push_iteration_info!(logger,1,iter;kwargs...);
end

   """
    struct PrintLogger <: Logger ;
    function PrintLogger(displaylevel)

When you use this logger, you will obtain printouts in stdout,
no other logging. The `displaylevel` parameter specifies
how much should be printed. The higher the value, the more
is printed. Zero means no printouts.
"""
struct PrintLogger <: Logger ;
    displaylevel::Int
end

function push_info!(logger::PrintLogger,level,v::String;continues::Bool=false)
    if (logger.displaylevel>=level)
        print(v);
        if (!continues)
            println();
        end
    end
end

function push_iteration_info!(logger::PrintLogger,level,iter;
                              err=Inf,λ=NaN,v=NaN,
                              continues::Bool=false)
    if (logger.displaylevel>=level)
        print("iter ",iter, " err:",err, " λ=",λ);
        if (!continues)
            println()
        end
    end
end

    """
    struct ErrorLogger <: Logger
    ErrorLogger(nof_eigvals=100,nof_iters=100,displaylevel=1)

Use this object if you want to save the error history in a method.
The `displaylevel` is interpreted as in `PrintLogger`.
The kwargs `nof_eigvals` and `nof_iters` is
used to preallocate memory used for saving.

When you use this logger, the error of `push_iteration_info!`-calls
will be stored in `logger.errs`.

"""
struct ErrorLogger{T} <: Logger  where {T};
    errs::Matrix{T}
    printlogger::Logger
end

function ErrorLogger(nof_eigvals=100,nof_iters=100,displaylevel=1)
    s=Matrix{Float64}(undef,nof_iters,nof_eigvals)
    s[:] .= NaN;
    printlogger=PrintLogger(displaylevel)
    return ErrorLogger{eltype(s)}(s,printlogger);
end

function push_iteration_info!(logger::ErrorLogger,level,iter;
                              err=Inf,λ=NaN,v=NaN,
                              continues::Bool=false)
    if (iter<=size(logger.errs,1))
        if (size(err,1)<=size(logger.errs,2))
            err_vec=err;
            # Make sure err_vec is a vector
            if (err_vec isa Number)
                err_vec=[err_vec]
            end
            logger.errs[iter,1:size(err,1)]=err_vec;
        else
            if (logger.printlogger.displaylevel>1)
                println("Warning: Insufficient space in logger matrix");
            end
        end

    end

    push_iteration_info!(logger.printlogger,level,iter;
                              err=err,λ=λ,v=v,
                              continues=continues)

end
function push_info!(logger::ErrorLogger,level::Int,
                    v::String;continues::Bool=false)
    push_info!(logger.printlogger,level,v,continues=continues)
end


text:
module TextWrapTest

using TextWrap
using Test

@testset "empty" begin

@test wrap("") == ""
@test wrap("  \n") == ""
@test wrap("\n  ") == ""
@test wrap("\t\n\t\n") == ""
@test wrap("", initial_indent=2) == "  "

end # testset

@testset "plain" begin

@test wrap("a b") == "a b"
@test wrap("a\n\nb\n\n") == "a b"
@test wrap("\n\na\n\nb\n\n") == "a b"
@test wrap(" a  ", initial_indent=3) == "   a"
@test wrap(" a  ", initial_indent=">  ") == ">  a"

global text = """
    Julia is a high-level, high-performance dynamic programming language
    for technical computing, with syntax that is familiar to users of
    other technical computing environments. It provides a sophisticated
    compiler, distributed parallel execution, numerical accuracy, and an
    extensive mathematical function library."""

@test wrap(text) == """
    Julia is a high-level, high-performance dynamic programming language
    for technical computing, with syntax that is familiar to users of
    other technical computing environments. It provides a sophisticated
    compiler, distributed parallel execution, numerical accuracy, and an
    extensive mathematical function library."""

@test wrap(text, width=30) == """
    Julia is a high-level, high-
    performance dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments. It
    provides a sophisticated
    compiler, distributed parallel
    execution, numerical accuracy,
    and an extensive mathematical
    function library."""

@test wrap(text, width=30, fix_sentence_endings=true, break_on_hyphens=false) == """
    Julia is a high-level,
    high-performance dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments.  It
    provides a sophisticated
    compiler, distributed parallel
    execution, numerical accuracy,
    and an extensive mathematical
    function library."""

@test wrap(text, width=30, initial_indent=2, subsequent_indent=0) == """
      Julia is a high-level, high-
    performance dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments. It
    provides a sophisticated
    compiler, distributed parallel
    execution, numerical accuracy,
    and an extensive mathematical
    function library."""

@test wrap(text, width=32, initial_indent=">   ", subsequent_indent="> ") == """
    >   Julia is a high-level, high-
    > performance dynamic
    > programming language for
    > technical computing, with
    > syntax that is familiar to
    > users of other technical
    > computing environments. It
    > provides a sophisticated
    > compiler, distributed parallel
    > execution, numerical accuracy,
    > and an extensive mathematical
    > function library."""

end # testset

@testset "whitespace" begin

tabtext = "aaaaaaa\tbbbbbbb\t\ncccccc\tddddd\teeee  \tfff\tgg\th\n"

@test wrap(tabtext, width=20, replace_whitespace=true,  expand_tabs=true) ==
    "aaaaaaa bbbbbbb\ncccccc  ddddd   eeee\nfff     gg      h"
@test wrap(tabtext, width=20, replace_whitespace=false, expand_tabs=true) ==
    "aaaaaaa bbbbbbb\ncccccc  ddddd   eeee\nfff     gg      h\n"
@test wrap(tabtext, width=20, replace_whitespace=true,  expand_tabs=false) ==
    "aaaaaaa bbbbbbb\ncccccc ddddd eeee\nfff gg h"
@test wrap(tabtext, width=20, replace_whitespace=false, expand_tabs=false) ==
    "aaaaaaa\tbbbbbbb\ncccccc\tddddd\teeee\nfff\tgg\th\n"

@test wrap("abc\ndefg h i j k", width=8, replace_whitespace=false) ==
           "abc\ndefg h i\nj k"
@test wrap("abc\ndefg h i j k", width=8, replace_whitespace=true) ==
           "abc defg\nh i j k"

@test wrap("\t.\n\t\n") == "."
@test wrap("\t..\n\t\n", fix_sentence_endings=true) == ".."
@test wrap("\tabc\n\t\n", width=2, expand_tabs=false) == "ab\nc"
@test wrap("\t \tabc\n\t\n", width=2) == "ab\nc"
@test wrap("\ta.\n\t\n", width=1, fix_sentence_endings=false) == "a\n."
@test wrap("\ta.\n\t\n", width=1, replace_whitespace=false) == "a\n.\n\n"
@test wrap("\ta.\tb\n", width=8, expand_tabs=true, replace_whitespace=true) == "a.\nb"
@test wrap("\ta.\tb\n", width=8, expand_tabs=true, replace_whitespace=false) == "a.\nb\n"
@test wrap("\ta.\tb\n", width=8, expand_tabs=false, replace_whitespace=false) == "a.\tb\n"
@test wrap("\ta.\tb\n", width=9, expand_tabs=false, replace_whitespace=false) == "a.\tb\n"
@test wrap("\ta.\tb\n", width=9, expand_tabs=true, replace_whitespace=false) == "a.      b\n"
@test wrap("\ta.\tb\n", width=9, expand_tabs=false, replace_whitespace=false, fix_sentence_endings=true) == "a.\tb\n"
@test wrap("\ta.\tb\n", width=9, expand_tabs=true, replace_whitespace=false, fix_sentence_endings=true) == "a.      b\n"

end # testset

@testset "long words" begin

longwordstext = """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word
    λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew word is the 25-letter-long (including vowels) וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=172) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word
    λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew word is the 25-letter-long (including vowels) וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=171) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτρ
    ιμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων. The longest Hebrew word is the 25-letter-
    long (including vowels) וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=80) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest
    English word that appears in a major dictionary. A 79 letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft,
    was named the longest published word in the German language by the 1996 Guinness
    Book of World Records, but longer words are possible. In his comedy
    Assemblywomen (c. 392 BC) Aristophanes coined the 171-letter word λοπαδοτεμαχοσε
    λαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφ
    οφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew word is the 25-letter-long (including vowels)
    וכשלאנציקלופדיותינו (u'chshelentsiklopedioténu), which means 'and when our
    encyclopedias will have....'"""

@test wrap(longwordstext, width=79) == """
    The 45-letter word pneumonoultramicroscopicsilicovolcanoconiosis is the longest
    English word that appears in a major dictionary. A 79 letter word, Donaudampfsc
    hiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft, was named
    the longest published word in the German language by the 1996 Guinness Book of
    World Records, but longer words are possible. In his comedy Assemblywomen (c.
    392 BC) Aristophanes coined the 171-letter word λοπαδοτεμαχοσελαχογαλεοκρανιολε
    ιψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλ
    εκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων. The longest
    Hebrew word is the 25-letter-long (including vowels) וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu), which means 'and when our encyclopedias will
    have....'"""

@test wrap(longwordstext, width=45) == """
    The 45-letter word
    pneumonoultramicroscopicsilicovolcanoconiosis
    is the longest English word that appears in a
    major dictionary. A 79 letter word, Donaudamp
    fschiffahrtselektrizitätenhauptbetriebswerkba
    uunterbeamtengesellschaft, was named the
    longest published word in the German language
    by the 1996 Guinness Book of World Records,
    but longer words are possible. In his comedy
    Assemblywomen (c. 392 BC) Aristophanes coined
    the 171-letter word λοπαδοτεμαχοσελαχογαλεοκρ
    ανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκ
    ατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυ
    ονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγ
    ανοπτερύγων. The longest Hebrew word is the
    25-letter-long (including vowels)
    וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu), which means 'and
    when our encyclopedias will have....'"""

@test wrap(longwordstext, width=44) == """
    The 45-letter word pneumonoultramicroscopics
    ilicovolcanoconiosis is the longest English
    word that appears in a major dictionary. A
    79 letter word, Donaudampfschiffahrtselektri
    zitätenhauptbetriebswerkbauunterbeamtengesel
    lschaft, was named the longest published
    word in the German language by the 1996
    Guinness Book of World Records, but longer
    words are possible. In his comedy
    Assemblywomen (c. 392 BC) Aristophanes
    coined the 171-letter word λοπαδοτεμαχοσελαχ
    ογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαρα
    βομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπερισ
    τεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιρ
    αιοβαφητραγανοπτερύγων. The longest Hebrew
    word is the 25-letter-long (including
    vowels) וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu), which means
    'and when our encyclopedias will have....'"""

@test wrap(longwordstext, width=19) == """
    The 45-letter word
    pneumonoultramicros
    copicsilicovolcanoc
    oniosis is the
    longest English
    word that appears
    in a major
    dictionary. A 79
    letter word, Donaud
    ampfschiffahrtselek
    trizitätenhauptbetr
    iebswerkbauunterbea
    mtengesellschaft,
    was named the
    longest published
    word in the German
    language by the
    1996 Guinness Book
    of World Records,
    but longer words
    are possible. In
    his comedy
    Assemblywomen (c.
    392 BC)
    Aristophanes coined
    the 171-letter word
    λοπαδοτεμαχοσελαχογ
    αλεοκρανιολειψανοδρ
    ιμυποτριμματοσιλφιο
    καραβομελιτοκατακεχ
    υμενοκιχλεπικοσσυφο
    φαττοπεριστεραλεκτρ
    υονοπτοκεφαλλιοκιγκ
    λοπελειολαγῳοσιραιο
    βαφητραγανοπτερύγων
    . The longest
    Hebrew word is the
    25-letter-long
    (including vowels)
    וכשלאנציקלופדיותינו
    (u'chshelentsiklope
    dioténu), which
    means 'and when our
    encyclopedias will
    have....'"""

@test wrap(longwordstext, width=18) == """
    The 45-letter word
    pneumonoultramicro
    scopicsilicovolcan
    oconiosis is the
    longest English
    word that appears
    in a major
    dictionary. A 79
    letter word, Donau
    dampfschiffahrtsel
    ektrizitätenhauptb
    etriebswerkbauunte
    rbeamtengesellscha
    ft, was named the
    longest published
    word in the German
    language by the
    1996 Guinness Book
    of World Records,
    but longer words
    are possible. In
    his comedy
    Assemblywomen (c.
    392 BC)
    Aristophanes
    coined the
    171-letter word λο
    παδοτεμαχοσελαχογα
    λεοκρανιολειψανοδρ
    ιμυποτριμματοσιλφι
    οκαραβομελιτοκατακ
    εχυμενοκιχλεπικοσσ
    υφοφαττοπεριστεραλ
    εκτρυονοπτοκεφαλλι
    οκιγκλοπελειολαγῳο
    σιραιοβαφητραγανοπ
    τερύγων. The
    longest Hebrew
    word is the
    25-letter-long
    (including vowels)
    וכשלאנציקלופדיותינ
    ו (u'chshelentsikl
    opedioténu), which
    means 'and when
    our encyclopedias
    will have....'"""

@test wrap(longwordstext, width=18, break_long_words=false) == """
    The 45-letter word
    pneumonoultramicroscopicsilicovolcanoconiosis
    is the longest
    English word that
    appears in a major
    dictionary. A 79
    letter word,
    Donaudampfschiffahrtselektrizitätenhauptbetriebswerkbauunterbeamtengesellschaft,
    was named the
    longest published
    word in the German
    language by the
    1996 Guinness Book
    of World Records,
    but longer words
    are possible. In
    his comedy
    Assemblywomen (c.
    392 BC)
    Aristophanes
    coined the
    171-letter word
    λοπαδοτεμαχοσελαχογαλεοκρανιολειψανοδριμυποτριμματοσιλφιοκαραβομελιτοκατακεχυμενοκιχλεπικοσσυφοφαττοπεριστεραλεκτρυονοπτοκεφαλλιοκιγκλοπελειολαγῳοσιραιοβαφητραγανοπτερύγων.
    The longest Hebrew
    word is the
    25-letter-long
    (including vowels)
    וכשלאנציקלופדיותינו
    (u'chshelentsiklopedioténu),
    which means 'and
    when our
    encyclopedias will
    have....'"""

end # testset

@testset "escape codes" begin

N = Base.text_colors[:normal]
U = Base.text_colors[:underline]
B = Base.text_colors[:bold]
R = Base.text_colors[:reverse]

r = N * Base.text_colors[:red]
b = N * U * Base.text_colors[:light_blue]
y = N * B * Base.text_colors[:light_yellow]
m = N * B * R * Base.text_colors[134]

etext = """
    Julia is a $(r)high-level$(N), $(b)$(B)high$(b)-performance$(N) dynamic programming language
    for technical computing, with syntax that is familiar to users of
    other technical computing environments. It provides a $(y)sophisticated$(N)
    compiler, distributed parallel execution, numerical accuracy, and an
    $(m)extensive mathematical function library$(N)."""

@test wrap(etext) == """
    Julia is a $(r)high-level$(N), $(b)$(B)high$(b)-performance$(N) dynamic programming language
    for technical computing, with syntax that is familiar to users of
    other technical computing environments. It provides a $(y)sophisticated$(N)
    compiler, distributed parallel execution, numerical accuracy, and an
    $(m)extensive mathematical function library$(N)."""

@test wrap(etext, width=30) == """
    Julia is a $(r)high-level$(N), $(b)$(B)high$(b)-
    performance$(N) dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments. It
    provides a $(y)sophisticated$(N)
    compiler, distributed parallel
    execution, numerical accuracy,
    and an $(m)extensive mathematical
    function library$(N)."""

@test wrap(etext, width=30, fix_sentence_endings=true, break_on_hyphens=false) == """
    Julia is a $(r)high-level$(N),
    $(b)$(B)high$(b)-performance$(N) dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments.  It
    provides a $(y)sophisticated$(N)
    compiler, distributed parallel
    execution, numerical accuracy,
    and an $(m)extensive mathematical
    function library$(N)."""

@test wrap(etext, width=30, initial_indent=2, subsequent_indent=0) == """
      Julia is a $(r)high-level$(N), $(b)$(B)high$(b)-
    performance$(N) dynamic
    programming language for
    technical computing, with
    syntax that is familiar to
    users of other technical
    computing environments. It
    provides a $(y)sophisticated$(N)
    compiler, distributed parallel
    execution, numerical accuracy,
    and an $(m)extensive mathematical
    function library$(N)."""

@test wrap(etext, width=32, initial_indent=">   ", subsequent_indent="> ") == """
    >   Julia is a $(r)high-level$(N), $(b)$(B)high$(b)-
    > performance$(N) dynamic
    > programming language for
    > technical computing, with
    > syntax that is familiar to
    > users of other technical
    > computing environments. It
    > provides a $(y)sophisticated$(N)
    > compiler, distributed parallel
    > execution, numerical accuracy,
    > and an $(m)extensive mathematical
    > function library$(N)."""

@test wrap(etext, width=35, recognize_escapes=false) == """
    Julia is a $(r)high-level$(N),
    $(b)$(B)high$(b)-
    performance$(N) dynamic programming
    language for technical computing,
    with syntax that is familiar to
    users of other technical computing
    environments. It provides a
    $(y)sophisticated$(N)
    compiler, distributed parallel
    execution, numerical accuracy, and
    an $(m)extensive
    mathematical function library$(N)."""

end # testset

@testset "print" begin

tmpf = tempname()
try
    open(tmpf, "w") do f
        print_wrapped(f)
    end
    @test read(tmpf, String) == ""

    open(tmpf, "w") do f
        print_wrapped(f, initial_indent=2, subsequent_indent=2)
    end
    @test read(tmpf, String) == "  "

    open(tmpf, "w") do f
        print_wrapped(f, "")
    end
    @test read(tmpf, String) == ""

    open(tmpf, "w") do f
        print_wrapped(f, "abc", "def")
    end
    @test read(tmpf, String) == "abcdef"

    open(tmpf, "w") do f
        println_wrapped(f, "abc", "def")
    end
    @test read(tmpf, String) == "abcdef\n"

    open(tmpf, "w") do f
        redirect_stdout(f) do
            print_wrapped("abc", "def")
        end
    end
    @test read(tmpf, String) == "abcdef"

    open(tmpf, "w") do f
        redirect_stdout(f) do
            println_wrapped("abc", "def")
        end
    end
    @test read(tmpf, String) == "abcdef\n"

    open(tmpf, "w") do f
        print_wrapped(f, text, width=30, fix_sentence_endings=true, break_on_hyphens=false)
    end
    @test read(tmpf, String) == """
        Julia is a high-level,
        high-performance dynamic
        programming language for
        technical computing, with
        syntax that is familiar to
        users of other technical
        computing environments.  It
        provides a sophisticated
        compiler, distributed parallel
        execution, numerical accuracy,
        and an extensive mathematical
        function library."""

    open(tmpf, "w") do f
        println_wrapped(f, text, width=30, fix_sentence_endings=true, break_on_hyphens=false)
    end
    @test read(tmpf, String) == """
        Julia is a high-level,
        high-performance dynamic
        programming language for
        technical computing, with
        syntax that is familiar to
        users of other technical
        computing environments.  It
        provides a sophisticated
        compiler, distributed parallel
        execution, numerical accuracy,
        and an extensive mathematical
        function library.
        """

    open(tmpf, "w") do f
        println_wrapped(f, width=30, fix_sentence_endings=true, break_on_hyphens=false)
    end
    @test read(tmpf, String) == "\n"

    open(tmpf, "w") do f
        bk_stdout = stdout
        try
            redirect_stdout(f)
            print_wrapped(text, width=30, fix_sentence_endings=true, break_on_hyphens=false)
        finally
            redirect_stdout(bk_stdout)
        end
    end
    @test read(tmpf, String) == """
        Julia is a high-level,
        high-performance dynamic
        programming language for
        technical computing, with
        syntax that is familiar to
        users of other technical
        computing environments.  It
        provides a sophisticated
        compiler, distributed parallel
        execution, numerical accuracy,
        and an extensive mathematical
        function library."""
finally
    isfile(tmpf) && rm(tmpf)
end

end # testset

@testset "argument checks" begin

@test_throws ArgumentError wrap("", width=0)
@test_throws ArgumentError wrap("", initial_indent=10, width=10)
@test_throws ArgumentError wrap("", subsequent_indent=10, width=10)
@test_throws ArgumentError wrap("", initial_indent="~~~~~~~~~~", width=10)
@test_throws ArgumentError wrap("", subsequent_indent="~~~~~~~~~~", width=10)

end # testset

end # module


text:
#
# Deleting
# --------

function deletefirst!(t::SegmentBTree{K, V, B}, first::K, last::K) where {K, V, B}
    return deletefirst!(t, Segment{K}(first, last))
end

function deletefirst!(t::SegmentBTree{K, V, B}, key::Tuple{K, K}) where {K, V, B}
    return deletefirst!(t, Segment{K}(key[1], key[2]))
end

function deletefirst!(t::SegmentBTree{K, V, B}, key::Segment{K}) where {K, V, B}
    result = _deletefirst!(t.root, key)
    if result.keyfound
        t.n -= 1
    end

    # if the root has only one child, promote the child
    if isa(t.root, InternalNode) && length(t.root) == 1
        t.root = t.root.children[1]
        t.root.parent = nothing
    end

    return t
end

# Indicate what steps are need to account for an updated child.
struct KeyFate
    value::UInt8
end

const KEYFATE_NONE         = KeyFate(0) # no changes
const KEYFATE_DELETE       = KeyFate(1) # delete the child node
const KEYFATE_DELETE_RIGHT = KeyFate(2) # delete the child's right sibling
const KEYFATE_UPDATE_LEFT  = KeyFate(3) # update the key separating the node
                                        # from its left sibling
const KEYFATE_UPDATE_RIGHT = KeyFate(4) # update the key separating the node
                                        # from its right sibling

# Information returned from a call to _delete
struct DeleteResult
    keyfound::Bool # true if the key was found

    # Indicate what steps are need to account for an updated child. One of:
    fate::KeyFate
end

# Delete key from the subtree if present. Return a DeleteResult object.
function _deletefirst!(t::InternalNode{K, V, B}, key::Segment{K}) where {K, V, B}
    i = findidx(t, key)
    j = i <= length(t) - 1 && key >= t.keys[i] ? i + 1 : i # child index

    ans = _deletefirst!(t.children[j], key)

    if !ans.keyfound
        return DeleteResult(false, KEYFATE_NONE)
    end

    if ans.fate == KEYFATE_UPDATE_LEFT
        t.maxend = nodemaxend(t)
        t.maxends[j] = t.children[j].maxend
        if j > 1
            t.maxends[j-1] = t.children[j-1].maxend
            t.keys[j-1] = minkey(t.children[j]);
            return DeleteResult(true, KEYFATE_NONE)
        else
            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    elseif ans.fate == KEYFATE_UPDATE_RIGHT
        t.maxend = nodemaxend(t)
        t.maxends[j] = t.children[j].maxend
        if j < length(t)
            t.maxends[j+1] = t.children[j+1].maxend
            t.keys[j] = minkey(t.children[j+1])
        end

        if j > 1
            t.keys[j - 1] = minkey(t.children[j]);
            return DeleteResult(true, KEYFATE_NONE)
        else
            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    elseif ans.fate == KEYFATE_DELETE || ans.fate == KEYFATE_DELETE_RIGHT
        deleteidx = ans.fate == KEYFATE_DELETE ? j : j + 1

        # deleteidx == 1 would only happen if we had only one child, but if
        # that were the case it would have been promoted.
        @assert deleteidx > 1

        splice!(t.keys, deleteidx - 1)
        splice!(t.children, deleteidx)
        splice!(t.maxends, deleteidx)
        t.maxend = nodemaxend(t)
        if ans.fate == KEYFATE_DELETE
            t.maxends[j-1] = t.children[j-1].maxend
        elseif ans.fate == KEYFATE_DELETE_RIGHT
            t.maxends[j] = t.children[j].maxend
        end

        # not underfull
        minsize = div(B, 2)
        if length(t) >= minsize
            return DeleteResult(true, deleteidx == 2 ? KEYFATE_UPDATE_LEFT : KEYFATE_NONE)
        end

        # borrow right
        if t.right !== nothing
            right = notnothing(t.right)
            if sameparent(right, t) && length(right) > minsize
                splice!(right.keys, 1)
                push!(t.children, splice!(right.children, 1))
                push!(t.maxends, splice!(right.maxends, 1))
                push!(t.keys, minkey(t.children[end]))
                t.maxend = max(t.maxend, nodemaxend(t.children[end]))
                right.maxend = nodemaxend(right)
                t.children[end].parent = t

                return DeleteResult(true, KEYFATE_UPDATE_RIGHT)
            end
        end

        # borrow left
        if t.left !== nothing
            left = notnothing(t.left)
            if sameparent(left, t) && length(left) > minsize
                insert!(t.children, 1, pop!(left.children))
                insert!(t.maxends, 1, pop!(left.maxends))
                pop!(left.keys)
                insert!(t.keys, 1, minkey(t.children[2]))
                t.maxend = max(t.maxend, nodemaxend(t.children[1]))
                left.maxend = nodemaxend(left)
                t.children[1].parent = t

                return DeleteResult(true, KEYFATE_UPDATE_LEFT)
            end
        end

        # merge with left
        if t.left !== nothing
            left = notnothing(t.left)
            if sameparent(left, t)
                merge!(left, t)
                left.right = t.right
                if t.right !== nothing
                    t.right = notnothing(t.right)
                    t.right.left = left
                end

                return DeleteResult(true, KEYFATE_DELETE)
            end
        end

        # merge with right
        if t.right !== nothing
            right = notnothing(t.right)
            if sameparent(right, t)
                merge!(t, right)
                if right.right !== nothing
                    right.right = notnothing(right.right)
                    right.right.left = t
                end
                t.right = right.right

                return DeleteResult(true, KEYFATE_DELETE_RIGHT)
            end
        end

        # Allow the root node to be underfull
        return DeleteResult(true, KEYFATE_NONE)
    else # KEYFATE_NONE
        t.maxend = nodemaxend(t)
        t.maxends[j] = t.children[j].maxend
        return DeleteResult(true, KEYFATE_NONE)
    end
end

function _deletefirst!(t::LeafNode{K, V, B}, key::Segment{K}) where {K, V, B}
    i = findidx(t, key)

    # do nothing if the key isn't present
    if i < 1 || i > length(t) ||
        first(t.keys[i]) != key.first || last(t.keys[i]) != key.last
        return DeleteResult(false, KEYFATE_NONE)
    end

    splice!(t, i)

    # This is the root node. Allow it to be empty.
    if isempty(t)
        return DeleteResult(true, KEYFATE_NONE)
    end

    minsize = div(B, 2)

    t.maxend = nodemaxend(t)

    # not underfull
    if length(t) >= minsize
        return DeleteResult(true, i == 1 ? KEYFATE_UPDATE_LEFT : KEYFATE_NONE)
    end

    # borrow right
    if t.right !== nothing
        right = notnothing(t.right)
        if sameparent(right, t) && length(right) > minsize
            push!(t, splice!(right, 1))
            t.maxend = max(t.maxend, last(t.keys[t.count]))
            right.maxend = nodemaxend(right)

            return DeleteResult(true, KEYFATE_UPDATE_RIGHT)
        end
    end

    # borrow left
    if t.left !== nothing
        left = notnothing(t.left)
        if sameparent(left, t) && length(left) > minsize
            insert!(t, 1, pop!(left))
            t.maxend = max(t.maxend, last(t.keys[1]))
            left.maxend = nodemaxend(left)

            return DeleteResult(true, KEYFATE_UPDATE_LEFT)
        end
    end

    # merge with left
    if t.left !== nothing
        left = notnothing(t.left)
        if sameparent(left, t)
            merge!(left, t)
            left.right = t.right
            if t.right !== nothing
                t.right = notnothing(t.right)
                t.right.left = left
            end
            return DeleteResult(true, KEYFATE_DELETE)
        end
    end

    # merge with right
    if t.right !== nothing
        right = notnothing(t.right)
        if sameparent(right, t)
            merge!(t, right)
            if right.right !== nothing
                right.right = notnothing(right.right)
                right.right.left = t
            end
            t.right = right.right
            return DeleteResult(true, KEYFATE_DELETE_RIGHT)
        end
    end

    # This must be the root node. Allow it to be underfull.
    return DeleteResult(true, KEYFATE_NONE)
end


text:
using TypeStability
using Compat
using Compat.Test

@testset "TypeStability.jl" begin

    include("StabilityAnalysisTests.jl")
    include("InlineCheckerTests.jl")
    include("Utils.jl")
end


text:
_checkweights(matrix) = (size(matrix, 1) == size(matrix, 2)) ? size(matrix, 1) : error("Connectivity (weight) matrix must be a square.")

@userplot RasterPlot

@recipe function f(h::RasterPlot)
    seriestype := :scatter
    if length(h.args) == 1 && isa(h.args[1], Dict)
        yguide := "Neuron ID"
        for (i, x) in pairs(h.args[1])
            @series begin
                if !haskey(plotattributes, :label) || isempty(plotattributes[:label])
                    label := ""
                else
                    label := plotattributes[:label][i]
                end
                x, fill(float(i), length(x))
            end
        end
    else
        for (i, x) in enumerate(h.args)
            yticks := :none
            @series begin
                if !haskey(plotattributes, :label) || isempty(plotattributes[:label])
                    label := ""
                else
                    label := plotattributes[:label][i]
                end
                float(x), fill(float(i), length(x))
            end
        end
    end
end


text:
# load the include -
include("Include.jl")

# load the training data -
path_to_training_data = joinpath(_PATH_TO_DATA, "Training-Synthetic-Thrombin-TF-1K.csv")
training_df = CSV.read(path_to_training_data, DataFrame)

# size of training set -
(R,C) = size(training_df)

# build the model structure -
path_to_model_file = joinpath(_PATH_TO_MODEL, "Coagulation.net")
model_buffer = read_model_file(path_to_model_file)

# build the default model structure -
model = build_default_model_dictionary(model_buffer)

# main simulation loop -
SF = 1e9
for i ∈ 1:10

    # build new model -
    dd = deepcopy(model)

    # setup static -
    sfa = dd["static_factors_array"]
    sfa[1] = training_df[i,:TFPI]       # 1 TFPI
    sfa[2] = training_df[i,:AT]         # 2 AT
    sfa[3] = (5e-12) * SF               # 3 TF
    sfa[6] = 0.005                      # 6 TRAUMA

    # grab the multiplier from the data -
    ℳ = dd["number_of_dynamic_states"]
    xₒ = zeros(ℳ)
    xₒ[1] = training_df[i, :II]         # 1 FII
    xₒ[2] = training_df[i, :VII]        # 2 FVII
    xₒ[3] = training_df[i, :V]          # 3 FV
    xₒ[4] = training_df[i, :X]          # 4 FX
    xₒ[5] = training_df[i, :VIII]       # 5 FVIII
    xₒ[6] = training_df[i, :IX]         # 6 FIX
    xₒ[7] = training_df[i, :XI]         # 7 FXI
    xₒ[8] = training_df[i, :XII]        # 8 FXII
    xₒ[9] = (1e-14)*SF                  # 9 FIIa
    dd["initial_condition_vector"] = xₒ

    # update α -
    α = dd["α"]
    α[1] = 0.061
    α[9] = 0.70

    # setup -
    G = dd["G"]
    idx = indexin(dd, "FVIIa")
    G[idx, 4] = 0.1

    # what is the index of TRAUMA?
    idx = indexin(dd, "AT")
    G[idx, 9] = 0.045

    # what is the index of TFPI?
    idx = indexin(dd, "TFPI")
    G[idx, 1] = -0.65

    # run the model -
    global (T,U) = evaluate(dd)
    X = hcat(U...)
    data = [T transpose(X)]

    # dump -
    path_to_sim_data = joinpath(_PATH_TO_TMP, "SIM-TF-NO-TM-SYN1K-$(i).dat")
    CSV.write(path_to_sim_data, Tables.table(data))
end

max_stars_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_issues_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation

max_forks_repo_name: varnerlab/BST-S-Sytem-Model-Coagulation


text:
module Lists

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export list, item, item_section, itemsection, item_label, itemlabel
export List, Item

register_normal_element("q__list", context = @__MODULE__)
register_normal_element("q__item", context = @__MODULE__)
register_normal_element("q__item__section", context = @__MODULE__)
register_normal_element("q__item__label", context = @__MODULE__)

# OOP API
mutable struct List
  args
  kwargs
end

List(args...; kwargs...) = List(args, kwargs)

Base.string(l::List) = list(l.args...; l.kwargs...)

mutable struct Item
  args
  kwargs
end

Item(args...; kwargs...) = Item(args, kwargs)

function Item(i::Item, args...; kwargs...)
  ix = Item()
  ix.args = (i.args..., args...)
  ix.kwargs = (i.kwargs..., kwargs...)

  ix
end

Base.string(i::Item) = item(i.args...; i.kwargs...)

# Functional API

function list(args...; kwargs...)
  q__list(args...; kw(kwargs)...)
end

# server side iteration
function list(data::Vector{T}, i::Item, args...; kwargs...) where {T}
  content = []
  for d in data
    push!(content, Item(i, [d]) |> string)
  end

  list(join(content), args...; kwargs...)
end

# client side iteration
function list(binding::Symbol, i::Item, args...; kwargs...)
  list(
    Item(i, ["{{__b__}}"], @recur("__b__ in $binding"), key! = "__b__")
  )
end

# list items

function item(args...; kwargs...)
  q__item(args...; kw(kwargs)...)
end

function item_section(args...; kwargs...)
  q__item__section(args...; kw(kwargs)...)
end
const itemsection = item_section

function item_label(args...; kwargs...)
  q__item__label(args...; kw(kwargs)...)
end
const itemlabel = item_label

end


text:
module Iterative

using CUDA
using KernelAbstractions
using IterativeSolvers
using Krylov
using LinearAlgebra
using SparseArrays
using TimerOutputs

using ..ExaPF: Precondition

export bicgstab, list_solvers

cuzeros = CUDA.zeros

"""
bicgstab according to

Van der Vorst, Henk A.
"Bi-CGSTAB: A fast and smoothly converging variant of Bi-CG for the solution
of nonsymmetric linear systems."
SIAM Journal on scientific and Statistical Computing 13, no. 2 (1992): 631-644.
"""
function bicgstab(A, b, P, xi, to::TimerOutput;
                  tol=1e-6, maxiter=size(A, 1), verbose=false)
    # parameters
    n    = size(b, 1)
    x0   = similar(b)
    x0 = P * b
    r0   = b - A * x0
    br0  = copy(r0)
    rho0 = 1.0
    alpha = 1.0
    omega0 = 1.0
    if A isa SparseArrays.SparseMatrixCSC
        v0 = zeros(Float64, n)
        p0 = zeros(Float64, n)
    else
        v0 = cuzeros(Float64, n)
        p0 = cuzeros(Float64, n)
    end

    ri     = copy(r0)
    rhoi   = copy(rho0)
    omegai = copy(omega0)
    vi = copy(v0)
    pi = copy(p0)
    xi .= x0

    y = similar(pi)
    s = similar(pi)
    z = similar(pi)
    t1 = similar(pi)
    t2 = similar(pi)
    pi1 = similar(pi)
    vi1 = similar(pi)
    xi1 = similar(pi)
    t = similar(pi)

    go = true
    iter = 1
    @timeit to "While loop" begin
        while go
            @timeit to "First stage" begin
                rhoi1 = dot(br0, ri) ; beta = (rhoi1/rhoi) * (alpha / omegai)
                pi1 .= ri .+ beta .* (pi .- omegai .* vi)
                y .= P * pi1
                vi1 .= A * y
                alpha = rhoi1 / dot(br0, vi1)
                s .= ri .- (alpha * vi1)
                z .= P * s
                t .= A * z
                t1 .= P * t
                t2 .= P * s
            end
            @timeit to "Main stage" begin
                omegai1 = dot(t1, t2) / dot(t1, t1)
                xi1 .= xi .+ alpha .* y .+ omegai1 .* z
            end

            @timeit to "End stage" begin
                anorm = norm((A * xi1) .- b)

                if verbose
                    println("\tIteration: ", iter)
                    println("\tAbsolute norm: ", anorm)
                end

                if anorm < tol
                    go = false
                    verbose && println("Tolerance reached at iteration $iter")
                end

                if maxiter == iter
                    @show iter
                    go = false
                    verbose && println("Not converged")
                end

                ri     .= s .- omegai1 .* t
                rhoi   = rhoi1
                pi     .= pi1
                vi     .= vi1
                omegai = omegai1
                xi     .= xi1
                iter   += 1
            end
        end
    end

    return xi, iter
end

list_solvers(::CPU) = ["bicgstab_ref", "bicgstab", "gmres", "dqgmres", "default"]

function ldiv!(
    dx::AbstractVector,
    J::SparseArrays.SparseMatrixCSC,
    F::AbstractVector,
    solver::String,
    preconditioner,
    timer::TimerOutputs.TimerOutput=nothing,
)
    if solver == "bicgstab_ref"
        @timeit timer "Preconditioner" P = Precondition.update(J, preconditioner, timer)
        @timeit timer "CPU-BICGSTAB" (dx[:], history) = IterativeSolvers.bicgstabl(P*J, P*F, log=true)
        n_iters = history.iters
    elseif solver == "bicgstab"
        @timeit timer "Preconditioner" P = Precondition.update(J, preconditioner, timer)
        @timeit timer "GPU-BICGSTAB" dx[:], n_iters = bicgstab(J, F, P, dx, timer, maxiter=10000)
    elseif solver == "gmres"
        @timeit timer "Preconditioner" P = Precondition.update(J, preconditioner, timer)
        @timeit timer "CPU-GMRES" (dx[:], history) = IterativeSolvers.gmres(P*J, P*F, restart=4, log=true)
        n_iters = history.iters
    elseif solver == "dqgmres"
        @timeit timer "Preconditioner" P = Precondition.update(J, preconditioner, timer)
        @timeit timer "GPU-DQGMRES" (dx[:], status) = Krylov.dqgmres(J, F, M=P, memory=4)
        n_iters = length(status.residuals)
    else
        @timeit timer "CPU-Default sparse solver" dx .= J\F
        n_iters = 0
    end
    return n_iters
end

list_solvers(::CUDADevice) = ["bicgstab", "dqgmres", "default"]

function ldiv!(
    dx::CuVector,
    J::CUDA.CUSPARSE.CuSparseMatrixCSR,
    F::CuVector,
    solver::String,
    preconditioner,
    timer=nothing,
)
    if solver == "bicgstab"
        @timeit timer "Preconditioner" P = Precondition.update(J, preconditioner, timer)
        @timeit timer "GPU-BICGSTAB" dx[:], n_iters = bicgstab(J, F, P, dx, timer, maxiter=10000)
    elseif solver == "dqgmres"
        @timeit timer "Preconditioner" P = Precondition.update(J, preconditioner, timer)
        @timeit timer "GPU-DQGMRES" (dx[:], status) = Krylov.dqgmres(J, F, M=P, memory=4)
        n_iters = length(status.residuals)
    else
        lintol = 1e-8
        @timeit timer "Sparse CUSOLVER" dx  = CUSOLVER.csrlsvqr!(J,F,dx,lintol,one(Cint),'O')
        n_iters = 0
    end
    return n_iters
end

end


text:
using Pkg
Pkg.activate(".")

if Threads.nthreads() < 2
    error("Multiple threads NOT found! Makes no sense to test this without them... Start as, for example: JULIA_NUM_THREADS=4 julia multithreaded_optimization.jl")
    exit(-1)
end

using BlackBoxOptim

# Functions to optimize. Should be thread-safe.
function rosenbrock(x)
    sleep(0.05) # So that there is some benefit in the thread switching...
    sum(i -> 100*abs2(x[i+1] - x[i]^2) + abs2(x[i] - 1), Base.OneTo(length(x)-1))
end

function rastrigin(x)
    D = length(x)
    10 * D + sum(abs2, x) - 10 * sum(xx -> cos(2π * xx), x)
end

# For multi-objective opt:
optfun(x) = (rosenbrock(x), rastrigin(x))

# Let's ensure they are really thread-safe:
for (fn, N) in [(rosenbrock, 100), (rastrigin, 1000000)]
    points = rand(N, 2)
    @time single_threaded = map(i -> fn(points[i, :]), 1:N)
    multi_threaded = zeros(Float64, N)
    @time Threads.@threads for i in 1:N
        multi_threaded[i] = fn(points[i, :])
    end
    @assert sum((single_threaded .- multi_threaded).^2) == 0
end

# Short runs to ensure things have compiled:
bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
    MaxTime = 0.1, TraceMode = :silent, Method = :dxnes, lambda = 10);
bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 2,
    MaxTime = 0.1, TraceMode = :silent, NThreads=Threads.nthreads()-1, Method = :dxnes, lambda = 10);
bboptimize(optfun; SearchRange = (-5.0, 5.0), NumDimensions = 2,
    MaxTime = 0.5, TraceMode = :silent, NThreads=Threads.nthreads()-1,
    Method = :borg_moea, ϵ=0.05,
    FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true));

# Same time given to single- and multi-threaded optimization runs:
MaxTime = 30.0

# Now run the single-threaded optimization:
@time res_single = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
    MaxTime = MaxTime, TraceMode = :silent, Method = :dxnes, lambda = 50)

# Now run the multi-threaded optimization:
@time res_multi = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
    MaxTime = MaxTime, TraceMode = :silent, NThreads=Threads.nthreads()-1,
    Method = :dxnes, lambda = 50)

println("Fitness (single-threaded): ", round(best_fitness(res_single), digits=4))
println("Fitness (multi-threaded):  ", round(best_fitness(res_multi), digits=4))
#@assert BlackBoxOptim.f_calls(res_single) < BlackBoxOptim.f_calls(res_multi)

@time res_single = bboptimize(optfun; SearchRange = (-5.0, 5.0), NumDimensions = 100,
    MaxTime = MaxTime, TraceMode = :silent,
    Method = :borg_moea, ϵ=0.05,
    FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true));

@time res_multi = bboptimize(optfun; SearchRange = (-5.0, 5.0), NumDimensions = 100,
    MaxTime = MaxTime, TraceMode = :silent, NThreads=Threads.nthreads()-1,
    Method = :borg_moea, ϵ=0.05,
    FitnessScheme=ParetoFitnessScheme{2}(is_minimizing=true));

# The whole idea with multi-threading is that we can evaluate more function calls in same time so:
@assert BlackBoxOptim.f_calls(res_single) < BlackBoxOptim.f_calls(res_multi)

# Now do repeated runs and compare if there is a real difference in fitness
using HypothesisTests
using Statistics

for method in [:adaptive_de_rand_1_bin_radiuslimited, :dxnes]
    fsingle, fmulti = Float64[], Float64[]
    for i in 1:10
        println("Single-threaded run $i")
        res_single = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
            Method = method, MaxTime = MaxTime, TraceMode = :silent, lambda = 50)
        push!(fsingle, best_fitness(res_single))
        println("Fitness (single-threaded): ", round(last(fsingle), digits=2))

        println("Multi-threaded run $i")
        res_multi = bboptimize(rosenbrock; SearchRange = (-5.0, 5.0), NumDimensions = 100,
            Method = method, MaxTime = MaxTime, TraceMode = :silent, lambda = 50,
            NThreads=Threads.nthreads()-1)
        push!(fmulti, best_fitness(res_multi))
        println("Fitness (multi-threaded): ", round(last(fmulti), digits=2))
    end

    pval = pvalue(MannWhitneyUTest(fsingle, fmulti))

    # Seems to give consistently better fitness in single-threaded runs on my 4-core
    # (as well as 2-core (but more expected there)) laptops for adaptive_de_rand_1_bin_radiuslimited.
    # Tried for MaxTime 5, 10, and 60 seconds with same results...
    println("For method: ", method)
    println("Mean Fitness (single-threaded): ", round(mean(fsingle), digits=2))
    println("Mean Fitness (multi-threaded):  ", round(mean(fmulti), digits=2))
    println("Significant difference: ", (pval < 0.01 ? "Yes" : "No"), " ($(round(pval, digits=5)))")
end


text:
b = SpinBasis(1//2)
up = spinup(b)
down = spindown(b)
# id()= identityoperator(b)

@test id() == identityoperator(b)
@test id(3) == identityoperator(b) ⊗ identityoperator(b) ⊗ identityoperator(b)
@test sigmax() == sigmax(SpinBasis(1//2))
@test X() == sigmax()
@test sigmay() == sigmay(SpinBasis(1//2))
@test Y() == sigmay()
@test sigmaz() == sigmaz(SpinBasis(1//2))
@test Z() == sigmaz()
@test sigmam() == sigmam(SpinBasis(1//2))
@test sigmap() == sigmap(SpinBasis(1//2))
@test hadamard() == DenseOperator(b, 1/sqrt(2)*[1 1; 1 -1])
@test H() == hadamard()
@test hadamard(2) == DenseOperator(b, 1/sqrt(2)*[1 1; 1 -1]) ⊗ DenseOperator(b, 1/sqrt(2)*[1 1; 1 -1])
@test sgate() == sparse(DenseOperator(b, [1 0; 0 im]))
@test S() == sgate()
@test sdggate() == sparse(DenseOperator(b, [1 0; 0 -im]))
@test Sdg() == sdggate()
@test tgate() == sparse(DenseOperator(b, [1 0; 0 (1+im)/sqrt(2)]))
@test T() == tgate()
@test tdggate() == sparse(DenseOperator(b, [1 0; 0 (1-im)/sqrt(2)]))
@test Tdg() == tdggate()
@test phaseshift(π/3) == sparse(DenseOperator(b, [1 0; 0 exp(im * π / 3)]))
@test U1(π/3) ≈ phaseshift(π/3)
@test phaseshift(π/3) ≈ singleunitary(0,0,π/3)
@test U2(π/3, π/6) ≈ singleunitary(π/2,π/3, π/6)
@test U3(π/2,π/3, π/6) ≈ singleunitary(π/2,π/3, π/6)
@test controlsgate(2, [1], 2, sigmax()) == cnot()
@test controlsgate(3, [2, 3], 1, sigmax()) ==  identityoperator(b) ⊗ (identityoperator(b^2) - sparse_spinalldown_dm(2)) + sigmax() ⊗ sparse_spinalldown_dm(2)
@test controlsgate(4, [3, 4], 2, sigmax()) == identityoperator(b) ⊗ controlsgate(3, [2, 3], 1, sigmax())
@test controlsgate(7, [1,2,5], 3, sigmax()) == controlsgate(7, [5,1,2], 3, sigmax())
@test cnot() == dm(up) ⊗ id()+ dm(down) ⊗ sigmax()
@test cnot() == CNOT()
@test cnot(2, 1, 2) == dm(up) ⊗ id()+ dm(down) ⊗ sigmax()
@test cnot(2, 2, 1) == id()⊗ dm(up) + sigmax() ⊗ dm(down)
@test cnot(3, 1, 2) == dm(up) ⊗ id()⊗ id()+ dm(down) ⊗ sigmax() ⊗ id()
@test cnot(3, 1, 3) == dm(up) ⊗ id()⊗ id()+ dm(down) ⊗ id()⊗ sigmax()
@test cnot(3, 2, 3) == identityoperator(b) ⊗ (dm(up) ⊗ id()+ dm(down) ⊗ sigmax())
@test cnot(3, 3, 2) == identityoperator(b) ⊗ (id()⊗ dm(up) + sigmax() ⊗ dm(down))
@test toffoli() == (identityoperator(b^2) - sparse_spinalldown_dm(2)) ⊗ identityoperator(b) + sparse_spinalldown_dm(2) ⊗ sigmax()
@test try toffoli(1, [1,2], 3) catch true end
@test try toffoli(3, [1,3], 3) catch true end
@test CCX() == toffoli()
@test swap() == sparse(dm(qubit("00")) + outer(qubit("01"), qubit("10")) + outer(qubit("10"), qubit("01")) + dm(qubit("11")))
@test SWAP() == swap()
@test swap(3, [1,2]) == swap() ⊗ id()
@test swap(3, [3,2]) == id() ⊗ swap()
@test swap(3, [1,3]) == cnot(3, 1, 3) * cnot(3, 3, 1) * cnot(3, 1, 3)
try swap(1, [1,2]) catch true end
try swap(2, [1,1]) catch true end
try swap(3, [1,2,3]) catch true end


text:
isapproxzero(a) = isapprox(a, 0.0; atol = 1e-8)

function findMaxAbsInColumn(A::Matrix, k)
    rows, _ = size(A)
    maxPosition = k
    for i in k+1:rows
        if abs(A[i, k]) > abs(A[maxPosition, k])
            maxPosition = i
        end
    end

    return maxPosition
end

function column_to_zeroes!(A, k, maxPosition)
    if maxPosition != k
        A[k, :], A[maxPosition, :] = A[maxPosition, :], A[k, :]
    end

    if isapproxzero(A[k,k])
        for i in k:lastindex(A,1)
            A[i, k] = 0
            return
        end
    end

    for i in k+1:(lastindex(A, 1))
        if isapproxzero(A[i,k])
            A[i,k] = 0
            continue
        end
        t = -A[i,k]/A[k,k]
        A[i, i:end] += t *(A[k, i:end])
        A[i,k] = 0
    end
end

function convert!(A)
    rows, columns = size(A)
    for k in 1:rows
        imax = findMaxAbsInColumn(A, k)
        column_to_zeroes!(A, k, imax)
    end
end


text:
using Revise
using LightField
using Plots; gr()
cd("/home/wikphi/LightField.jl/")
lf = setup("configexample.toml")
@time result = propagate(lf)

normal = result[:,:,:,:,1] ./ maximum(result[:,:,:,:,1])
using FFTW
normal[normal .< (0.005*maximum(normal))] .= 0
heatmap(normal[:,:,13,27])