text:
function ueberpruefeMatrixVektor2x2(A,m)

Akorrekt = [1 0 1 0;
            0 1 0 1;
            1 1 0 0;
            0 0 1 1];

mkorrekt = [5 8  6 7];

# homogene Lösungen berechnen
homogenkorrekt = nullspace(Akorrekt);

# inhomogene Lösung berechnen
inhomogenkorrekt = pinv(Akorrekt)*mkorrekt';

homogen = nullspace(A);
inhomogen = pinv(A)*m';

u1 = homogen[:,1];
# u2 = homogen(:,2);
# u3 = homogen(:,3);
# u4 = homogen(:,4);
p = inhomogen;

Vergleich = Akorrekt'*mkorrekt';
VergleichSuS = A'*m';

if norm(Vergleich) == norm(VergleichSuS)
    println("Eure Matrix und euer Vektor sind korrekt!")
    println("Die Lösung eures Gleichungssystems lautet wie folgt, wobei r eine beliebige relle Zahl ist.")
    q = round.(p[:,1];digits=2);
    v1 = round.(u1[:,1];digits=2);
    print("Für f1 lautet die Lösung: ", q[1], " + r * ", v1[1],"\n")
    print("Für f2 lautet die Lösung: ", q[2], " + r * ", v1[2],"\n")
    print("Für f3 lautet die Lösung: ", q[3], " + r * ", v1[3],"\n")
    print("Für f4 lautet die Lösung: ", q[4], " + r * ", v1[4],"\n")
    #show(stdout,"text/plain",round.(q;digits=2))
    #show(stdout,"text/plain",round.(v1;digits=2))
        #println("p = $p")
    #println("u1 = $u1")
elseif norm(m) == norm([15 12 14 13])
    println("Formt euer Gleichungssystem um, so dass alle Variablen rechts und Zahlen links stehen. Stellt dann die Matrix-Vektorschreibweise auf.")
else
    println("Schaut euch noch einmal den Verlauf der Strahlen an und überdenkt eure Matrix.")
end

end


text:
# AdaM: https://arxiv.org/pdf/1412.6980v8.pdf
function adamUpdate(params, grads, adamParams, encLimit=3, decLimit=2, b1=.95, b2=.999, e=1e-8)
  adamParams["t"] += 1

  for paramType in ["w_encoder", "b_encoder", "w_decoder", "b_decoder"]
    if contains(paramType,"encoder") idx_limit = encLimit else idx_limit = decLimit end
    for param_idx in 1:idx_limit
      # mean term
      adamParams["m"][paramType][param_idx] = b1*adamParams["m"][paramType][param_idx] + (1-b1)*grads[paramType][param_idx]
      m_hat = adamParams["m"][paramType][param_idx] ./ (1-b1.^adamParams["t"])

      # variance term
      adamParams["v"][paramType][param_idx] = b2*adamParams["v"][paramType][param_idx] + (1-b2)*grads[paramType][param_idx].^2
      v_hat = adamParams["v"][paramType][param_idx] ./ (1-b2.^adamParams["t"])

      # update model param
      params[paramType][param_idx] -= adamParams["lr"] * m_hat./(sqrt(v_hat) + e)
    end
  end

  return params
end

ext: jl


text:
@inline function add_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yhi, ylo = y
    yhi, ylo = two_sumof3(x, yhi, ylo)
    return yhi, ylo
end

@inline function sub_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yhi, ylo = y
    yhi, ylo = two_diffof3(x, yhi, ylo)
    return yhi, ylo
end

@inline function mul_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    hi, lo = y
    hihi, hilo = two_prod(x, hi)
    lohi, lolo = two_prod(x, lo)
    hi, lo = two_sumof4(hihi, hilo, lohi, lolo)
    return hi, lo
end

@inline function dvi_fpdd_dd(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yinv = inv_dd_dd(y)
    zhi, zlo = mul_fpdd_dd(x, yinv)
    return zhi, zlo
end

@inline function dvi_fpdd_dd_fast(x::T, y::Tuple{T,T}) where {T<:IEEEFloat}
    yinv = inv_dd_dd_fast(y)
    zhi, zlo = mul_fpdd_dd(x, yinv)
    return zhi, zlo
end


text:
"""
Functions for area of various geometric shapes

Contributed By:- [Ashwani Rathee](https://github.com/ashwani-rathee)
"""

"""
    surfarea_cube(side)

Finds surface area of a cube

# Example

```julia
surfarea_cube(1)  # returns 6
surfarea_cube(3)  # returns 54
surfarea_cube(-1) # returns DomainError
```
"""
function surfarea_cube(side)
    if side < 0
        throw(DomainError(side, "surfarea_cube() only works for non-negative values"))
    end
    return 6(side^2)
end

"""
    surfarea_sphere(side)

Finds surface area of a sphere

# Example

```julia
surfarea_sphere(5)  # returns 314.1592653589793
surfarea_sphere(1)  # returns 12.566370614359172
surfarea_sphere(-1) # returns DomainError
```
"""
function surfarea_sphere(radius)
    if radius < 0
        throw(DomainError(radius, "surfarea_sphere() only works for non-negative values"))
    end
    return 4pi*(radius^2)
end

"""
    area_rectangle(length, width)

Finds area of the rectangle

# Example

```julia
area_rectangle(10,20) # returns 200
area_rectangle(-1,-2) # returns DomainError
area_rectangle(1,-2)  # returns DomainError
area_rectangle(-1,2)  # returns DomainError
```
"""
function area_rectangle(length, width)
    if length < 0 || width < 0
        throw(DomainError("area_rectangle() only works for non-negative values"))
    end
    return length * width
end

"""
    area_square(side)

Finds area of the area_square

# Example

```julia
area_square(10) # returns 100
area_square(-1) # returns DomainError
```
"""
function area_square(side)
    if side < 0
        throw(DomainError("area_square() only works for non-negative values"))
    end
    return side^2
end

"""
    area_triangle(base, height)

Finds area of the right angled triangle with base height

# Example

```julia
area_triangle(10,10) # returns 50.0
area_triangle(-1,-2) # returns DomainError
area_triangle(1,-2)  # returns DomainError
area_triangle(-1,2)  # returns DomainError
```
"""
function area_triangle(base, height)
    if base < 0 || height < 0
        throw(DomainError("area_triangle() only accepts non-negative values"))
    end
    return (base * height) / 2
end

"""
    area_heron_triangle(side1, side2, side3)

Finds area of a triangle using heron's formula

# Example

```julia
area_heron_triangle(5,12,13) # returns 30.0
area_heron_triangle(-1,-2,1) # returns DomainError
area_heron_triangle(1,-2,1)  # returns DomainError
area_heron_triangle(-1,2,1)  # returns DomainError
```
"""
function area_heron_triangle(side1, side2, side3)
    if side1 < 0 || side2 < 0 || side3 < 0
        throw(DomainError("area_triangle_three_sides() only accepts non-negative values"))
    elseif side1 + side2 < side3 || side1 + side3 < side2 || side2 + side3 < side1
        throw(DomainError("Given three sides do not form a triangle"))
    end
    semi_perimeter = (side1 + side2 + side3) / 2
    area = sqrt(
        semi_perimeter
        * (semi_perimeter - side1)
        * (semi_perimeter - side2)
        * (semi_perimeter - side3)
    )
    return area
end

"""
    area_parallelogram(base, height)

Finds area of the parallelogram

# Example

```julia
area_parallelogram(10,20) # returns 200
area_parallelogram(-1,-2) # returns DomainError
area_parallelogram(1,-2)  # returns DomainError
area_parallelogram(-1,2)  # returns DomainError
```
"""
function area_parallelogram(base, height)

    if base < 0 || height < 0
        throw(DomainError("area_parallelogram() only accepts non-negative values"))
    end
    return base * height
end

"""
area_trapezium(base1,base2,height)

Finds area of the trapizium

# Example

```julia
area_trapezium(10, 20, 30)  # returns 450.0
area_trapezium(-1, -2, -3)  # returns DomainError
area_trapezium(-1, 2, 3)    # returns DomainError
area_trapezium(1, -2, 3)    # returns DomainError
area_trapezium(1, 2, -3)    # returns DomainError
area_trapezium(-1, -2, 3)   # returns DomainError
area_trapezium(1, -2, -3)   # returns DomainError
area_trapezium(-1, 2, -3)   # returns DomainError
```
"""
function area_trapezium(base1,base2,height)
    if base1 < 0 || base2 < 0 || height < 0
        throw(DomainError("area_trapezium() only accepts non-negative values"))
    end
    return 1 / 2 * (base1 + base2) * height
end

"""
    area_circle(radius)

Finds area of the circle

# Example

```julia
area_circle(20) # returns 1256.6370614359173
area_circle(-1) # returns DomainError
```
"""
function area_circle(radius)
    if radius < 0
        throw(DomainError("area_circle() only accepts non-negative values"))
    end
    return pi * radius ^ 2
end

"""
    area_ellipse(radius_x, radius_y)

Finds area of the ellipse

# Example

```julia
area_ellipse(10, 10) # returns 314.1592653589793
area_ellipse(10, 20) # returns 628.3185307179587
area_ellipse(1, -2) # returns DomainError
area_ellipse(-1, 2) # returns DomainError
area_ellipse(-1, -2) # returns DomainError
```
"""
function area_ellipse(radius_x, radius_y)
    if radius_x < 0 || radius_y < 0
        throw(DomainError("area_ellipse() only accepts non-negative values"))
    end
    return pi * radius_x * radius_y
end

"""
    area_rhombus(diagonal_1, diagonal_2)

Finds area of the rhombus

# Example

```julia
area_rhombus(10, 20) # returns 100.0
area_rhombus(-1,-2) # returns DomainError
area_rhombus(1,-2)  # returns DomainError
area_rhombus(-1,2)  # returns DomainError
```
"""
function area_rhombus(diagonal_1, diagonal_2)
    if diagonal_1 < 0 || diagonal_2 < 0
        throw(DomainError("area_rhombus() only accepts non-negative values"))
    end
    return 1 / 2 * diagonal_1 * diagonal_2
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
"""
    MinimumIntermodes <: AbstractImageBinarizationAlgorithm
    MinimumIntermodes()

    binarize([T,] img, f::MinimumIntermodes)
    binarize!([out,] img, f::MinimumIntermodes)

Under the assumption that the image histogram is bimodal the histogram is
smoothed using a length-3 mean filter until two modes remain. The binarization
threshold is then set to the minimum value between the two modes.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Arguments

The function argument is described in more detail below.

##  `img::AbstractArray`

The image that needs to be binarized. The image is automatically converted
to `Gray` in order to construct the requisite graylevel histogram.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(img, MinimumIntermodes())
```

# Reference

1. C. A. Glasbey, “An Analysis of Histogram-Based Thresholding Algorithms,” *CVGIP: Graphical Models and Image Processing*, vol. 55, no. 6, pp. 532–537, Nov. 1993. [doi:10.1006/cgip.1993.1040](https://doi.org/10.1006/cgip.1993.1040)
2. J. M. S. Prewitt and M. L. Mendelsohn, “THE ANALYSIS OF CELL IMAGES*,” *Annals of the New York Academy of Sciences*, vol. 128, no. 3, pp. 1035–1053, Dec. 2006. [doi:10.1111/j.1749-6632.1965.tb11715.x](https://doi.org/10.1111/j.1749-6632.1965.tb11715.x)
"""
struct MinimumIntermodes <: AbstractImageBinarizationAlgorithm end

function (f::MinimumIntermodes)(out::GenericGrayImage, img::GenericGrayImage)
    edges, counts = build_histogram(img,  256)
    t = find_threshold(HistogramThresholding.MinimumIntermodes(), counts[1:end], edges)
    @simd for i in CartesianIndices(img)
        out[i] = img[i] < t ? 0 : 1
    end
    out
end

(f::MinimumIntermodes)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
using Luna
import PyPlot: plt

#= In this example we simulate the effect of pump pulse dispersion on soliton self-compression
and RDW emission. To do this we use the "ϕ" and "propagator" keyword arguments to
prop_capillary to apply fixed third-order dispersion and variable glass insertion,
respectively.
By using the default settings otherwise, this example runs a mode-averaged simulation. This
is less accurate than multi-mode due to the presence of ionisation in the propagation, but
much faster, enabling a quick run of multiple simulations like this. =#

a = 75e-6 # core radius
flength = 0.6 # HCF length
gas = :Ne # gas species
pressure = 1.5 # gas pressure in bar

λ0 = 800e-9
τfwhm = 6e-15
energy = 100e-6

glassprop! = Fields.propagator_material(:SiO2)

thicknesses = [-2e-3, -1e-3, 0, 1e-3]

outputs = [prop_capillary(a, flength, gas, pressure; trange=1e-12, λlims=(150e-9, 4e-6),
                          λ0, τfwhm, energy,
                          ϕ=[0, 0, 0, -350e-45],
                          propagator=(Eω, grid) -> glassprop!(Eω, grid.ω, di, λ0))
            for di in thicknesses]
##
λ, _ = Processing.getIω(outputs[1], :λ)
Iλs = [Processing.getIω(oi, :λ)[2] for oi in outputs]

plt.figure()
for (d, Iλ) in zip(thicknesses, Iλs)
    plt.semilogy(λ*1e9, 1e-3*Iλ[:, end], label="$(d*1e3) mm")
end
plt.xlim(130, 1100)
plt.ylim(1e-5, 10)
plt.xlabel("Wavelength (nm)")
plt.ylabel("Spectral energy density (μJ/nm)")
plt.legend()


text:
function MultivariateSummaryStatistics(arg0::jint)
    return MultivariateSummaryStatistics((jint,), arg0)
end

function MultivariateSummaryStatistics(arg0::jint, arg1::jboolean)
    return MultivariateSummaryStatistics((jint, jboolean), arg0, arg1)
end

function add_value(obj::MultivariateSummaryStatistics, arg0::Vector{jdouble})
    return jcall(obj, "addValue", void, (Vector{jdouble},), arg0)
end

function clear(obj::MultivariateSummaryStatistics)
    return jcall(obj, "clear", void, ())
end

function equals(obj::MultivariateSummaryStatistics, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_covariance(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getCovariance", RealMatrix, ())
end

function get_dimension(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getDimension", jint, ())
end

function get_geometric_mean(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getGeometricMean", Vector{jdouble}, ())
end

function get_max(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getMax", Vector{jdouble}, ())
end

function get_mean(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getMean", Vector{jdouble}, ())
end

function get_min(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getMin", Vector{jdouble}, ())
end

function get_n(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getN", jlong, ())
end

function get_standard_deviation(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getStandardDeviation", Vector{jdouble}, ())
end

function get_sum(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getSum", Vector{jdouble}, ())
end

function get_sum_log(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getSumLog", Vector{jdouble}, ())
end

function get_sum_sq(obj::MultivariateSummaryStatistics)
    return jcall(obj, "getSumSq", Vector{jdouble}, ())
end

function hash_code(obj::MultivariateSummaryStatistics)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::MultivariateSummaryStatistics)
    return jcall(obj, "toString", JString, ())
end


text:
module TestExamples1DEulerGravity

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "1d")

@testset "Compressible Euler with self-gravity" begin
  @trixi_testset "elixir_eulergravity_eoc_test.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_eulergravity_eoc_test.jl"),
      l2   = [0.0002170799126638106, 0.0002913792848717502, 0.0006112320856262327],
      linf = [0.0004977401033188222, 0.0013594223337776157, 0.002041891084400227])
  end
end

end # module


text:
# The libbz2 Interfaces
# =====================

const WIN32 = Sys.iswindows() && Sys.WORD_SIZE == 32

mutable struct BZStream
    next_in::Ptr{UInt8}
    avail_in::Cint
    total_in_lo32::Cint
    total_in_hi32::Cint

    next_out::Ptr{UInt8}
    avail_out::Cint
    total_out_lo32::Cint
    total_out_hi32::Cint

    state::Ptr{Cvoid}

    bzalloc::Ptr{Cvoid}
    bzfree::Ptr{Cvoid}
    opaque::Ptr{Cvoid}
end

function BZStream()
    return BZStream(
        C_NULL, 0, 0, 0,
        C_NULL, 0, 0, 0,
        C_NULL,
        C_NULL, C_NULL, C_NULL)
end

# Action code
const BZ_RUN              = Cint(0)
const BZ_FLUSH            = Cint(1)
const BZ_FINISH           = Cint(2)

# Return code
const BZ_OK               = Cint( 0)
const BZ_RUN_OK           = Cint( 1)
const BZ_FLUSH_OK         = Cint( 2)
const BZ_FINISH_OK        = Cint( 3)
const BZ_STREAM_END       = Cint( 4)
const BZ_SEQUENCE_ERROR   = Cint(-1)
const BZ_PARAM_ERROR      = Cint(-2)
const BZ_MEM_ERROR        = Cint(-3)
const BZ_DATA_ERROR       = Cint(-4)
const BZ_DATA_ERROR_MAGIC = Cint(-5)
const BZ_IO_ERROR         = Cint(-6)
const BZ_UNEXPECTED_EOF   = Cint(-7)
const BZ_OUTBUFF_FULL     = Cint(-8)
const BZ_CONFIG_ERROR     = Cint(-9)

# Compressor
# ----------

function compress_init!(stream::BZStream,
                        blocksize100k::Integer,
                        verbosity::Integer,
                        workfactor::Integer)
    if WIN32
        return ccall(
            ("BZ2_bzCompressInit@16", libbz2),
            stdcall,
            Cint,
            (Ref{BZStream}, Cint, Cint, Cint),
            stream, blocksize100k, verbosity, workfactor)
    else
        return ccall(
            (:BZ2_bzCompressInit, libbz2),
            Cint,
            (Ref{BZStream}, Cint, Cint, Cint),
            stream, blocksize100k, verbosity, workfactor)
    end
end

function compress_end!(stream::BZStream)
    if WIN32
        return ccall(
            ("BZ2_bzCompressEnd@4", libbz2),
            stdcall,
            Cint,
            (Ref{BZStream},),
            stream)
    else
        return ccall(
            (:BZ2_bzCompressEnd, libbz2),
            Cint,
            (Ref{BZStream},),
            stream)
    end
end

function compress!(stream::BZStream, action::Integer)
    if WIN32
        return ccall(
            ("BZ2_bzCompress@8", libbz2),
            stdcall,
            Cint,
            (Ref{BZStream}, Cint),
            stream, action)
    else
        return ccall(
            (:BZ2_bzCompress, libbz2),
            Cint,
            (Ref{BZStream}, Cint),
            stream, action)
    end
end

# Decompressor
# ------------

function decompress_init!(stream::BZStream, verbosity::Integer, small::Bool)
    if WIN32
        return ccall(
            ("BZ2_bzDecompressInit@12", libbz2),
            stdcall,
            Cint,
            (Ref{BZStream}, Cint, Cint),
            stream, verbosity, small)
    else
        return ccall(
            (:BZ2_bzDecompressInit, libbz2),
            Cint,
            (Ref{BZStream}, Cint, Cint),
            stream, verbosity, small)
    end
end

function decompress_end!(stream::BZStream)
    if WIN32
        return ccall(
            ("BZ2_bzDecompressEnd@4", libbz2),
            stdcall,
            Cint,
            (Ref{BZStream},),
            stream)
    else
        return ccall(
            (:BZ2_bzDecompressEnd, libbz2),
            Cint,
            (Ref{BZStream},),
            stream)
    end
end

function decompress!(stream::BZStream)
    if WIN32
        return ccall(
            ("BZ2_bzDecompress@4", libbz2),
            stdcall,
            Cint,
            (Ref{BZStream},),
            stream)
    else
        return ccall(
            (:BZ2_bzDecompress, libbz2),
            Cint,
            (Ref{BZStream},),
            stream)
    end
end

# Error
# -----

struct BZ2Error <: Exception
    code::Cint
end

function bzerror(stream::BZStream, code::Cint)
    @assert code < 0
    throw(BZ2Error(code))
end


text:
using SemiclassicalOrthogonalPolynomials, ClassicalOrthogonalPolynomials, BandedMatrices, Test
import ClassicalOrthogonalPolynomials: orthogonalityweight, Weighted, associated, plotgrid

@testset "Two Band" begin
    @testset "TwoBandWeight" begin
        w = TwoBandWeight(0.6,1,2,3)
        @test w[0.7] == w[-0.7] == 0.7^6 * (0.7^2-0.6^2)^2 * (1-0.7^2)
        @test_throws BoundsError w[0.4]
        @test_throws BoundsError w[-0.4]
        @test_throws BoundsError w[1.1]
        @test copy(w) == w
        @test_broken sum(w)
    end
    @testset "Chebyshev case" begin
        ρ = 0.5
        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
        @test T[0.6,1:5] ≈ [T[0.6,k] for k=1:5]

        U = TwoBandJacobi(ρ, 1/2, 1/2, -1/2)
        @test U[0.6,1:5] ≈ [U[0.6,k] for k=1:5]

        @test copy(T) == T
        @test U ≠ T
        @test orthogonalityweight(T) == TwoBandWeight(ρ, -1/2, -1/2, 1/2)

        # bug
        @test !issymmetric(jacobimatrix(T)[1:10,1:10])
    end

    @testset "Associated" begin
        ρ = 0.5
        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
        Q = associated(T)
        x = axes(Q,1)
        @test 0 in x
        @test (Q[0.6,1:101]' * (Q[:,Base.OneTo(101)] \ exp.(x))) ≈ exp(0.6)
        @test (Q[-0.6,1:101]' * (Q[:,Base.OneTo(101)] \ exp.(x))) ≈ exp(-0.6)

        @test (Q * (Q \ exp.(x)))[0.6] ≈ exp(0.6)
    end

    @testset "Hilbert" begin
        ρ = 0.5
        w = TwoBandWeight(ρ, -1/2, -1/2, 1/2)
        T = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
        Q = associated(T)
        t = axes(w,1)
        x = axes(Q,1)
        H = inv.(x .- t')
        @test iszero(H*w)
        @test sum(w) ≈ π

        B = Q \ H * Weighted(T)
        @test B isa BandedMatrix

        @test Q[0.6,:]'* (B * (T\exp.(t))) ≈ -4.701116657130821 # Mathematica

        @test_broken H*TwoBandWeight(ρ, 1/2, 1/2, -1/2)
    end

    @testset "plotgrid" begin
        ρ = 0.5
        P = TwoBandJacobi(ρ, -1/2, -1/2, 1/2)
        @test all(x -> ρ ≤ abs(x) ≤ 1, plotgrid(P[:,1:5]))
    end

    @testset "associated transform error" begin
        a,b = 0.5, 0.9
        Vp = x -> 4x^3 - 20x
        Q = associated(TwoBandJacobi((a/b), -1/2, -1/2, 1/2))
        x = axes(Q,1)
        @test norm((Q[:,Base.OneTo(30)] \ Vp.(b*x))[1]) ≤ 1E-13
    end
end


text:
using ITensors,
      Test

@testset "Physics Sites" begin

  N = 10

  @testset "Spin Half sites" begin
    s = spinHalfSites(N)

    @test state(s[1],"Up") == s[1](1)
    @test state(s[1],"Dn") == s[1](2)
    @test_throws ArgumentError state(s[1],"Fake")

    Sz5 = op(s,"Sz",5)
    @test hasinds(Sz5,s[5]',s[5])

    @test_throws ArgumentError op(s, "Fake", 2)
    @test array(op(s,"S+",3),s[3]',s[3])  ≈ [ 0.0  1.0; 0.0  0.0]
    @test array(op(s,"S-",4),s[4]',s[4])  ≈ [ 0.0  0.0; 1.0  0.0]
    @test array(op(s,"Sx",2),s[2]',s[2])  ≈ [ 0.0  0.5; 0.5  0.0]
    @test array(op(s,"iSy",2),s[2]',s[2]) ≈ [ 0.0  0.5;-0.5  0.0]
    @test array(op(s,"Sy",2),s[2]',s[2])  ≈ [0.0  -0.5im; 0.5im  0.0]
    @test array(op(s,"Sz",2),s[2]',s[2])  ≈ [ 0.5  0.0; 0.0 -0.5]
    @test array(op(s,"projUp",2),s[2]',s[2])  ≈ [ 1.0  0.0; 0.0 0.0]
    @test array(op(s,"projDn",2),s[2]',s[2])  ≈ [ 0.0  0.0; 0.0 1.0]
    @test array(op(s,"Up",2),s[2])  ≈ [1.0,0.0]
    @test array(op(s,"Dn",2),s[2])  ≈ [0.0,1.0]
    @test array(op(s,"F",2),s[2]',s[2])  ≈ [ -1.0  0.0; 0.0 1.0]
  end

  @testset "Spin One sites" begin
    s = spinOneSites(N)

    @test state(s[1],"Up") == s[1](1)
    @test state(s[1],"0")  == s[1](2)
    @test state(s[1],"Dn") == s[1](3)
    @test_throws ArgumentError state(s[1],"Fake")

    Sz5 = op(s,"Sz",5)
    @test hasinds(Sz5,s[5]',s[5])

    @test_throws ArgumentError op(s, "Fake", 2)
    @test array(op(s,"S+",3),s[3]',s[3]) ≈ [ 0 √2 0; 0 0 √2; 0 0 0]
    @test array(op(s,"S-",3),s[3]',s[3]) ≈ [ 0 0 0; √2 0 0; 0.0 √2 0]
    @test array(op(s,"Sx",3),s[3]',s[3]) ≈ [ 0 1/√2 0; 1/√2 0 1/√2; 0 1/√2 0]
    @test array(op(s,"iSy",3),s[3]',s[3]) ≈ [ 0 1/√2 0; -1/√2 0 1/√2; 0 -1/√2 0]
    @test array(op(s,"Sy",3),s[3]',s[3]) ≈ [ 0 -1/√2im 0; +1/√2im 0 -1/√2im; 0 +1/√2im 0]
    @test array(op(s,"Sz",2),s[2]',s[2]) ≈ [1.0 0 0; 0 0 0; 0 0 -1.0]
    @test array(op(s,"Sz2",2),s[2]',s[2]) ≈ [1.0 0 0; 0 0 0; 0 0 +1.0]
    @test array(op(s,"Sx2",2),s[2]',s[2]) ≈ [0.5 0 0.5;0 1.0 0;0.5 0 0.5]
    @test array(op(s,"Sy2",2),s[2]',s[2]) ≈ [0.5 0 -0.5;0 1.0 0;-0.5 0 0.5]
    @test array(op(s,"projUp",2),s[2]',s[2]) ≈ [1.0 0 0;0 0 0;0 0 0]
    @test array(op(s,"projZ0",2),s[2]',s[2]) ≈ [0 0 0;0 1.0 0;0 0 0]
    @test array(op(s,"projDn",2),s[2]',s[2]) ≈ [0 0 0;0 0 0;0 0 1.0]
    @test array(op(s,"XUp",2),s[2]) ≈ [0.5,im*√2,0.5]
    @test array(op(s,"XZ0",2),s[2]) ≈ [im*√2,0,-im*√2]
    @test array(op(s,"XDn",2),s[2]) ≈ [0.5,-im*√2,0.5]
  end

  @testset "Electron sites" begin
    s = electronSites(N)

    @test state(s[1],"0")    == s[1](1)
    @test state(s[1],"Up")   == s[1](2)
    @test state(s[1],"Dn")   == s[1](3)
    @test state(s[1],"UpDn") == s[1](4)
    @test_throws ArgumentError state(s[1],"Fake")

    Nup5 = op(s,"Nup",5)
    @test hasinds(Nup5,s[5]',s[5])

    @test_throws ArgumentError op(s, "Fake", 2)
    Nup3 = array(op(s,"Nup",3),s[3]',s[3])
    @test Nup3 ≈ [0. 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 1]
    Ndn3 = array(op(s,"Ndn",3),s[3]',s[3])
    @test Ndn3 ≈ [0. 0 0 0; 0 0 0 0; 0 0 1 0; 0 0 0 1]
    Ntot3 = array(op(s,"Ntot",3),s[3]',s[3])
    @test Ntot3 ≈ [0. 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 2]
    Cup3 = array(op(s,"Cup",3),s[3]',s[3])
    @test Cup3 ≈ [0. 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0]
    Cdagup3 = array(op(s,"Cdagup",3),s[3]',s[3])
    @test Cdagup3 ≈ [0. 0 0 0; 1 0 0 0; 0 0 0 0; 0 0 1 0]
    Cdn3 = array(op(s,"Cdn",3),s[3]',s[3])
    @test Cdn3 ≈ [0. 0 1 0; 0 0 0 1; 0 0 0 0; 0 0 0 0]
    Cdagdn3 = array(op(s,"Cdagdn",3),s[3]',s[3])
    @test Cdagdn3 ≈ [0. 0 0 0; 0 0 0 0; 1 0 0 0; 0 1 0 0]
    F3 = array(op(s,"F",3),s[3]',s[3])
    @test F3 ≈ [1. 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 1]
    Fup3 = array(op(s,"Fup",3),s[3]',s[3])
    @test Fup3 ≈ [1. 0 0 0; 0 -1 0 0; 0 0 1 0; 0 0 0 -1]
    Fdn3 = array(op(s,"Fdn",3),s[3]',s[3])
    @test Fdn3 ≈ [1. 0 0 0; 0 1 0 0; 0 0 -1 0; 0 0 0 -1]
    Sz3 = array(op(s,"Sz",3),s[3]',s[3])
    @test Sz3 ≈ [0. 0 0 0; 0 0.5 0 0; 0 0 -0.5 0; 0 0 0 0]
    Sx3 = array(op(s,"Sx",3),s[3]',s[3])
    @test Sx3 ≈ [0. 0 0 0; 0 0 0.5 0; 0 0.5 0 0; 0 0 0 0]
    Sp3 = array(op(s,"S+",3),s[3]',s[3])
    @test Sp3 ≈ [0. 0 0 0; 0 0 1 0; 0 0 0 0; 0 0 0 0]
    Sm3 = array(op(s,"S-",3),s[3]',s[3])
    @test Sm3 ≈ [0. 0 0 0; 0 0 0 0; 0 1 0 0; 0 0 0 0]
    Sem = array(op(s,"Emp",3),s[3])
    @test Sem ≈ [1.0; 0.0; 0.0; 0.0]
    Sup = array(op(s,"Up",3),s[3])
    @test Sup ≈ [0.0; 1.0; 0.0; 0.0]
    Sdn = array(op(s,"Dn",3),s[3])
    @test Sdn ≈ [0.0; 0.0; 1.0; 0.0]
    Supdn = array(op(s,"UpDn",3),s[3])
    @test Supdn ≈ [0.0; 0.0; 0.0; 1.0]
  end

  @testset "tJ sites" begin
    s = tJSites(N)

    @test state(s[1],"0")    == s[1](1)
    @test state(s[1],"Up")   == s[1](2)
    @test state(s[1],"Dn")   == s[1](3)
    @test_throws ArgumentError state(s[1],"Fake")

    @test_throws ArgumentError op(s, "Fake", 2)
    Nup_2 = op(s,"Nup",2)
    @test Nup_2[2,2] ≈ 1.0
    Ndn_2 = op(s,"Ndn",2)
    @test Ndn_2[3,3] ≈ 1.0
    Ntot_2 = op(s,"Ntot",2)
    @test Ntot_2[2,2] ≈ 1.0
    @test Ntot_2[3,3] ≈ 1.0
    Cup3 = array(op(s,"Cup",3),s[3]',s[3])
    @test Cup3 ≈ [0. 1 0; 0 0 0; 0 0 0]
    Cdup3 = array(op(s,"Cdagup",3),s[3]',s[3])
    @test Cdup3 ≈ [0 0 0; 1. 0 0; 0 0 0]
    Cdn3 = array(op(s,"Cdn",3),s[3]',s[3])
    @test Cdn3 ≈ [0. 0. 1; 0 0 0; 0 0 0]
    Cddn3 = array(op(s,"Cdagdn",3),s[3]',s[3])
    @test Cddn3 ≈ [0 0 0; 0. 0 0; 1 0 0]
    FP3 = array(op(s,"FP",3),s[3]',s[3])
    @test FP3 ≈ [1.0 0. 0; 0 -1.0 0; 0 0 -1.0]
    Fup3 = array(op(s,"Fup",3),s[3]',s[3])
    @test Fup3 ≈ [1.0 0. 0; 0 -1.0 0; 0 0 1.0]
    Fdn3 = array(op(s,"Fdn",3),s[3]',s[3])
    @test Fdn3 ≈ [1.0 0. 0; 0 1.0 0; 0 0 -1.0]
    Sz3 = array(op(s,"Sz",3),s[3]',s[3])
    @test Sz3 ≈ [0.0 0. 0; 0 0.5 0; 0 0 -0.5]
    Sx3 = array(op(s,"Sx",3),s[3]',s[3])
    @test Sx3 ≈ [0.0 0. 0; 0 0 1; 0 1 0]
    Sp3 = array(op(s,"Splus",3),s[3]',s[3])
    @test Sp3 ≈ [0.0 0. 0; 0 0 1.0; 0 0 0]
    Sm3 = array(op(s,"Sminus",3),s[3]',s[3])
    @test Sm3 ≈ [0.0 0. 0; 0 0 0; 0 1.0 0]
    Up3 = array(op(s,"Up",3),s[3])
    @test Up3 ≈ [0.0; 1.0; 0]
    Dn3 = array(op(s,"Dn",3),s[3])
    @test Dn3 ≈ [0.0; 0.0; 1.0]
    Em3 = array(op(s,"Emp",3),s[3])
    @test Em3 ≈ [1.0; 0.0; 0.0]
  end

end

const MySite = makeTagType("MySite")

@testset "Custom Site Tag Type" begin

  function ITensors.op(::MySite,s::Index,opname::AbstractString)
    Op = ITensor(s,s')
    if opname=="MyOp"
      Op[s(1),s'(1)] = 11
      Op[s(1),s'(2)] = 12
      Op[s(2),s'(1)] = 21
      Op[s(2),s'(2)] = 22
    end
    return Op
  end

  function ITensors.state(::MySite,statename::AbstractString)
    if statename == "One"
      return 1
    elseif statename == "Two"
      return 2
    end
  end

  i = Index(2,"MySite")

  expectedOp = ITensor(i,i')
  expectedOp[i(1),i'(1)] = 11
  expectedOp[i(1),i'(2)] = 12
  expectedOp[i(2),i'(1)] = 21
  expectedOp[i(2),i'(2)] = 22
  @test norm(op(i,"MyOp")-expectedOp) < 1E-10

  @test state(i,"One") == i(1)
  @test state(i,"Two") == i(2)
end


text:
module ComputeManager
using Cairn.Client
using Excavator: Job
function run(job::t)
    # TODO keep track of ComputeNodes

    # start
        # run the job
end
immutable MapJob
    dfs::Vector{Cairn.Client.t}
    executable::Job.Key
    files::Cairn.RPC.KeySpan
    results_dir::Job.Key
    done::Job.Key # Special 0-replication file!

end

immutable ReduceJob
    dfs::Vector{Cairn.Client.t}
    executable::Job.Key
    files::Cairn.RPC.KeySpan
    results_dir::Job.Key
end


text:
"""
    channelify(m::AbstractArray{CT,N}) where {CT <: Colorant, N}
Channelify function.

Input should be in ``WHN`` order ``\\implies (*, batch)``.
Output will be in the order ``(*, channels, batch)``.

# Examples
```jldoctest; setup = :(using ImageCore, DiffImages)
julia> input_size = (16, 16, 2)
(16, 16, 2)

julia> size(channelify(HSV.(rand(BGR,input_size...))))
(16, 16, 3, 2)
```
"""
function channelify(m::AbstractArray{CT,N}) where {CT <: Colorant, N}
    e = eltype(m)
    m = channelview(m)
    if e <: AbstractGray
        m = reshape(m, 1, size(m)...)
    end
    t = ntuple(identity, ndims(m))
    m = permutedims(m, (t[2:end-1]...,1,t[end]))
    return m
end

@adjoint function channelview(x::AbstractArray{T,N}) where {T, N}
    e = eltype(x)
    y = channelview(x)
    function pullback(Δ)
        return (collect(colorview(e,Δ)),)
    end
    return (y, pullback)
end

@adjoint function colorview(T, x)
    y = colorview(T,x)
    function pullback(Δ)
        return (nothing, channelview(Δ))
    end
    return (y, pullback)
end

# adjoint for (::Colorant{T,N})(x::Real)
for f in (:HSV,:AHSV,:HSVA,
          :Gray,:AGray,:GrayA,
          :HSL,:AHSL,:HSLA,
          :RGB,:ARGB,:RGBA,
          :BGR,:ABGR,:BGRA,
          :XYZ,:AXYZ,:XYZA,
          :xyY,:AxyY,:xyYA,
          :Lab,:ALab,:LabA,
          :Luv,:ALuv,:LuvA,
          :LCHab,:ALCHab,:LCHabA,
          :LCHuv,:ALCHuv,:LCHuvA,
          :DIN99,:ADIN99,:DIN99A,
          :LMS,:ALMS,:LMSA,
          :YIQ,:AYIQ,:YIQA)
    @eval @adjoint function $f(x::Real...)
        y = $f(x::Real...)
        function pull(Δ...)
            return (Δ...,)
        end
        return (y, pull)
    end
end

# Constructor adjoint
# function ChainRules.rrule(::Type{HSL{T}}, x, y, z) where T <: AbstractFloat
#     β = HSL{T}(x, y, z)
#     function Lab_pullback(Δ)
#         @show Δ x y z β
#         return (ChainRules.NoTangent(),Δ.h,Δ.s,Δ.l)
#     end
#     return (β, Lab_pullback)
# end

"""
    colorify(color::Type{CT}, m::AbstractArray) where CT <: Colorant
Colorify function.

Expecting an input of the type ``(*, channels, batch)``.
Converts the array to the `color` specified.

# Examples
```jldoctest; setup = :(using ImageCore, DiffImages)
julia> input_size = (25, 25, 3, 7)
(25, 25, 3, 7)

julia> size(colorify(HSV, rand(input_size...)))
(25, 25, 7)
```
"""
function colorify(color::Type{CT}, m::AbstractArray) where CT <: Colorant
    t = ntuple(identity, ndims(m))
    m = permutedims(m, (t[end-1],t[1:end-2]...,t[end]))
    if color <: AbstractGray
        m = dropdims(m; dims=1)
    end
    m = colorview(color, m)
    return m
end
# TODO: adjoints for `colorview(T, gray1, gray2...)` are not added yet.


text:
"Generic PDDL action definition."

GenericAction(term::Term, precond::Term, effect::Term) =
    GenericAction(term.name, get_args(term), Symbol[], precond, effect,
        Compound(:cont_precond, Term[]), Compound(:cont_effect, Term[]),
        Compound(:dynamics, Term[]))

Base.:(==)(a1::GenericAction, a2::GenericAction) = (a1.name == a2.name &&
    Set(a1.args) == Set(a2.args) && Set(a1.types) == Set(a2.types) &&
    a1.precond == a2.precond && a1.effect == a2.effect)

get_name(action::GenericAction) = action.name

get_argvars(action::GenericAction) = action.args

get_argtypes(action::GenericAction) = action.types

get_precond(action::GenericAction) = action.precond

function get_precond(action::GenericAction, args)
    subst = Subst(k => v for (k, v) in zip(action.args, args))
    return substitute(action.precond, subst)
end

get_effect(action::GenericAction) = action.effect

get_action(domain::GenericDomain, name::Symbol) = domain.actions[name]

function get_effect(action::GenericAction, args)
    subst = Subst(k => v for (k, v) in zip(action.args, args))
    return substitute(action.effect, subst)
end


text:
"""
    Basic implementation of Weibull function with 2-parameters
    κ: shape parameter, κ > 0
    λ: scale parameter, λ > 0
    t: current time (generally calculated with wb_t function)

"""
weibull_f(λ,κ,t) = t >= 0 ? (κ / λ) * ((t/ λ)^(κ - 1))*(ℯ^-(t/λ)^κ) : 0

"""
 Weibull scaling function
 ttf: expected time-to-failure
 γ: factor to tune max y (yₙ) value, e.g. for max y=1, with ttf =100,
 λ=0.5, κ=05 => γ=8

 TODO: check if γ might be calculated from ttf, λ and κ
"""
wb_t(γ,ttf,t) = (t%ttf)*γ/ttf

"""
 Basic logistic function
 http://reliawiki.org/index.php/Logistic

"""

log_f(b,k,t) = 1 / ( 1+ b*(ℯ^(-k*t)))

"""
Add random normally distributed  noise to a given value.
Ensures that always:  0 ⋜  v + ϵ ⋜ y, for a random noise ϵ.

    v: original value without noise
    ϵₛ: noise standard deviation
    ϵ₀: noise μ
    y:  max value the noised value can take
"""
function nnoised(v,ϵₛ,ϵ₀,y)
    ϵ = 0
    #force only noised values within range
    while true
        ϵ = rand(Normal(ϵ₀,ϵₛ),1)[1]
        (( v + ϵ ) > y || ( v + ϵ ) < 0)  || break
    end
    return v + ϵ
end

"""
Add random normally distributed noise between 0 and y=1.
See nnoised(v,ϵₛ,ϵ₀,y).
"""

nnoised(v,ϵₛ) = nnoised(v,ϵₛ,0,1)

"""
    Calculates remaining ttf(time-to-failure) for a given expected ttf
    and current t
    ttf: expected time-to-failure
    t: current time

"""
exp_t(ttf,t) = ttf - (t%ttf)

"""
    Basic implementation of exponential deterioration function
    according to:

    * T. Wang, J. Yu, D. Siegel, and J. Lee, “2008 International Conference on Prognostics and Health Management, PHM2008,”
      2008 Int. Conf. Progn. Heal. Manag. PHM 2008, 2008.

      * A. Salvador Palau, M. Dhada, and A. Parlikad, “Multi-Agent System architectures for collaborative prognostics,”
      J. Intell. Manuf., 2019.

    a: max expected value of condition/health indicator (i.e. 1 if normalised)
    b: curvature parameter
    t: remaining time-to-failure (generally calculated with exp_t function)

"""
exp_f(a,b,t) = a*(1.0 - ℯ^(-(b)*t))

"""
    Scales weibull input t value considering shifting
    for long term behaviour

    ttf: expected time-to-failure
    γ: factor to tune max y (yₙ) value, e.g. for max y=1, with ttf =100,
    λ=0.5, κ=05 => γ=8
    d: units to shift to the right (near TTF)

"""
wb_ts(γ,d,ttf,t) = ((t+d)%ttf)*γ/ttf

"""
    Scales exponential function considering shifting for long term
    behaviour.

    ttf: expected time-to-failure
    t: value to scale
    d: units to shift to the right (near TTF)
"""
exp_ts(d,ttf,t) = ttf - (t%ttf) -d

"""
    Scales logistic function considering shifting for long term
    behaviour.

    ttf: expected time-to-failure
    t: value to scale
    d: units to shift to the right (near TTF)
"""

log_ts(d,ttf,t) = (t+d)%ttf

"""
Introduces a variation on the x value according to the phase (ϕ)
Used concretely for long-term TTF, we assume:
    * An asset has expected TTF (x)
    * After TTF the asset can be repaired or maintained and then will
      operate for a number of cycles/phases.
    * The current cycle/phase is given by ϕ
    * After every phase ϕ, the asset reduces its TTF by Δᵩ
"""
xᵩ(x,ϕ,Δᵩ) =  x * (1 -  (ϕ * Δᵩ))

"""
    Determines criteria for valid values to TTF for
    the exponential function.

        t₀: unused, defined to keep uniformity of function
            (TODO: check if this can be removed)
        t: value to evaluate
"""
exp_c(t₀,t) = t > 0

"""
    Determines criteria for valid values to TTF for
    the weibull function.

        t₀: initial value for the cycle/phase
        t: value to evaluate
"""

wb_c(t₀,t) = t >= t₀

"""
    Determines criteria for valid values to TTF for
    the logistics function.

        t₀: initial value for the cycle/phase
        t: value to evaluate
"""
log_c(t₀,t) = t >= t₀

"""
Linear degradation
m: degradation slope
b: expected life
t: current time
"""
lineal_d(m,b,t) = begin
    v = -m * t + b
    return v > 0 ? v : 0
end

"""
    Generates the time series until failure of an asset, given
    the provided parameters:

    ttf: expected time-to-failure
    n: maximum number of cycles the asset can last
    Δᵩ: Delta of number of time steps that are reduced as the asset
        reaches a new cycle/phase
    ϵₛ: Standard deviation of the gaussian noise
    f: function to generate values until ttf
    params: params required by the f
    tsf: function to scale t for the given function
    tse: extra parameter required by tsf
    f_c: function that defines criteria for a valid value in the time series.
"""
function values_f(ttf,n,Δᵩ,ϵₛ,f,params,tsf,tse_params,f_c)
    ϕ = 0
    t₀ = 0
    vs = []
    for i=0:n-1
        t = i
        ϕ₀ = ϕ
        ϕ = floor(t/ttf)
        rttf = xᵩ(ttf,ϕ,Δᵩ)
        diff = ttf - rttf
        t = tsf(tse_params...,diff,ttf,t)
        t₀ = i == 1 ? t : ϕ > ϕ₀ ? t : t₀
        remt = ϕ > ϕ₀ ? (n - i) >= rttf && rttf > 0 : true
        if f_c(t₀,t) && remt
            #push!(vs_f,(t,t₀,ϕ,ϕ₀,rmng,rttf,(n - i),f(params...,t)))
            push!(vs,nnoised(f(params...,t),ϵₛ))
        end
    end
    l = length(vs) + 1
    for i=l:n
        push!(vs,0.0)
    end
    return vs
end

"""
 It illustrates how to use the other TTE functions
"""
function generate_ttf_series()
    n=1000
    ttf = 200
    λ = 1.0 #0.5
    γ = 6.0 #16
    k = 1.0 #0.5
    #p = plot(title="",legend=false)#,xlims=[0,1000])
    ϵₛ = 0.05
    a = 1.0
    b = 0.05
    Δᵩ = 0.05#0.009
    maxᵩ = 9

    vs = []
    funs = [
            (exp_f,(a,b),exp_ts,(),exp_c),
            (weibull_f,(λ,k),wb_ts,(γ),wb_c),
            (log_f,(50,0.1),log_ts,(),log_c)
            ]

    for fs in funs
        push!(vs,values_f(ttf,n,Δᵩ,ϵₛ,fs...))
    end

    # for d=1:length(vs)
    #     p = plot!(p, vs[d])
    # end
    # for c=ttf:ttf:n
    #     p = plot!(p,[c,c],[0,1.0], c=:red, ls=:dash, w=2)
    # end
    # p
end

"""
    It adds downtime to a series
        ot: downtime (time steps)

"""

function add_downtime_to_series(ttf,ot,s_series)

    #log_info(" RECEIVED adding downtime ==>$(s_series)<==")

    # split s_series in phases/cycles of ttf
    phased_s = [s_series[i:min(i + ttf - 1, end)] for i in 1:ttf:length(s_series)]

    ph_ot = Vector{Float64}()

    for ph in phased_s
        if length(ph) < ttf
            push!(ph_ot,ph...)
        else
            push!(ph_ot,vcat(ph,zeros(Float64,ot))...)
        end
    end

    # log_info(" adding downtime $(size(ph_ot)) --- $(phased_s)")
    return ph_ot
end

"""
 It generates the sensor time series for a given asset
    ttf:time-to-failure
    ot: off-time
"""
function generate_sensor_series(ttf,n,Δᵩ,ϵₛ,ot,funs)
    vs = zeros(Float64,length(funs),n)
    #log_info("received generate sensor ==> $ot ")
    dim_dtvs = n + Int64(floor(n/ttf)) * ot  #n%ttf > 0 && n >= ttf ? n+ot+1 : n+ot

    dtvs = zeros(Float64,length(funs),dim_dtvs)

    for i=1:length(funs)
        vs[i,:] = values_f(ttf,n,Δᵩ,ϵₛ,funs[i]...)
    end

    for si=1:length(funs)
        s = vs[si,:]
        dtvs[si,:] = add_downtime_to_series(ttf,ot,s)
    end

    return dtvs
end

function generate_rul_series(ttf,Δᵩ,n,ot)
    rul = [  i%xᵩ(ttf,floor((i-1)/ttf),Δᵩ) == 0 ? 0 : round(xᵩ(ttf,floor((i-1)/ttf),Δᵩ) - i%xᵩ(ttf,floor((i-1)/ttf),Δᵩ)) for i=1:n ]
    #rul = [ xᵩ(ttf,floor(i/ttf),Δᵩ) - i%xᵩ(ttf,floor(i/ttf),Δᵩ) + 1 for i=1:n ]
    # print("PRE RUL: $(size(rul))")
    rul = add_downtime_to_series(ttf,ot,rul)
    # print("RUL: $(size(rul))")
    return rul
end


text:
"""
    mehrotra!(mechanism; opts)

    interior-point solver for simulation-step feasibility problem

    mechanism: Mechanism
    opts: SolverOptions
"""
function mehrotra!(mechanism::Mechanism; opts=SolverOptions())
	reset!.(mechanism.contacts, scale=1.0) # resets the values of s and γ to the scaled neutral vector; TODO: solver option
	reset!.(mechanism.joints,   scale=1.0) # resets the values of s and γ to the scaled neutral vector; TODO: solver option

	status = :failed
    mechanism.μ = 0.0
	μtarget = 0.0
	no_progress = 0
    undercut = opts.undercut
    α = 1.0

	initialize!.(mechanism.contacts) # TODO: redundant with resetVars--remove
    set_entries!(mechanism) # compute the residual

    bvio = bilinear_violation(mechanism) # does not require to apply set_entries!
    rvio = residual_violation(mechanism) # does not require to apply set_entries!

	opts.verbose && solver_header()
    for n = Base.OneTo(opts.max_iter)
        opts.verbose && solver_status(mechanism, α, rvio, bvio, n, μtarget, undercut)

        ((rvio < opts.rtol) && (bvio < opts.btol)) && (status=:success; break)
		(n == opts.max_iter) && (opts.verbose && (@warn "failed mehrotra"))

        # affine search direction
		μ = 0.0
		pull_residual!(mechanism)               # store the residual inside mechanism.residual_entries
        ldu_factorization!(mechanism.system)    # factorize system, modifies the matrix in place
        pull_matrix!(mechanism)                 # store the factorized matrix inside mechanism.matrix_entries
        ldu_backsubstitution!(mechanism.system) # solve system, modifies the vector in place

		αaff = cone_line_search!(mechanism; τort=0.95, τsoc=0.95, scaling=false) # uses system.vector_entries which holds the search drection
		ν, νaff = centering!(mechanism, αaff)
		σcentering = clamp(νaff / (ν + 1e-20), 0.0, 1.0)^3

		# corrected search direction
		μtarget = max(σcentering * ν, opts.btol / undercut)
		mechanism.μ = μtarget
		correction!(mechanism) # update the residual in mechanism.residual_entries

		push_residual!(mechanism)               # cache residual + correction
        push_matrix!(mechanism)                 # restore the factorized matrix
        ldu_backsubstitution!(mechanism.system) # solve system

		τ = max(0.95, 1 - max(rvio, bvio)^2) # τ = 0.95
		α = cone_line_search!(mechanism; τort=τ, τsoc=min(τ, 0.95), scaling=false) # uses system.vector_entries which holds the corrected search drection

		# steps taken without making progress
		rvio_, bvio_ = line_search!(mechanism, α, rvio, bvio, opts)

        # evaluate progress
		made_progress = (!(rvio_ < opts.rtol) && (rvio_ < 0.8rvio)) || (!(bvio_ < opts.btol) && (bvio_ < 0.8bvio)) # we only care when progress is made while the tolerance is not met
		made_progress ? no_progress = max(no_progress - 1, 0) : no_progress += 1
		rvio, bvio = rvio_, bvio_
		(no_progress >= opts.no_progress_max) && (undercut *= opts.no_progress_undercut)

		# update solution
        update!.(mechanism.bodies)
        update!.(mechanism.joints)
        update!.(mechanism.contacts)

		# recompute Jacobian and residual
        set_entries!(mechanism)
    end

    return status
end

function solver_status(mechanism::Mechanism, α, rvio, bvio, n, μtarget, undercut)
    fv = full_vector(mechanism.system)
    Δvar = norm(fv, Inf)
    fM = full_matrix(mechanism.system)
    fΔ = fM \ fv
    Δalt = norm(fΔ, Inf)
    res = norm(fv, Inf)
	println(
        n,
        "   ", scn(bvio, digits=0),
        "   ", scn(rvio, digits=0),
        "   ", scn(α, digits=0),
        "   ", scn(μtarget, digits=0),
        "   ", scn(res, digits=0),
        "   ", scn(Δvar, digits=0),
        # "   ucut", scn(undercut),
        )
end

function solver_header()
	println("                                                 ")
	println("n    bvio    rvio     α       μ     |res|∞   |Δ|∞")
	println("–––––––––––––––––––––––––––––––––––––––––––––––––")
end


text:
# julia --trace-compile=pms_trace.jl
# PackageCompiler.restore_default_sysimage()

# DO NOT compile GameOne! Julia will not compile local dev pkg changes??

import Pkg

Pkg.add(url="https://botm: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@github.com/dustyirwin/GameOne.jl")
Pkg.instantiate()

using PackageCompiler

project_symbols = [
    :Colors
    :DataStructures
    :GameOne
    :HTTP
    :ImageFiltering
    :ImageMagick
    :ImageTransformations
    :Images
    :JSON
    :Pluto
    :PlutoUI
    :Reexport
    :Revise
    :Rotations
    :ShiftedArrays
    :Dates
    :Logging
    :Pkg
    :Random
    :SHA
    :Serialization
]

create_sysimage(project_symbols,
    precompile_statements_file = "BotM-trace.jl",
    sysimage_path = "botm.so"
)

#=
PackageCompiler.audit_app(".")  # passing!

create_app(
    "../PlaymatSimulator.jl",
    "../PS.compiled",
    filter_stdlibs = false,
    #precompile_statements_file="tmp/trace.jl",
    incremental = true,
    #force=true
)

# UUID for GameOne b831c3d0-16f1-4eaa-ad23-1d0424ad4597
# UUID for PlaymatSimulator e5f19dc6-78e9-41ab-ba95-fec0547b5fb9
using Dates
using Random
using UUIDs
uuid4(MersenneTwister(Int(floor(datetime2unix(now())))))

=#

max_stars_repo_name: dustyirwin/PlaymatSimulator

max_issues_repo_name: dustyirwin/PlaymatSimulator

max_forks_repo_name: dustyirwin/PlaymatSimulator


text:
__precompile__() ## Required in order to segfault

module BigsSegfault
const α = parse(BigFloat,"1.01")
f = (t,u) -> (α*u)
## Other test
function fmaker(α2=1.01)
  f(t,u) = α2*u
  return f
end
export f,α,fmaker
end # module


text:
"""
    log1mexp(x)

Compute `log(1 - exp(x))` accurately. See
"**Accurately Computing log(1 - exp(.)) – Assessed by Rmpfr**" by Martin Mächler
(2012) for details:
https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
"""
function log1mexp(x::Real)
  if x > -log(2)
    log(-expm1(x))
  else
    log1p(-exp(x))
  end
end

"""
    logexpm1(x)

Compute `log(exp(x) - 1)` accurately. Thresholds were obtained similarly to
"**Accurately Computing log(1 - exp(.)) – Assessed by Rmpfr**" by Martin Mächler
(2012):
https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
"""
function logexpm1(x::Real)
  if x >= 37
    x
  elseif x >= 19
    x - exp(-x)
  else
    log(expm1(x))
  end
end

"""
    log1pexp(x)

Compute `log(1 + exp(x))` accurately. See
"**Accurately Computing log(1 - exp(.)) – Assessed by Rmpfr**" by Martin Mächler
(2012) for details:
https://cran.r-project.org/web/packages/Rmpfr/vignettes/log1mexp-note.pdf
"""
function log1pexp(x::Real)
  if x <= -37
    exp(x)
  elseif x <= 18
    log1p(exp(x))
  elseif x <= 33.3
    x + exp(-x)
  else
    x
  end
end


text:
module LargeMatrixEquations

export kpik, rksm, lp_lradi, lp_para

"""
    `kpik(A,B,E=1;<keyword arguments>)`

Julia code for K-PIK (Krylov-plus-inverted-Krylov)

Based on kpik.m avalible from V. Simoncini's website
(http://www.dm.unibo.it/~simoncin/software.html)
and includes many of the same comments and much of the same description.

Approximately solve

     A X E + E X A' + BB' = 0

by means of the extended Krylov subspace method

ARGUMENTS

    'A' : coeff matrix, A < 0

    'B' : factor of rhs,   nxk matrix with k << n

NAMED ARGUMENTS

    'E' : coeff matrix, spd, Defult: 1

    'LE' : lower triang factor of coeff matrix, Defult: 1 Note: This is an
    optional argument, if not provided it will be set to cholfact(E)[:L]

    'm' : max space dimension, Defult: 100

    'tol' : stopping tolerance based on the backwards error,
    with stopping criterion

      ||LE\\A X LE  + LE' X A'/LE'-LE\\BB'/LE'||
    -----------------------------------------------   < tol
    ||LE\\BB'/LE'|| + ||E^{-1}|| ||A|| ||LE'X LE ||

    computed in a cheap manner. Defult: 1e-9

Output

    'Z' : solution factor s.t.  X = Z Z'

    'er2' : history of scaled residual, as above

Comments

    1. The projected solution is computed at each iteration
    As an alternative, a periodic computation could be considered.

    2. This code performs a factorization of A. As an alternative,
    iterative solves could be considered.

    3. To solve

     A X + X A' + BB' = 0

     Use kpik(A,B) as E is set to 1 by Defult

If you use this code, please cite the following article:

    V. Simoncini
    A new iterative method for solving large-scale Lyapunov matrix equations,
    SIAM J.  Scient. Computing, v.29, n.3 (2007), pp. 1268-1288.
"""
function kpik(A,B,E=1;LE=1,m=100,tol=1e-9,tolY=1e-12,infoV=true)

  @assert(isdefined(:vecnorm),"Your julia version is too old. vecnorm() not defined")

  infoV && tic()

  #### Check if we solve the general lyapunov equation and user did not provide LE
  if E != 1 && LE == 1
    #### If this is the case, calculate LE
    if issparse(E)
      LE = try
            sparse(cholfact(E,perm=1:size(E,1))[:L])
           catch
            error("E must be SPD")
          end
    else
      LE = try
            cholfact(E)[:L]
          catch
            error("E must be SPD")
          end
    end
  end

  const rhs=LE\B
  const nrmb=vecnorm(rhs)^2
  const nrma=vecnorm(A)
  const sqrt2=sqrt(2)
  er2=zeros(m,1)

  const n,sh=size(rhs)

  Y=[]
  odds=[]
  er2=[]

  if E == 1
    const condestE=1
    const singE=condestE/vecnorm(eye(n))
  elseif (vecnorm(E-speye(n))>1e-14)
    const condestE=cond(full(E))
    const singE=condestE/vecnorm(E)
  else
    const singE=1
  end

  if norm(A-A',1)<1e-14
     if issparse(A)
       cfA = cholfact(-A,perm=1:size(A,1))
       const LA = sparse(cfA[:L])
       const UA = -LA'
     else
       const UA = chol(-A)
       const LA = -UA'
     end
     infoV && println("A sym. Completed Chol factorization\n")
     const k_max =2
   else
     luA=lufact(full(A),Val{false})
     const LA = luA[:L]
     const UA = luA[:U]

     infoV && println("A nonsym. Completed LU factorization\n")
     const k_max = m
   end

   const s=2*sh
   #rhs1=LE'*(UA\(LA\(LE*rhs)))

   # Orthogonalize [B,A^{-1}B] with an economy-size QR
   rr = [ rhs LE'*(UA\(LA\(LE*rhs))) ]

   # Julia qr decomposition is always "economy size"
   U,beta=qr(rr)
   U = U[1:n,1:s]

   const ibeta=inv(beta[1:s,1:s])
   beta = beta[1:sh,1:sh]
   const beta2=beta*beta'
   H=zeros((m+1)*s,m*s)
   T=zeros((m+1)*s,m*s)
   L=zeros((m+1)*s,m*s)
   local js, j, rho
   for j=1:m
     jms=(j-1)*s+1
     j1s=(j+1)*s
     js=j*s
     js1=js+1
     jsh=(j-1)*s+sh

     # Expand the basis
     # multiply by A
     Up = zeros(n,s)
     Up[1:n,1:sh] = LE\(A*(LE'\U[:,jms:jsh]))
     # solve with A

     Up[1:n,sh+1:s] = LE'*(UA\(LA\(LE*U[1:n,jsh+1:js])))

     # orthogonalize the new basis block wrt all the previous ones by modified gram
     for l=1:2
        k_min=max(1,j-k_max)
        for kk=k_min:j
            k1=(kk-1)*s+1
            k2=kk*s
            coef= U[1:n,k1:k2]'*Up
            H[k1:k2,jms:js] = H[k1:k2,jms:js]+ coef
            Up = Up - U[:,k1:k2]*coef
        end
      end

      if j<=m
        Up,H[js1:j1s,jms:js] = qr(Up)
        hinv=inv(H[js1:j1s,jms:js])
      end

      ###############################################################
      # Recover the columns of T=U'*A*U (projection of A onto the space) from
      # the colums of H.
      # REMARK: we need T as coefficient matrix of the projected problem.
      Iden=speye(js+s)

      if (j==1)
        L[1:j*s+sh,(j-1)*sh+1:j*sh] = [H[1:s+sh,1:sh]/ibeta[1:sh,1:sh] eye(s+sh,sh)/ibeta[1:sh,1:sh]]*ibeta[1:s,sh+1:s]
      else
        L[1:j*s+s,(j-1)*sh+1:j*sh] = L[1:j*s+s,(j-1)*sh+1:j*sh] + H[1:j*s+s,jms:jms-1+sh]*rho
      end

      odds = [odds; jms:(jms-1+sh)]   # store the odd block columns
      evens = 1:js
      flag = trues(size(evens))
      flag[odds] = false
      evens = evens[flag]
      T[1:js+s,odds]=H[1:js+s,odds]   #odd columns

      T[1:js+sh,evens]=L[1:js+sh,1:j*sh]   #even columns
      L[1:j*s+s,j*sh+1:(j+1)*sh] = ( Iden[1:j*s+s,(js-sh+1):js]- T[1:js+s,1:js]*H[1:js,js-sh+1:js])*hinv[sh+1:s,sh+1:s]
      rho = hinv[1:sh,1:sh]\hinv[1:sh,sh+1:s]

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

      # Solve the projected problem by Bartels-Stewart
      # Do "type lyap" from command window if interested

      Y = lyap((T[1:js,1:js]),eye(j*s,sh)*beta2*eye(j*s,sh)')

      # safeguard to preserve symmetry
      Y = (Y+Y')/2

      # Compute the residual norm. See the article by Valeria

      cc = [H[js1:j1s,js-s+1:js-sh] L[js1:j1s,(j-1)*sh+1:j*sh]]

      nrmx = vecnorm(Y)

      er2=[er2;sqrt2*vecnorm(cc*Y[js-s+1:js,:])/(nrmb+singE*nrma*nrmx)]

      infoV && println("KPIK It: $j -- Current Backwards Error: $(er2[j])")

      (er2[j]<tol) ? break : U = [U Up]

    end

    # Done
    # reduce solution rank if needed
    sY,uY=eig(Y)
    id=sortperm(sY)
    sY=sort(sY)

    sY=flipdim(sY,1)
    uY=uY[:,id[end:-1:1]]
    is = 0
    for ii in 1:size(sY)[1]
      if abs(sY[ii])>tolY
        is = is+1
      end
    end

    Y0 = uY[:,1:is]*diagm(sqrt(sY[1:is]))
    Z = LE'\(U[1:n,1:js]*Y0)
    er2=er2[1:j]

    if infoV
      println("its  Back. Error            space dim. CPU Time")
      println("$j    $(er2[j])  $js          $(toq())")
    end

    return Z, er2

end

"""
  `rksm(A,B,E=1;LE=1,s1=NaN,emax=NaN,m=100,tol=1e-9,ch=true,tolY=1e-12,infoV=true)`

Based on rksm.m on Valeria Simoncini's website
(http://www.dm.unibo.it/~simoncin/software.html)
and includes many of the same comments and much of the same description

Approximately Solve

    A X E + E X A' + BB' = 0

by the Rational Krylov subspace method
(Galerkin condition onto the Rational Krylov subspace)
This code performs system solves with (A-s E)

Input:

    A, E  coeff. matrices. A<0,  E is spd

    LE   Cholesky lower factor of E

    B     rhs factor

    m       max space dimension allowed

    tol     stopping tolerance based on the backwards error,
            with stopping criterion

                ||LE\A X LE  + LE' X A'/LE'-LE\BB'/LE'||
            ----------------------------------------------  < tol
            ||LE\BB'/LE'|| + ||E^{-1}|| ||A|| ||LE'X LE ||

            computed in a cheap manner

    s1,smax   estimates for real spectral interval associated with
              field of values of (A,E)

    ch      ch=true  complex poles  ch=false real poles

    tolY    truncation tolerance for final solution, e.g., tolY=1e-12

    infoV   If true, print out backwards error every iteration and a summary
            once the program finishes.

Output:

    Z    factor of approximate solution  X = Z Z'
         If ch=true,  Z may be complex

    nrmrestot  history of residuals

Hints:

    1. Before the call permute entries of A, E and B so as to
    limit fill-in in the system solves

    2. Provide "comfortable" (loose bounds) estimates s1, emax

When using this code, please cite the following reference:

    V. Druskin and V. Simoncini,
    Adaptive rational Krylov subspaces for large-scale dynamical systems
    Systems & Control Letters, 60 (2011), pp. 546-560.

#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
#FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
#COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
#IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
#CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
function rksm(A,B,E=1;LE=1,s1=NaN,emax=NaN,m=100,tol=1e-9,ch=false,tolY=1e-12,
                                                          infoV=true)

@assert(isdefined(:vecnorm),"Your julia version is too old. vecnorm() not defined")

infoV && tic()

const symm = norm(A-E-(A-E)',1) < 1e-14
# If not symmetric, we may have complex poles.
# Therefore we must use complex types
symm ? (const typ = Float64) : (const typ = Complex{Float64})

const n=size(A,1)
const p=size(B,2)

E==1 ? ee = eye(n) : ee = E
# If user does not give smallest generalized eigenvalue, calculate it
if (isnan(s1))
  (symm) && (s1 = eigmax(inv(full(ee))*full(-A)))
  (!symm) && (s1 = maximum(abs(eigvals(inv(full(ee))*full(-A)))))
end

# If user does not give largest generalized eigenvalue, calculate it
if (isnan(emax))
  (symm) && (emax = eigmin(inv(full(ee))*full(-A)))
  (!symm) && (emax = minimum(abs(eigvals(inv(full(ee))*full(-A)))))
end

if (E != 1 && LE == 1)
  if issparse(E)
    LE = try
          LE = sparse(cholfact(E,perm=1:size(E,1))[:L])
        catch
          error("E must be Symmetric Positive Definite")
        end
  else
    LE = try
          cholfact(E)[:L]
        catch
          error("E must be Symmetric Positive Definite")
        end
 end
end

const Iden=eye(p)
const O=0*Iden
const uno=ones(1,p)
const Lres = convert(Array{typ,2},LE\B)

V,irr=qr(Lres)
rr=inv(irr)
const nrmb=vecnorm(inv(rr))^2
const beta=V'*Lres
const beta2=beta*beta'
s=typ[]

VV=zeros(n,p*(m+2))
(typ == Float64) && (VV=VV*1.)
(typ == Complex{Float64}) && (VV=VV*1.*im)
VV[1:n,1:p]=V

H=zeros(p*(m+2),p*(m+1))
(typ == Float64) && (H=H*1.)
(typ == Complex{Float64}) && (H=H*1.*im)

nrmrestot=[]
const nrma=vecnorm(A)

if E == 1
  const condestE=1
  const singE=condestE/vecnorm(eye(n))
elseif (vecnorm(E-speye(n))>1e-14)
  const condestE=cond(full(E))
  const singE=condestE/vecnorm(full(E))
else
  const singE=1
end

newAv=LE\(A*(LE'\V))

K=V'*newAv
push!(s,s1)
eH=eig(K)
eHpoints = sort([s1,emax])
snew=newpolei(eHpoints,eH[1],s1*uno',typ)

if real(snew)<0
   snew=-real(snew)+im*imag(snew)
end
push!(s,snew)

# additional steps
cmplxflag=false

i=0

local Y, j1s
while i < m
  local i1, j1, js
  i=i+1

  paired=0
  while (paired==0)

    i1=i+1
    w=LE*V
    #E==1 ? wrk = (A-snew*diagm(ones(n),0))\w : wrk = (A-snew*E)\w
    wrk = (A-snew*speye(n)*E)\w
    wrk= LE'*wrk

    # Gram-Schmidt step
    jms=(i-1)*p+1
    j1s=(i+1)*p
    js=i*p
    js1=js+1
    for it=1:2
      for kk=1:i
        k1=(kk-1)*p+1
        k2=kk*p
        gamma=VV[1:n,k1:k2]'*wrk
        H[k1:k2,jms:js] = H[k1:k2,jms:js]+ gamma
        wrk = wrk - VV[:,k1:k2]*gamma
      end
    end
    V, H[js1:j1s,jms:js]=qr(wrk)
    if (cmplxflag)
      snew=conj(snew)
      push!(s,snew)
      cmplxflag=false
      newAv=LE\(A*(LE'\V))
      g = VV[1:n,1:js]'*newAv
      g1 = g
      g2 = V'*(LE\(A*(LE'\VV[1:n,1:js])))
      g3 = V'*(LE\(A*(LE'\V)))
      K = [K g1; g2 g3]
      VV[1:n,js+1:j1s]=V
      i=i+1
    else
      paired=1
    end
  end

  ih1=i1
  ih=i
  newAv=LE\(A*(LE'\V))
  g = VV[1:n,1:js]'*newAv

  if (symm)
    K=(K+K')/2
  end

  rhs2=speye(ih*p,p)*beta2*speye(ih*p,p)'
  Y = lyap(K,rhs2)
  nrmx = vecnorm(Y)

  # computed residual   (exact, in exact arithmetic)
  u1=newAv-VV[1:n,1:js]*g
  d=-VV[1:n,1:js]*(Y*(H[1:ih*p,1:ih*p]'\[zeros(p*(ih-1),p);Iden])*H[p*ih+1:p*ih1,p*ih-p+1:p*ih]')

  U=[-V*s[end]  d u1 ]
  extra,rr=qr(full(U))

  nrmres=vecnorm(rr*([O Iden O; Iden O Iden; O Iden O ])*rr')/(nrmb+singE*nrma*nrmx)
  push!(nrmrestot, nrmres)

  infoV && println("RKSM It: $i -- Current backwards error: $nrmres")
  (nrmres<tol) && break
  # New poles and zeros
  eH=eig(K)[1]
  symm ? sort!(eH) : eH[sortperm(abs(eH))]

  eHorig=eH

   if (ch)
     # Complex poles. Compute set for next complex pole of r_m

      if (countnz(imag(eH))>0 && length(eH)>2) # Roots lambdas come from convex hull too
        eH=full([eH;-emax])
        eH=convhull(eH)
        ieH=length(eH)
        missing=ih*p-ieH
        while missing>0                        # include enough points from the border
          neweH=(eH[1:ieH-1]+eH[2:ieH])/2
          missing=ih*p-length(eH)
          eH=[eH;neweH]
        end
        eHpoints=-eH
        eH=eHorig
      else     # if all real eigs, no convex hull possible
        eHpoints = sort([s1; emax.';-real(eH)])
      end

   else   # Real poles s from real set. Compute complex roots of r_m via Ritz convex hull
       if (countnz(imag(eH))>0 && length(eH)>2)    # Roots lambdas come from convex hull too
         eH=full([eH;-s1;-emax.'])
         eH=convhull(eH)
         ieH=length(eH)
         missing=ih*p-ieH
         while missing>0 # include enough points from the border
           neweH=(eH[1:ieH-1]+eH[2:ieH])/2
           eH=[eH;neweH]
           missing=ih*p-length(eH)
         end
         eH=eH[1:ih*p]
       end
        eHpoints = sort([s1; emax.';-real(eH)])
        eH=eHorig
   end

   gs=kron(s[2:i+1].',uno)'

   snew = newpolei(eHpoints,eH,gs,typ)
   (real(snew)<0) && (snew=-real(snew)+im*imag(snew)); #safeguard strategy

   # If pole is complex, include its conjugate
   (imag(snew) !=0) && (cmplxflag=true)

   push!(s,snew)

   g1 = g
   g2 = V'*(LE\(A*(LE'\VV[1:n,1:js])))
   g3 = V'*(LE\(A*(LE'\V)))
   K = [K g1; g2 g3]
   VV[1:n,js+1:j1s]=V

end

# Done
# Reduce rank of solution, if needed
sY,uY=eig(Y)
symm ? id = sortperm(sY) : id = sortperm(abs(sY))

sY=sY[id]
sY=flipdim(sY,1)
uY=uY[:,id[end:-1:1]]
is = 0
for ii in 1:size(sY)[1]
  if abs(sY[ii])>tolY
    is = is+1
  end
end

Y0 = uY[:,1:is]*diagm(sqrt(sY[1:is]))
Z = LE'\(VV[:,1:size(Y0,1)]*Y0)

if infoV
  RKStotal_time=toq()
  println("Space dim $j1s  Solution rank $is  time $RKStotal_time")
end

return Z,nrmrestot

end

#################################################
#
#  ratfun and newpolei function definitions. These
#  are auxiliary routines for rksm based on Simoncini's
#  corresponding routines.
#
##################################################

function ratfun(x,eH,s,typ)
r = Array(typ,length(x))

for j=1:length(x)
  xj = x[j]*ones(length(s))
  r[j]=abs(prod( (xj-s)./(xj-eH) ))
end
return r
end

##################################
function newpolei(eHpoints,eH,s,typ)
snew=Array(typ,length(eHpoints)-1)

for j=1:length(eHpoints)-1

    sval = linspace(real(eHpoints[j]),real(eHpoints[j+1]),200)
    (typ == Complex{Float64}) && (sval = sval + im*linspace(imag(eHpoints[j]),imag(eHpoints[j+1]),200))

    jx = indmax(abs(ratfun(sval,eH,s,typ)))

    snew[j]=sval[jx]
end
jx=indmax(abs(ratfun(snew,eH,s,typ)))
snew=snew[jx]
return snew
end

"""
  `convexhull(pnts)`
Based on code from (https://github.com/intdxdt/convexhull.jl)
Function to compute 2-D convex hull
"""
function convhull(pnts)
    const T = eltype(pnts) #get point type
    const N = length(pnts) #number of pnts
    # sort the points lexicographically.
    # copy points into mutable container

    #trivial case 0 or 1 point
    length(pnts) <= 1 && (return pnts)

    #sort points lexicographically
    sort!(pnts, lt=lt2d)

    #function to orient boundry using 2D cross product of OA and OB vectors,
    #i.e. z-component of their 3D cross product.
    #Returns a positive value, if OAB makes a counter-clockwise turn,
    #negative for clockwise turn, and zero if the points are collinear.
    orient(b, pnt) = (real(b[end]-b[end-1]) * imag(pnt - b[end-1])) -
                      (imag(b[end-1] - b[end-1]) * real(pnt - b[end-1]))

    #build lower hull
    lower = T[]
    for i in 1:N
        pnt = pnts[i]
        while length(lower) >= 2 && orient(lower, pnt) <= 0
            pop!(lower)
        end
        push!(lower, pnt)
    end

    #build upper hull
    upper = T[]
    for i in N:-1:1
        pnt = pnts[i]
        while length(upper) >= 2 && orient(upper, pnt) <= 0
            pop!(upper)
        end
        push!(upper, pnt)
    end

    # omit end points :
    #   - reapeated at the beginning of lower & upper
    length(upper) > 0 && pop!(upper)
    length(lower) > 0 && pop!(lower)

    # two reapeated points
    if (length(upper)==1 && length(lower)==1 && isequal(upper, lower))
      pop!(upper)
    end

    # concat lower and upper hull.
    hull =[lower; upper]

    #close ring
    length(hull) > 1 && push!(hull, hull[1])

    return hull
end

#less than comparator
function lt2d(a, b)
    dx, dy =  real(a - b), imag(a - b)
    dx != 0 &&  (return <(dx, 0))
    return <(dy, 0)
end

############################################################
#
#     Please note: The following functions lp_lradi, lp_nrmu, lp_para,
#                  lp_arn, lp_mnmx are based on the functions
#                  with the same name in LYAPACK 1.0 (lp_arn
#                  is based on lp_arn_m.m and lp_arn_p.m). The
#                  structure of these files is closely followed
#                  and even includes many of the same comments.
#
#############################################################

"""
Based on lp_lradi.m from LYAPACK 1.0

Low Rank Cholesky Factor Alternating Direction Implicit (LRCF-ADI) Method for
solving the stable Lyapunov equation:

  For tp = :B

    F * X + X * F' = -B * B',

  for tp = :C

    F' * X + X * F = -B' * B,

where F = A-Bf*Kf'.

The routine works in two modes (depending on the choice of zk)

    For zk = :Z, this routine deliveres a low rank factor Z,
    such that Z * Z' approximates X

    for zk = :K, it only generates the product K_out = Z * Z'* K
    without forming Z itself.

Input:

    A         matrix A
    B         matrix B (n-x-m if tp=:B or m-x-n if tp=:C,
              where  m << n !)
    p         l-vector with ADI shift parameters. Complex parameters must
              appear as conjugate complex pairs in consecutive order.
              If more than length(p) iterations need to be done, the
              parameters p[i] are reused in a cyclic manner. The user can calculate
              these values using lp_para.
    Bf        "feedback" matrix Bf
              Named argument with defult value [], so if not provided by user,
              the code solves
              A * X + X * A' = -B * B' or A' * X + X * A = -B'* B
              depending on choice of tp.
    Kf        "feedback" matrix Kf
              amed argument with defult value [], so if not provided by user,
              the code solves
              A * X + X * A' = -B * B' or A' * X + X * A = -B' * B
              depending on choice of tp.
    tp        (= :B or :C) named argument that determines the type of Lyapunov equation. Defult value is :B
    rc        (= :R or :C) named argument that determines whether the low rank factor Z
              must be real (:R) or if complex factors are allowed (:C)
              If p contains complex parameters, then the low rank factor
              Z is complex, too, although Z * Z' is real. A real low rank
              factor is determined from the complex data, if rc = :R.
              However, this requires some additional computation.
              If zk = :K, rc is ignored. Defult value is :C.
    K         named argument n-x-r matrix (where r should be small: r << n !).
              Defult value is [].
    max_it    stopping criterion: maximal number of LRCP-ADI steps
              Set to [] or +Inf (default) to avoid this criterion. Named argument
              With defult value of 100.
    min_res   stopping criterion: minimal relative residual norm. The
              iteration is stopped when res(i+1) <= min_res (See Remarks).
              Set min_res = [] or 0 (default) to avoid this criterion. Note:
              If min_res<=0 and with_rs=:N, the (often expensive) calculation
              of the residual norm is avoided, but, of course, res is not
              provided on exit. Named argument with defult value of 1e-6.
    with_rs   (= :S or :N) if with_rs = "S", the iteration is stopped,
              when the routine detects a stagnation of the residual norms,
              which is most likely the case, when roundoff errors rather
              than the approximation error start to dominante the residual
              norm. This implies that the residual norms are computed (which
              can be expensive). This criterion works quite well in practice,
              but is not absolutely sure. Use with_rs = "S" only if you
              want to compute the Lyapunov solution as accurate as possible
              (for a given machine precision).
              If with_rs = :N, this criterion is not used. Named argument with
              defult value :S.
    min_in    stopping criterion: This value limits the "minimal increase"
              in the matrix Z by the "new" columns . The iteration is
              terminated as soon as

                || Z_nc ||F^2 < min_in * || Z_new ||F^2

              holds for a certain number of consecutive iteration steps.
              Here, Z_nc are the currently computed "new" colums, which
              appended to the old iterate Z_old deliver the new iterate
              Z_new  = [ Z_old  Z_nc ]. Set min_in = 0 to avoid it.
              Default value is eps, the machine precision. min_in = []
              has the same effect. Note that this is an "adaptive"
              stopping criterion which does not require the
              (often expensive) computation of the residual norm. Named argument
              with defult value 0.
    infoV      Bool: Return infoVrmation given during the
              iteration. Named argument with default value true.

Output:

    Z         Z * Z' is a low rank approximation of X
              (Note that Z can be complex if rc=:C!)
    flag      the criterion which had stopped the iteration:
               = "I": max_it,
               = "R": min_res,
               = "S": with_rs,
               = "N": min_in,
    res       the relative residual norms attained in the course of
              the iterations (res[i+1] is the norm after the i-th step
              of the iteration!). See note in min_res.

Remarks:

    1. The eigenvalues of F must have negative real parts.

    2. The values in res correspond to the following "relative" norms

      tp = :B

        res(i+1) = ||F * Z_i * Z_i' + Z_i * Z_i' * F' + B * B'||F/|| B * B'||F

      tp = :C

        res(i+1) = ||F' * Z_i * Z_i' + Z_i * Z_i' * F + B' * B||F/||B' * B||F

    3. Note that all stopping criteria are checked only after a step
    with a real parameter or a "double step" with a pair of conjugate
    complex parameters. This ensures that Z*Z' is real, even if Z is
    not.

References:

    The algorithm is a slight modifivation of that proposed in:

    1. J.Li, F.Wang, and J.White.
    An efficient Lyapunov equation based approach for generating
    reduced-order models of interconnect.
    Proceedings of the 36th IEEE/ACM Design Automation Conference,
    New Orleans, LA, 1999.

    Another (though more expensive) low rank algorithm is proposed in:

    2. T. Penzl.
    A cyclic low rank Smith method for large sparse Lyapunov equations.
    To appear in SIAM Sci. Comp.

    See also:

    3. P. Benner, J. Li, and T. Penzl
    Numerical solution of large Lyapunov equations, Riccati equations,
    and linear-quadratic optimal control problems.
    In preparation.

    4. T. Penzl.
    LYAPACK (Users' Guide - Version 1.0).
    1999.

Internal remarks:

    Input data not completely checked!

    The procedure to generate real factors in case of complex parameters is
    different from that in [3]!

    The matrices SMi (i = 1:length(p)) for the "Sherman-Morrison trick"
    (only used if Bf and Kf are nonzero) are computed a priori, which
    is good in view of computation if parameters p(i) are used cyclically,
    but may be sometimes bad in view of memory demand, in particular, when
    length(p) is large.

    The stopping criterion related to the input parameter with_rs
    corresponds to the stagnation of the residual curve caused by
    round-off errors. Its performance depends on the constants stcf and
    min_rs. The iteration is stopped, as soon as

    (ra-rb) * (i-stcf+1) / ((r(1)-ra) * stcf) < min_rs

    and r(1)-ra > 0 and i >= stcf hold for stcf consecutive iteration
    steps. Here res(i+1) is the residual norm after the i-th LRCF-ADI step,
    r(i+1) = log(res(i+1)), ra = min(r(1:i-stcf+1)), rb = min(r(i-stcf+2:i+1)).

    stcf is also the number of consecutive steps, for which the criterion
    w.r.t. min_in must be fulfilled.
"""
function lp_lradi(A,B;p=[],Bf=[],Kf=[],K=[],max_it=100,tp=:B,rc=:C,
  min_res=1e-6,with_rs=:S,min_in=0,infoV=true)

tic()
isempty(p) && (p=lp_para(A))
const stcf = 10
const min_rs = .1

(tp!=:B && tp!=:C) && error("tp must be either :B or :C.")
(rc!=:R && rc!=:C) && error("rc must be either :R or :C.")

const tpB = (tp == :B)

const with_min_rs = (with_rs==:S)
const with_norm = (min_res>0)||with_min_rs
const make_real = (rc==:R)

const with_min_in = min_in>0

const with_BK = !isempty(Bf)

const l = length(p)

if tpB
  const n,m = size(B)
else
  const m,n = size(B)
end

const Ide = speye(size(A,1))
LP_L = Array(Complex{Float64},size(A)...,l)
LP_U = Array(Complex{Float64},size(A)...,l)
for i = 1:l
  LP_L[:,:,i],LP_U[:,:,i],~=lu(full(A+p[i]*Ide),Val{false})
end

if with_BK
  SM=Array(Complex{Float64},size(Bf)...,l)
  Im = speye(size(Bf,2))
  if tpB
    # SMi = TM*inv(I-Kf'*TM),
    # where TM = inv(F+p(i)*I)*Bf
    # (These are the columns of the LR terms for the
    # rank corrections of the "inverse"
    # in the Sherman-Morrison formulae.)

    for i = 1:l
      TM = LP_U[:,:,i]\(LP_L[:,:,i]\Bf)
      SM[:,:,i] = TM/(Im-Kf'*TM)
    end
  else  # tp==:C
    # SMi = TM*inv(I-Bf'*TM),
    # where TM = inv(F.'+p(i)*I)*Kf
    # (These are the columns of the LR terms for the
    # rank corrections of the "inverse"
    # in the Sherman-Morrison formulae.)
    for i = 1:l
      TM = LP_L[:,:,i].'\(LP_U[:,:,i].'\Kf)
      SM[:,:,i] = TM/(Im-Bf'*TM)
    end
  end
end

# Initialize QR factorization for norm computation
if with_norm
  res0,nrmQ,nrmR,nrmbs = lp_nrmu(A,B,Bf,Kf,tpB,[],[],[],[])
  res = [1.]
  res_log = [log(res0)]
  # Vector containing log of residual norms
  # corresponds to r(:) in prolog.
end

flag = "I"
if with_min_in
  nrmF_Z_2 = 0
  nrmF_rec = +Inf*ones(stcf,1)
  # Current squared Frobenius norm of Z
  # Records the values of
  # ||Z_nc||_F^2 / ||Z_new||_F^2 (see prolog)
  # for the last stcf iteration steps.
end

i_p = 1;                   # Pointer to i-th entry in p(:)
is_compl = imag(p[1])!=0;  # is_compl = (current parameter is complex.)
is_first = true;           # is_first = (current parameter is the first
                           #            of a pair, PROVIDED THAT is_compl.)

local V, Z
for i = 1:max_it       # The iteration itself

  if i==1

    if tpB
      # V = last columns of Cholesky factor Z
      # Initialize:
      # V := sqrt(-2*real(p[1]))*inv(F+p[1]*I)*B

      V=LP_U[:,:,1]\(LP_L[:,:,1]\B)
      with_BK && (V = V+SM[:,:,1]*(Kf'*V))
      V = sqrt(-2*real(p[1]))*V

    else #( tp = :C )
      # Initialize:
      # V := sqrt(-2*real(p[1]))*inv(F.'+p[1]*I)*B'

      V = LP_L[:,:,1].'\(LP_U[:,:,1].'\B')
      with_BK && (V = V+SM[:,:,1]*(Bf'*V))
      V = sqrt(-2*real(p[1]))*V
    end

    Z = V;             # Note:  Z*Z' = current ADI iterate

  else # i > 1

    p_old = p[i_p]

    i_p = i_p+1
    i_p>l && (i_p = 1)    # update current parameter index

    if is_compl && is_first
      is_first = false
      if i_p==1
        error("Parameter sequence ends in the 'middle' of a complex pair!")
      end
      if p[i_p]!=conj(p_old)
        error("Parameters p[i] must be either real or pairs of conjugate complex numbers.")
      end
    else
      is_compl = imag(p[i_p])!=0
      is_first = true
    end

    if tpB
      # Evaluate
      #   V := sqrt(real(p[i_p])/real(p_old))*...
      #        (V-(p[i_p]+conj(p_old))*inv(F+p[i_p]*I)*V)

      TM = LP_U[:,:,i_p]\(LP_L[:,:,i_p]\V)
      with_BK && (TM = TM+SM[:,:,i_p]*(Kf'*TM))
      TM = V-(p[i_p]+conj(p_old))*TM
      V = sqrt(real(p[i_p])/real(p_old))*TM

    else   #  tp==:C

      # Evaluate
      #   V := sqrt(real(p[i_p])/real(p_old))*...
      #        (V-(p[i_p]+conj(p_old))* inv(F.'+p[i_p]*I)*V)

      TM = LP_L[:,:,i_p].'\(LP_U[:,:,i_p].'\V)
      with_BK && (TM = TM+SM[:,:,i_p]*(Bf'*TM))
      TM = V-(p[i_p]+conj(p_old))*TM
      V = sqrt(real(p[i_p])/real(p_old))*TM
    end

    # Form new iterate Z.

    !is_compl && (V = real(V))

    Z = [Z V]

    # Make last 2*m columns real.
    if make_real && is_compl && !is_first
      for j = (i-1)*m+1:i*m
        U1,S1,V1 = svd([real(Z[:,j-m]) real(Z[:,j]) imag(Z[:,j-m]) imag(Z[:,j])])
        S1 = diagm(S1)
        U2,S2,V2 = svd(V1[1:2,1:2]'*V1[1:2,1:2]+V1[3:4,1:2]'*V1[3:4,1:2],thin=false)
        TMP = U1[:,1:2]*S1[1:2,1:2]*U2*diagm(sqrt(S2))
        Z[:,j-m] = TMP[:,1]
        Z[:,j] = TMP[:,2]
      end

    end

  end

  if with_norm
    # Compute residual norm

    resnrm,nrmQ,nrmR,nrmbs = lp_nrmu( A, B, Bf, Kf, tpB, V, nrmQ, nrmR, nrmbs)
    push!(res_log, log(resnrm))
    akt_res = resnrm/res0
    push!(res, akt_res)

    infoV && println("LRCF-ADI step $i -- norm. residual = $akt_res")

    # After pair of complex parameters or
    # a real parameter, check stopping criteria
    # based on residual norm.
    if  !(is_compl && is_first)

      if akt_res <= min_res
        flag = "R"
        break
      end

      if with_min_rs && i>=2*stcf
        ra = minimum(res_log[1:i-stcf+1])
        rb = minimum(res_log[i-stcf+2:i+1])
        if res_log[1]-ra > 0 && (ra-rb)*(i-stcf+1)/((res_log[1]-ra)*stcf) < min_rs
          flag = "S"
          break
        end
      end

    end

  end
  # Check stopping criteria based on increase in ||Z_i||_F.
  if with_min_in
    nrmF_V_2 = sum(sum(abs(V).*abs(V)));    # Note the "abs"; V is complex.
    nrmF_Z_2 = nrmF_Z_2 + nrmF_V_2
    nrmF_rec[1:stcf-1] = nrmF_rec[2:stcf]
    nrmF_rec[stcf] = nrmF_V_2/nrmF_Z_2
    if !(is_compl && is_first) && i>stcf && all(nrmF_rec .< min_in)
      flag = "N"
      break
    end
  end
end

infoV && println("Computational Time:  ", toq())

Z, flag, res
end

"""
Based on lp_nrmu.m from LYAPACK 1.0

Helper function for lp_lradi

Using updated QR factorizations, this routine computes efficiently
a sequence of norms which correspond to either of the following types:

for tp = :B

    nrm = || F*Z_i*Z_i' + Z_i*Z_i'* F' + B*B' ||F,

for tp = :C

    nrm = || F' * Z_i * Z_i' + Z_i * Z_i' * F + B' * B ||F.

Here, F = A-Bf*Kf'.

The matrices Z_i must have much more rows than columns and they
have to obey the recursion

    Z_i = [ Z_{i-1}  V ].

Calling sequence:

    nrm, nrmQ, nrmR, nrmbs = ...
    lp_nrmu(  A, B, Bf, Kf, tpB, V, nrmQ, nrmR, nrmbs )

Input:

    tpB        (= true or false) the type of the norm

    Bf        real matrix Bf
              Set Bf = [] if not existing or zero!

    Kf        real matrix Kf
              Set Kf = [] if not existing or zero!

    B         n-x-m or m-x-n matrix B (must be real)

    V         n-x-r matrix V (may be complex)

    nrmQ,

    nrmR,

    nrmbk     variables for internal use (they contain the data of
              the updated QR factorization).

Output:

    nrm       the current value of the Frobenius norm nrm_i.

    nrmQ,

    nrmR,

    nrmbk     variables for internal use (they must be output
              parameters).

Remarks:

    Using this routine for computing the norms within an iteration
    can save a lot of computation compared to repeatedly calling
    'lp_nrm'.

    This routine must be used as follows:

    1. Before the iteration starts, it must be invoked with the
       parameters V = nrmQ = nrmR = nrmbl = []. This initializes the
       variables nrmQ, nrmR, nrmbl.

    2. The routine must be called in each step of the iteration:

Example:

    ...

    nrm,nrmQ,nrmR,nrmbs = lp_nrmu(A,B,[],[],tpB,[],[],[],[])

    ...

    Z = zeros(n,0)

      ...

      for i = 1:100

        ...

        V = ...

        ...

        Z = [ Z V ];    # iteration:  Z_i = [ Z_i-1 V ]

        ...

        nrm,nrmQ,nrmR,nrmbs = lp_nrmu(A,B,[],[],tpB,V,nrmQ,nrmR,nrmbs)

        ...

      end
      ...
"""
function lp_nrmu( A, B, Bf, Kf, tpB, V, nrmQ, nrmR, nrmbs )

const with_BK = length(Bf)>0

const n= size(A,1);                    # Get system order.

if isempty(V)
  # The routine is called for the first time (before the iteration
  # 'outside' is started (i.e. V = zeros(n,0) or []!)) in order to
  # initialize the QR factorization nrmQ*nrmR correctly.

  if tpB
    nrmbs = [size(B,2)]
    nrmQ,nrmR = qr(B)
  else
    nrmbs = [size(B,1)]
    nrmQ,nrmR = qr(B')
  end
  nrm = vecnorm(nrmR*nrmR')

else
  # The routine is called during the iteration 'outside'.
  # The QR factorization is updated and the norm is computed.

  # Update of the QR factorization.
  push!(nrmbs,size(V,2))
  lw = size(nrmQ,2)
  lz = size(V,2)

  if tpB
    TM = A*V
    with_BK && (TM = TM-Bf*(Kf'*V))
    Z = [ TM V ]
  else
    TM = A'*V
    with_BK && (TM = TM-Kf*(Bf'*V))
    Z = [ TM V ]
  end

  for j = 1:2*lz
    a = Array(Complex{Float64},lw,1)
    #a = complex(zeros(lw,1))
    t = Z[:,j]
    for k = 1:lw
      u = nrmQ[:,k]
      alpha = dot(u,t)
      t = t-alpha*u
      a[k] = alpha[1]
    end
    beta = norm(t)
    nrmQ = [nrmQ t/beta]
    #nrmQ = push!(nrmQ, t/beta...)
    nrmR = [ nrmR a ; zeros(1,lw) beta ]
    lw = lw+1
  end

  # Computation of  nrmR * [permutation matrix] * nrmR'
  RT = copy(nrmR)
  ie2 = nrmbs[1]
  for j = 2:length(nrmbs)
    ia1 = ie2+1
    ie1 = ie2+nrmbs[j]
    ia2 = ie1+1
    ie2 = ie1+nrmbs[j]
    TMP = RT[:,ia1:ie1]
    RT[:,ia1:ie1] = RT[:,ia2:ie2]
    RT[:,ia2:ie2] = TMP
  end
  nrm = vecnorm(nrmR*RT')
end
nrm, nrmQ, nrmR, nrmbs
end

"""
    `lp_para(A,Bf,Kf,l0,kp,km,b0)`

Based on lp_para.m from LYAPACK 1.0

Estimation of suboptimal ADI shift parameters for the matrix

F = A - Bf * Kf'.

Calling sequence:

p = lp_para(A,Bf,Kf,l0,kp,km,b0)

Input:

    Bf        matrix Bf
              Set Bf = [] if not existing or zero!
    Kf        matrix Kf
              Set Kf = [] if not existing or zero!
    l0        desired number of shift parameters (kp+km > 2*l0)
              (The algorithm delivers either l0 or l0+1 parameters!)
    kp, km    numbers of Arnoldi steps w.r.t. F and inv(F),
              respectively (kp, km << order of A)
    b0         Arnoldi start vector (optional; chosen at random if not
              provided).

Output:

    p         an l0- or l0+1-vector of suboptimal ADI parameters

Remarks:

    Typical values are l0 = 10..40, kp = 20..80, km = 10..40.
    The harder the problem is the large values are necessary.
    Larger values mostly result in a faster convergence, but also in a
    larger memory requirement.
    However, for "well-conditioned" problems small values of l0 can
    lead to the optimal performance.

References:

    1. T. Penzl.
    LYAPACK (Users' Guide - Version 1.0).
    1999.

Input data not completely checked!
"""
function lp_para(A,Bf=[],Kf=[],l0=NaN,kp=NaN,km=NaN,b0=[])

err_code = 0
n=size(A,1); # Get system order.

if isnan(l0) && isnan(kp) && isnan(km)
  if n > 51
    kp=50;
    km=25;
    l0=20;
  elseif n > 15
    kp = n - 5
    km = Int(round(kp/2))
    l0 = km-2
  else
    kp = n-2
    km = Int(round(kp/2))
    l0 = km-1
  end
else
  kp >= n && error("kp must be smaller than n!")
  km >= n && error("km must be smaller than n!")
  2*l0 >= kp+km && error("2*l0 must be smaller than kp+km!")
end

isempty(b0) && (b0 = randn(n,1))
b0 = (1/norm(b0))*b0

rwp = []
rwm = []
rw = []
Hp = []
Hm = []

if kp > 0
  Hp,V = lp_arn(A,Bf,Kf,kp,b0,:p)
  rwp = eigvals(Hp[1:kp,1:kp]);                  # =: R_+
  rw = rwp
end

if km > 0
  Hm,V = lp_arn(A,Bf,Kf,km,b0,:m)
  rwm = ones(km,1)./eigvals(Hm[1:km,1:km]);      # =: 1 / R_-
  push!(rw,rwm...);                           # =: R
end

if any(real(rw) .>= zeros(size(rw)))
  err_code = 1
  println("These are the Ritz values computed by the Arnoldi process w.r.t. F: $rwp
  These are the Ritz values computed by the Arnoldi process w.r.t. inv(F): $rwm

  ####################################################################
  WARNING in 'lp_para': NON-STABLE RITZ VALUES DETECTED!!!

  This is quite a serious problem, that can be caused by
  (i)   non-stable matrices F (Be sure that F is stable. ADI like
        methods only work for stable or antistable problems. If your
        Lyapunov equation is antistable, multiply it by -1.)
  (ii)  matrices F that are stable but have an indefinite symmetric
        part (This is THE weak point of this algorithm. Try to work
        with the 'reduced' Ritz values, i.e., the unstable values are
        simply removed. This is not an elegant measure but it may work.
        However, the convergence of ADI can be poor. This measure is
        taken automatically. Another measure might be to enlarge the
        values of kp or km, and run the program again.
  (iii) matrices F with a negative definite, but ill-conditioned
        symmetric part (This is quite unlikely. The problem is
        caused by round-off errors).

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

  NOTE: The unstable Ritz values will be ignored in the further computation!!! ")
  rw0 = copy(rw)
  rw = []
  for j = 1:length(rw0)
    if real(rw0[j])<0
      push!(rw,rw0[j])
    end
  end
end

p = lp_mnmx(rw,l0)

p
end

"""
Based on lp_arn_m.m and lp_arn_p.m from LYAPACK 1.0 and includes many of the
same comments much of the same description.

helper function for lp_para

Arnoldi method w.r.t. inv(F), where F = A-Bf*Kf'.

Calling sequence:

    H,V = lp_arn_m(A,Bf,Kf,k,r)

Input:

    A         Matrix A

    Bf        matrix Bf
              Set Bf = [] if not existing or zero!

    Kf        matrix Kf

              Set Kf = [] if not existing or zero!

    k         number of Arnoldi steps (usually k << n, where
              n is the order of the system)

    r         initial n-vector
              (optional - chosen by random, if []).

    pm        %TODO write description

Output:

    H         matrix H ((k+1)-x-k matrix, upper Hessenberg)

    V         matrix V (n-x-(k+1) matrix, orthogonal columns).

Method:

    The ("inverse") Arnoldi method produces matrices V and H such that

    V(:,1) in span{r},
    V' * V = eye(k+1),
      inv(F) * V(:,1:k) = V * H.

Remark:

    This implementation does not check for (near-)breakdown!

"""
function lp_arn(A,Bf,Kf,k,r,pm)

const with_BK = !isempty(Bf)

const n = size(A,1)                 # Get system order.
k >= n-1 && error("k must be smaller than the order of A!")
isempty(r) && (r = randn(n,1))

V = zeros(n,k+1)
H = zeros(k+1,k)

V[:,1] = (1.0/norm(r))*r

beta = 0

LP_L,LP_U,~ = lu(full(A),Val{false})

if with_BK && pm==:m
  # SM = inv(F)*Bf*inv(I-Kf'*inv(F)*Bf)
  # (This is the main part of the term needed for the low
  # rank correction in the Sherman-Morrison formula.)
  Im = speye(size(Bf,2))
  TM = LP_U\(LP_L\Bf)
  SM = TM/(Im-Kf'*TM)
end

for j = 1:k

  if j > 1
    H[j,j-1] = beta
    V[:,j] = (1.0/beta)*r
  end

  if (pm == :m)
    w = LP_U\(LP_L\V[:,j])
    with_BK && (w = w + SM*(Kf'*w)); # LR correction by SM formula
  else
    w = A*V[:,j]
    with_BK && (w = w-Bf*(Kf'*V[:,j]))
  end
  r = w

  for i = 1:j
    H[i,j] = dot(V[:,i],w)
    r = r-H[i,j]*V[:,i]
  end

  beta = norm(r)
  H[j+1,j] = beta

end

V[:,k+1] = (1.0/beta)*r

H,V
end

"""
Based on lp_mnmx.m from LYAPACK 1.0 and includes many of the same comments and
much of the same description.

Suboptimal solution of the ADI minimax problem. The delivered parameter
set is closed under complex conjugation.

Calling sequence:

    p = lp_mnmx(rw,l0)

Input:

    rw        a vector containing numbers in the open left half plane, which
              approximate the spectrum of the corresponding matrix, e.g.,
              a set of Ritz values. The set must be closed w.r.t. complex
              conjugation

    l0        desired number of shift parameters (length(rw) >= l0)
              (The algorithm delivers either l0 or l0+1 parameters!).

Output:

    p         an l0- or l0+1-vector of suboptimal ADI parameters

"""
function lp_mnmx(rw,l0)

if length(rw)<l0
  error("length(rw) must be at least l0.")
end

max_rr = +Inf;                       # Choose initial parameter (pair)
p0 = 0
for i = 1:length(rw)
  max_r,~ = lp_s(rw[i],rw)
  if max_r < max_rr
    p0 = rw[i]
    max_rr = max_r
  end
end

imag(p0)!=0 ? (p = [ p0; conj(p0) ]) : (p = [p0])

max_r,i = lp_s(p,rw);         # Choose further parameters.

while size(p,1) < l0

  p0 = rw[i]
  imag(p0)!=0 ? push!(p,p0,conj(p0)) : push!(p,p0)

  max_r,i = lp_s(p,rw)

end

p
end

"""
Based on lp_s.m from LYAPACK 1.0 and includes much of the same description

Helper function for lp_para

Computation of the maximal magnitude of the rational ADI function over
a discrete subset of the left complex half plane.

Calling sequence:

    max_r,ind = lp_s(p,set)

Input:

    p        vector of ADI parameters

    set      vector representing the discrete set.

Output:

    max_r    maximal magnitude of the rational ADI function over set

    ind      index - maximum is attained for set(ind).
"""
function lp_s(p,set)

max_r = -1
ind = 0

for i = 1:length(set)

  x = set[i]

  rr = 1
  for j = 1:length(p)

    rr = rr*abs(p[j]-x)/abs(p[j]+x)

  end

  if rr > max_r

    max_r = rr
    ind = i

  end

end
max_r,ind
end

# end module matrixEqs
end

max_stars_repo_name: garrettthomaskth/LargeMatrixEquations

max_issues_repo_name: garrettthomaskth/LargeMatrixEquations

max_forks_repo_name: garrettthomaskth/LargeMatrixEquations


text:
struct KDTreeIndex{I<:DataFreeTree} <: AbstractUniqueIndex
    I::I
end

function KDTreeIndex(A::AbstractArray)
    dftree = DataFreeTree(KDTree, A)
    return KDTreeIndex(dftree)
end

Base.summary(::KDTreeIndex) = "KDTreeIndex"

"""Answer using outer bbox of the index."""
function Base.in(x, a::AcceleratedArray{<:Any,<:Any,<:Any,<:KDTreeIndex})
    tree = injectdata(a.index.I, a.parent)  # yields a KDTree
    in(x, tree.hyper_rec)
end

"""Inside Hyperectangle check."""
function Base.in(coordinate::Vector{T}, bbox::HyperRectangle{T}) where {T}
    length(coordinate) <= length(bbox.mins) || error("Coordinate has more dimensions than the bounding box.")
    for (i, ax) in enumerate(coordinate)
        bbox.mins[i] <= ax <= bbox.maxes[i] || return false
    end
    true
end

function Base.findfirst(pred::Base.Fix2{typeof(in),HyperSphere{N,T}}, points::AcceleratedArray{<:Any,<:Any,<:Any,<:KDTreeIndex}) where {N} where {T}
    findall(pred, points)[1]  # this will error on empty sets
end

function Base.findall(pred::Base.Fix2{typeof(in),HyperSphere{N,T}}, points::AcceleratedArray{<:Any,<:Any,<:Any,<:KDTreeIndex}) where {N} where {T}
    tree = injectdata(points.index.I, points.parent)  # yields a KDTree
    hypersphere = pred.x
    inrange(tree, hypersphere.center, hypersphere.r, false)
end

function HyperSphere(coords::Vector{T}, radius::T) where {T<:AbstractFloat}
    HyperSphere(SVector{length(coords)}(coords), radius)
end


text:
using Random, LinearAlgebra, Distributions, PDMats, PyPlot, ProgressMeter, Optim, StatsFuns
using HCubature, QuadGK

struct LineDensity

         joint_a_b::AbstractMvNormal
        marginal_b::AbstractMvNormal
    cond_a_given_b::Function
          bsamples::Array{Array{T,1},1} where T<:Real
        numsamples::Int
                 D::Int

end

function LineDensity(; posterior::MvNormal, numsamples=250)

    joint_a_b = posterior # just another name

    D = round(Int, length(joint_a_b)/2)

    marginal_b = MvNormalPartition(joint_a_b, collect(D+1:2*D))

    cond_a_given_b = MvNormalConditional(joint_a_b, collect(1:D),collect(D+1:2D))

    bsamples = [rand(marginal_b) for i=1:numsamples]

    LineDensity(joint_a_b, marginal_b, cond_a_given_b, bsamples, numsamples, D)

end

function Base.show(io::IO, a::LineDensity)

    print(io, "Line density.",
              "\n\t Number of samples ", a.numsamples,
              "\n\t Dimension is  ", a.D)
end

#
# Calculates the integral: ∫ N(t⋅x | μ, Σ) dt
# (see also function loganalyticalintegraloft)
#
function analyticalintegraloft(x, μ, Σ)

    # see https://en.wikipedia.org/wiki/Gaussian_integral#Generalizations

    D = length(μ) ; @assert(size(Σ, 1) == size(Σ, 2) == D == length(x))

    a =  0.5*x'*(Σ\x)

    b =  x'*(Σ\μ)

    c = -0.5*μ'*(Σ\μ)

    return 1/sqrt((2π)^D*det(Σ)) * sqrt(π/a) * exp(b^2/(4*a)+c)

end

#
# Calculates: log ∫ N(t⋅x | μ, Σ) dt
# (see also function analyticalintegraloft)
#
function loganalyticalintegraloft(x, μ, Σ)

    # see https://en.wikipedia.org/wiki/Gaussian_integral#Generalizations

    D = length(μ) ; @assert(size(Σ, 1) == size(Σ, 2) == D == length(x))

    a =  0.5*x'*(Σ\x)

    b =  x'*(Σ\μ)

    c = -0.5*μ'*(Σ\μ)

    return - 0.5*D*log(2π) - 0.5*logdet(Σ) + log(sqrt(π/a)) + (b^2/(4*a)+c)

end

function (linepdf::LineDensity)(x)

    sum(map(b -> integralconddensityatx(x, b, linepdf.cond_a_given_b), linepdf.bsamples))

end

function logl(linepdf::LineDensity, x)

    logsumexp(map(b -> logintegralconddensityatx(x, b, linepdf.cond_a_given_b), linepdf.bsamples))

end

function numintegralconddensityatx(x, b, cond_a_b)

    # This serves as numerical verification for function integralconddensityatx
    aux = (x-b)/norm(x-b)

    f(t) = pdf(cond_a_b(b), t * aux)

    quadgk(f, 0.01, 100)[1]

end

function integralconddensityatx(x, b, cond_a_b)

    aux = (x-b)/norm(x-b)

    g = cond_a_b(b)

    analyticalintegraloft(aux, mean(g), cov(g))

end

function logintegralconddensityatx(x, b, cond_a_b)

    aux = (x-b)/norm(x-b)

    g = cond_a_b(b)

    loganalyticalintegraloft(aux, mean(g), cov(g))

end

function example(lpdf::LineDensity)

    if lpdf.D > 3
        @printf("Visualisation only available for D=2 or D=3\n")
        return nothing
    end

    figure()

    c = zeros(lpdf.D, lpdf.numsamples)

    for (index, b) in enumerate(lpdf.bsamples)

        r =  rand()*9 + 1

        a = rand(lpdf.cond_a_given_b(b))

        c[:,index] = a * r + b

        plotvector(b)

        if lpdf.D == 3
            plot([b[1];b[1]+r*a[1]], [b[2];b[2]+r*a[2]], [b[3];b[3]+r*a[3]], "--g")
        else
            plot([b[1];b[1]+r*a[1]], [b[2];b[2]+r*a[2]], "--g")
        end

    end

    if lpdf.D == 3
        plot3D(c[1,:], c[2,:], c[3,:], "ro", alpha=0.9)
    else
        plot(c[1,:], c[2,:], "ro", alpha=0.9)
    end

end


text:
# I need to better understand the general case before designing a nice dispatch system...

"""
An abstract types which dictates the problem to be solved. For example, all problem with two spatial dimensions, all of which use 2D effective wavenumbers, we have the type TwoDimensions{T} <: PhysicalSetup{T}. In particular, for planar interface, we then have PlaneWaves2D{T} <: TwoDimensions{T}.
"""
abstract type ProblemSetup{T<:AbstractFloat, Dim<:Integer} end

abstract type PlaneWaves3D{T} <: ProblemSetup{T,3} end

"""Extract the dimension of the space that this physical property lives in"""
dim(p::ProblemSetup{T,Dim}) where {Dim,T} = Dim

"""
A setup with an incident plane wave, from the left, on a plate region filled with particles.
"""
struct Plate2D{T} <: ProblemSetup{T,2}
    """The plate is the region defined by x in the interval [x1, x2]"""
    x1::T
    x2::T
    """The background medium"""
    medium::Medium{T}
    """the species of the particles"""
    species::Vector{Specie{T}}
    """maximum order of hankel functions"""
    hankel_order::Int
    """minimal particle distance multiplier"""
    radius_multiplier::T
    """orientation of wave-vector (k*cos(θin), k*sin(θin))"""
    θin::T
end

"""
A setup with an incident plane wave, from the left, on a plate region filled with particles.
"""
struct Sphere{T} <: ProblemSetup{T,3}
    centre::T
    radius::T
    """The background medium"""
    medium::Medium{T}
    """the species of the particles"""
    species::Vector{Specie{T}}
    """maximum order of hankel functions"""
    hankel_order::Int
    """minimal particle distance multiplier"""
    radius_multiplier::T
    """orientation of wave-vector (k*cos(θin), k*sin(θin))"""
    θin::T
end


text:
#files.jl: Tools to help with file manipulations

#==Helper Functions
===============================================================================#
#Compute hash key from file contents to verify file equivalence.
function hashcontents(filepath::String)
	if !isfile(filepath)
		return sha2_256("")
	end
	open(filepath) do f
		return sha2_256(f)
	end
end

#SHA-based file check:
function filesmatch_SHA(src::String, dest::String)
	return hashcontents(src) == hashcontents(dest)
end

#Fast file matching algorithm.  Only checks dates & sizes.
#DEPRECATED
#   Cannot get copy operation to preserve time stamps.
function filesmatch_fast(src::String, dest::String)
	tstamp_tol = 10 #seconds
	result = false

	try
		#Argh... Windows cp -p function does not copy time down to msec
		result = isfile(dest) && abs(mtime(src) - mtime(dest)) <= tstamp_tol && filesize(src) == filesize(dest)
	finally
		return result
	end
end

#==Helper Functions
===============================================================================#
#Generate output file name from source filename
#(Strip out leading numbers, spaces, -, & _ from source basename)
function cleannamme(src::String)
	pat = r"^[0-9|\-|_| ]*(.*)$"

	#Get simplified filename:
	m = match(pat, src)
	result = strip(m.captures[1])
	if length(result) < 1
		result = src #Don't simplify
	end
	return result
end

#Last line

max_stars_repo_name: ma-laforge/AudioStick

max_issues_repo_name: ma-laforge/AudioStick

max_forks_repo_name: ma-laforge/AudioStick