text:
using LinearAlgebra
using KrylovKit
using Random

safesign(x::Number) = iszero(x) ? one(x) : sign(x)
"""
    qrpos(A)

Returns a QR decomposition, i.e. an isometric `Q` and upper triangular `R` matrix, where `R`
is guaranteed to have positive diagonal elements.
"""
qrpos(A) = qrpos!(copy(A))
function qrpos!(A)
    F = qr!(A)
    Q = Matrix(F.Q)
    R = F.R
    phases = safesign.(diag(R))
    rmul!(Q, Diagonal(phases))
    lmul!(Diagonal(conj!(phases)), R)
    return Q, R
end

"""
    lqpos(A)

Returns a LQ decomposition, i.e. a lower triangular `L` and isometric `Q` matrix, where `L`
is guaranteed to have positive diagonal elements.
"""
lqpos(A) = lqpos!(copy(A))
function lqpos!(A)
    F = qr!(Matrix(A'))
    Q = Matrix(Matrix(F.Q)')
    L = Matrix(F.R')
    phases = safesign.(diag(L))
    lmul!(Diagonal(phases), Q)
    rmul!(L, Diagonal(conj!(phases)))
    return L, Q
end

function cellones(Ni,Nj,D)
    Cell = Array{Array{Float64,2},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        Cell[i,j] = Matrix{Float64}(I, D, D)
    end
    return Cell
end

function ρmap(ρ,Ai,J)
    Nj = size(Ai,1)
    for j = 1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        ρ = ein"dc,csb,dsa -> ab"(ρ,Ai[jr],conj(Ai[jr]))
        # @tensor ρ[a,b] := ρ[a',b']*Ai[jr][b',s,b]*conj(Ai[jr][a',s,a])
    end
    return ρ
end

function initialA(M, D)
    Ni, Nj = size(M)
    T = eltype(M[1,1])
    A = Array{Array{Float64,3},2}(undef, Ni, Nj)
    for j in 1:Nj, i in 1:Ni
        d = size(M[i,j], 4)
        A[i,j] = rand(T, D, d, D)
    end
    return A
end
"""
    getL!(A,L; kwargs...)

````
┌ A1─A2─    ┌      L ─
ρ │  │    = ρ   =  │
┕ A1─A2─    ┕      L'─
````

ρ=L'*L, return L, where `L`is guaranteed to have positive diagonal elements.
L = cholesky!(ρ).U
If ρ is not exactly positive definite, cholesky will fail
"""
function getL!(A,L; kwargs...)
    Ni,Nj = size(A)
    for j = 1:Nj, i = 1:Ni
        _,ρs,_ = eigsolve(ρ->ρmap(ρ,A[i,:],j), L[i,j]'*L[i,j], 1, :LM; ishermitian = false, maxiter = 1, kwargs...)
        ρ = real(ρs[1] + ρs[1]')
        ρ ./= tr(ρ)
        F = svd!(ρ)
        Lo = lmul!(Diagonal(sqrt.(F.S)), F.Vt)
        _, L[i,j] = qrpos!(Lo)
    end
    return L
end

"""
    getAL(A,L)

Given an MPS tensor `A` and `L` ，return a left-canonical MPS tensor `AL`, a gauge transform `R` and
a scalar factor `λ` such that ``λ AR R = L A``
"""
function getAL(A,L)
    Ni,Nj = size(A)
    AL = Array{Array{Float64,3},2}(undef, Ni, Nj)
    Le = Array{Array{Float64,2},2}(undef, Ni, Nj)
    λ = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        D, d, = size(A[i,j])
        Q, R = qrpos!(reshape(L[i,j]*reshape(A[i,j], D, d*D), D*d, D))
        AL[i,j] = reshape(Q, D, d, D)
        λ[i,j] = norm(R)
        Le[i,j] = rmul!(R, 1/λ[i,j])
    end
    return AL, Le, λ
end

function getLsped(Le, A, AL; kwargs...)
    Ni,Nj = size(A)
    L = Array{Array{Float64,2},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        _, Ls, _ = eigsolve(X -> ein"dc,csb,dsa -> ab"(X,A[i,j],conj(AL[i,j])), Le[i,j], 1, :LM; ishermitian = false, kwargs...)
        _, L[i,j] = qrpos!(real(Ls[1]))
    end
    return L
end

"""
    leftorth(A,L=cellones(size(A,1),size(A,2),size(A[1,1],1)); tol = 1e-12, maxiter = 100, kwargs...)

Given an MPS tensor `A`, return a left-canonical MPS tensor `AL`, a gauge transform `L` and
a scalar factor `λ` such that ``λ AL L = L A``, where an initial guess for `L` can be
provided.
"""
function leftorth(A,L=cellones(size(A,1),size(A,2),size(A[1,1],1)); tol = 1e-12, maxiter = 100, kwargs...)
    L = getL!(A,L; kwargs...)
    AL, Le, λ= getAL(A,L;kwargs...)
    numiter = 1
    while norm(L.-Le) > tol && numiter < maxiter
        L = getLsped(Le, A, AL; kwargs...)
        AL, Le, λ = getAL(A, L; kwargs...)
        numiter += 1
    end
    L = Le
    return AL, L, λ
end

"""
    rightorth(A,L=cellones(size(A,1),size(A,2),size(A[1,1],1)); tol = 1e-12, maxiter = 100, kwargs...)

Given an MPS tensor `A`, return a gauge transform R, a right-canonical MPS tensor `AR`, and
a scalar factor `λ` such that ``λ R AR^s = A^s R``, where an initial guess for `R` can be
provided.
"""
function rightorth(A,L=cellones(size(A,1),size(A,2),size(A[1,1],1)); tol = 1e-12, maxiter = 100, kwargs...)
    Ni,Nj = size(A)
    Ar = Array{Array{Float64,3},2}(undef, Ni, Nj)
    Lr = Array{Array{Float64,2},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        Ar[i,j] = permutedims(A[i,j],(3,2,1))
        Lr[i,j] = permutedims(L[i,j],(2,1))
    end
    AL, L, λ = leftorth(Ar,Lr; tol = tol, kwargs...)
    R = Array{Array{Float64,2},2}(undef, Ni, Nj)
    AR = Array{Array{Float64,3},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        R[i,j] = permutedims(L[i,j],(2,1))
        AR[i,j] = permutedims(AL[i,j],(3,2,1))
    end
    return R, AR, λ
end

"""
    LRtoC(L,R)

```
 ── Cᵢⱼ ──  =  ── Lᵢⱼ ── Rᵢⱼ₊₁ ──
```
"""
function LRtoC(L,R)
    Ni, Nj = size(L)
    C = Array{Array{Float64,2},2}(undef, Ni, Nj)
    for j in 1:Nj,i in 1:Ni
        jr = j + 1 - (j + 1 > Nj) * Nj
        C[i,j] = L[i,j] * R[i,jr]
    end
    return C
end

"""
    FLm = FLmap(ALi, ALip, Mi, FL, J)

ALip means ALᵢ₊₁
```
  ┌──        ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁   ──   ...
  │          │     │        │
 FLm   =   FLᵢⱼ ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...
  │          │     │        │
  ┕──        ┕──  ALᵢ₊₁ⱼ ─ ALᵢ₊₁ⱼ₊₁ ──   ...
```
"""
function FLmap(ALi, ALip, Mi, FL, J)
    Nj = size(ALi,1)
    FLm = copy(FL)
    for j=1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        FLm = ein"abc,cde,bfhd,afg -> ghe"(FLm,ALi[jr],Mi[jr],conj(ALip[jr]))
    end
    return FLm
end

"""
    FRm = FRmap(ARi, ARip, Mi, FR, J)

ARip means ARᵢ₊₁
```
 ──┐       ... ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐
   │                │          │      │
──FRm  =   ... ──── Mᵢⱼ₋₁  ── Mᵢⱼ  ──FRᵢⱼ
   │                │          │      │
 ──┘       ... ─ ARᵢ₊₁ⱼ₋₁ ─ ARᵢ₊₁ⱼ  ──┘
```
"""
function FRmap(ARi, ARip, Mi, FR, J)
    Nj = size(ARi,1)
    FRm = copy(FR)
    for j=1:Nj
        jr = J-(j-1) + (J-(j-1) < 1)*Nj
        FRm = ein"abc,eda,hfbd,gfc -> ehg"(FRm,ARi[jr],Mi[jr],conj(ARip[jr]))
    end
    return FRm
end

function FLint(AL, M)
    Ni,Nj = size(AL)
    FL = Array{Array{Float64,3},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        D = size(AL[i,j],1)
        dL = size(M[i,j],1)
        FL[i,j] = rand(Float64, D, dL, D)
    end
    return FL
end

function FRint(AR, M)
    Ni,Nj = size(AR)
    FR = Array{Array{Float64,3},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        D = size(AR[i,j],1)
        dR = size(M[i,j],3)
        FR[i,j] = rand(Float64, D, dR, D)
    end
    return FR
end

"""
    λL, FL = leftenv(AL, M, FL = FLint(AL,M); kwargs...)

Compute the left environment tensor for MPS A and MPO M, by finding the left fixed point
of AL - M - conj(AL) contracted along the physical dimension.
```
 ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁   ──   ...         ┌──
 │     │        │                        │
FLᵢⱼ ─ Mᵢⱼ  ── Mᵢⱼ₊₁    ──   ...  = λLᵢⱼ FLᵢⱼ
 │     │        │                        │
 ┕──  ALᵢ₊₁ⱼ ─ ALᵢ₊₁ⱼ₊₁ ──   ...         ┕──
```
"""
leftenv(AL, M, FL = FLint(AL,M); kwargs...) = leftenv!(AL, M, copy(FL); kwargs...)
function leftenv!(AL, M, FL; kwargs...)
    Ni,Nj = size(AL)
    λL = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        λLs, FL1s, _= eigsolve(X->FLmap(AL[i,:], AL[ir,:], M[i,:], X, j), FL[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                FL[i,j] = real(FL1s[1])
                λL[i,j] = real(λLs[1])
            else
                FL[i,j] = real(FL1s[2])
                λL[i,j] = real(λLs[2])
            end
        else
            FL[i,j] = real(FL1s[1])
            λL[i,j] = real(λLs[1])
        end
    end
    return λL, FL
end

"""
    λR, FR = rightenv(AR, M, FR = FRint(AR,M); kwargs...)

Compute the right environment tensor for MPS A and MPO M, by finding the left fixed point
of AR - M - conj(AR) contracted along the physical dimension.
```
   ... ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐          ──┐
            │          │      │            │
   ... ──── Mᵢⱼ₋₁  ── Mᵢⱼ  ──FRᵢⱼ  = λRᵢⱼ──FRᵢⱼ
            │          │      │            │
   ... ─ ARᵢ₊₁ⱼ₋₁ ─ ARᵢ₊₁ⱼ  ──┘          ──┘
```
"""
rightenv(AR, M, FR = FRint(AR,M); kwargs...) = rightenv!(AR, M, copy(FR); kwargs...)
function rightenv!(AR, M, FR; kwargs...)
    Ni,Nj = size(AR)
    λR = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        λRs, FR1s, _= eigsolve(X->FRmap(AR[i,:], AR[ir,:], M[i,:], X, j), FR[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λRs) > 1 && norm(abs(λRs[1]) - abs(λRs[2])) < 1e-12
            @show λRs
            if real(λRs[1]) > 0
                FR[i,j] = real(FR1s[1])
                λR[i,j] = real(λRs[1])
            else
                FR[i,j] = real(FR1s[2])
                λR[i,j] = real(λRs[2])
            end
        else
            FR[i,j] = real(FR1s[1])
            λR[i,j] = real(λRs[1])
        end
    end
    return λR, FR
end

"""
    ACm = ACmap(ACij, FLj, FRj, Mj, II)

```
                                ┌─────── ACᵢⱼ ─────┐
                                │        │         │
┌─────── ACm  ─────┐      =     FLᵢⱼ ─── Mᵢⱼ ───── FRᵢⱼ
│        │         │            │        │         │
                                FLᵢ₊₁ⱼ ─ Mᵢ₊₁ⱼ ──  FRᵢ₊₁ⱼ
                                │        │         │
                                .        .         .
                                .        .         .
                                .        .         .
```
"""
function ACmap(ACij, FLj, FRj, Mj, II)
    Ni = size(FLj,1)
    ACm = copy(ACij)
    for i=1:Ni
        ir = II+i-1 - (II+i-1 > Ni)*Ni
        ACm = ein"abc,cde,bhfd,efg -> ahg"(FLj[ir],ACm,Mj[ir],FRj[ir])
    end
    return ACm
end

"""
    Cmap(Cij, FLjp, FRj, II)

```
                    ┌────Cᵢⱼ ───┐
                    │           │
┌──── Cm ───┐   =   FLᵢⱼ₊₁ ──── FRᵢⱼ
│           │       │           │
                    FLᵢ₊₁ⱼ₊₁ ── FRᵢ₊₁ⱼ
                    │           │
                    .           .
                    .           .
                    .           .
```
"""
function Cmap(Cij, FLjp, FRj, II)
    Ni = size(FLjp,1)
    Cm = copy(Cij)
    for i=1:Ni
        ir = II+i-1 - (II+i-1 > Ni)*Ni
        Cm = ein"abc,cd,dbe -> ae"(FLjp[ir],Cm,FRj[ir])
    end
    return Cm
end

"""
    ACenv(AC, FL, M, FR;kwargs...)

Compute the up environment tensor for MPS `FL`,`FR` and MPO `M`, by finding the up fixed point
        of `FL - M - FR` contracted along the physical dimension.
```
┌─────── ACᵢⱼ ─────┐
│        │         │
FLᵢⱼ ─── Mᵢⱼ ───── FRᵢⱼ
│        │         │
FLᵢ₊₁ⱼ ─ Mᵢ₊₁ⱼ ──  FRᵢ₊₁ⱼ  =  λACᵢⱼ ┌──── ACᵢⱼ ───┐
│        │         │                │      │      │
.        .         .
.        .         .
.        .         .
```
"""
ACenv(AC, FL, M, FR; kwargs...) = ACenv!(copy(AC), FL, M, FR; kwargs...)
function ACenv!(AC, FL, M, FR; kwargs...)
    Ni,Nj = size(AC)
    λAC = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        λACs, ACs, = eigsolve(X->ACmap(X, FL[:,j], FR[:,j], M[:,j], i), AC[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λACs) > 1 && norm(abs(λACs[1]) - abs(λACs[2])) < 1e-12
            @show λACs
            if real(λACs[1]) > 0
                AC[i,j] = real(ACs[1])
                λAC[i,j] = real(λACs[1])
            else
                AC[i,j] = real(ACs[2])
                λAC[i,j] = real(λACs[2])
            end
        else
            AC[i,j] = real(ACs[1])
            λAC[i,j] = real(λACs[1])
        end
    end
    return λAC, AC
end

"""
    Cenv(C, FL, FR;kwargs...)

Compute the up environment tensor for MPS `FL` and `FR`, by finding the up fixed point
    of `FL - FR` contracted along the physical dimension.
```
┌────Cᵢⱼ ───┐
│           │
FLᵢⱼ₊₁ ──── FRᵢⱼ
│           │
FLᵢ₊₁ⱼ₊₁ ── FRᵢ₊₁ⱼ   =  λCᵢⱼ ┌──Cᵢⱼ ─┐
│           │                │       │
.           .
.           .
.           .
```
"""
Cenv(C, FL, FR; kwargs...) = Cenv!(copy(C), FL, FR; kwargs...)
function Cenv!(C, FL, FR; kwargs...)
    Ni,Nj = size(C)
    λC = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        jr = j + 1 - (j==Nj) * Nj
        λCs, Cs, = eigsolve(X->Cmap(X, FL[:,jr], FR[:,j], i), C[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λCs) > 1 && norm(abs(λCs[1]) - abs(λCs[2])) < 1e-12
            @show λCs
            if real(λCs[1]) > 0
                C[i,j] = real(Cs[1])
                λC[i,j] = real(λCs[1])
            else
                C[i,j] = real(Cs[2])
                λC[i,j] = real(λCs[2])
            end
        else
            C[i,j] = real(Cs[1])
            λC[i,j] = real(λCs[1])
        end
    end
    return λC, C
end

function ACCtoAL(ACij,Cij)
    D,d, = size(ACij)
    QAC, _ = qrpos(reshape(ACij,(D*d, D)))
    QC, _ = qrpos(Cij)
    reshape(QAC*QC', (D, d, D))
end

function ACCtoAR(ACij,Cijr)
    D,d, = size(ACij)
    _, QAC = lqpos(reshape(ACij,(D, d*D)))
    _, QC = lqpos(Cijr)
    reshape(QC'*QAC, (D, d, D))
end

"""
    itoir(i,Ni,Nj)

````
i -> (i,j) -> (i,jr) -> ir
````
"""
function itoir(i,Ni,Nj)
    Liner = LinearIndices((1:Ni,1:Nj))
    Cart = CartesianIndices((1:Ni,1:Nj))
    Index = Cart[i]
    i,j = Index[1],Index[2]
    jr = j - 1 + (j==1)*Nj
    Liner[i,jr]
end

function ALCtoAC(AL,C)
    Ni,Nj = size(AL)
    ACij = [ein"asc,cb -> asb"(AL[i],C[i]) for i=1:Ni*Nj]
    reshape(ACij,Ni,Nj)
end

"""
    ACCtoALAR(AL, C, AR, M, FL, FR; kwargs...)

QR factorization to get `AL` and `AR` from `AC` and `C`

````
──ALᵢⱼ──Cᵢⱼ──  =  ──ACᵢⱼ──  = ──Cᵢ₋₁ⱼ ──ARᵢⱼ──
  │                  │                  │
````
"""
function ACCtoALAR(AL, C, AR, M, FL, FR; kwargs...)
    Ni,Nj = size(AL)
    AC = ALCtoAC(AL,C)
    _, AC = ACenv(AC, FL, M, FR; kwargs...)
    _, C = Cenv(C, FL, FR; kwargs...)

    ALij = [ACCtoAL(AC[i],C[i]) for i=1:Ni*Nj]
    AL = reshape(ALij,Ni,Nj)
    ARij = [ACCtoAR(AC[i],C[itoir(i,Ni,Nj)]) for i=1:Ni*Nj]
    AR = reshape(ARij,Ni,Nj)
    return AL, C, AR
end

"""
Compute the error through all environment `AL,C,FL,M,FR`

````
        ┌── AC──┐
        │   │   │           ┌── AC──┐
MAC1 =  FL─ M ──FR  =  λAC  │   │   │
        │   │   │

        ┌── AC──┐
        │   │   │           ┌──C──┐
MAC2 =  FL─ M ──FR  =  λAC  │     │
        │   │   │
        ┕── AL─

── MAC1 ──    ≈    ── AL ── MAC2 ──
    │                 │
````
"""
function error(AL,C,FL,M,FR)
    Ni,Nj = size(AL)
    AC = Array{Array{Float64,3},2}(undef, Ni,Nj)
    err = 0
    for j = 1:Nj,i = 1:Ni
        AC[i,j] = ein"asc,cb -> asb"(AL[i,j],C[i,j])
        MAC = ACmap(AC[i,j], FL[:,j], FR[:,j], M[:,j], i)
        MAC -= ein"asd,cpd,cpb -> asb"(AL[i,j],conj(AL[i,j]),MAC)
        err += norm(MAC)
    end
    return err
end

"""
    obs_env()

If `Ni,Nj>1` and `Mij` are different bulk tensor, the up and down environment are different. So to calculate observable, we must get ACup and ACdown, which is easy to get by overturning the `Mij`. Then be cautious to get the new `FL` and `FR` environment.
"""
function obs_env() end

"""
    λL, FL = obs2x2FL(AL, M, FL = FLint(AL,M); kwargs...)

This function is designed specifically for 2x2 longitudinally symmetric cell to get correct `FL` environment.
```
 ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁   ──   ...         ┌──
 │     │        │                        │
FLᵢⱼ ─ Mᵢⱼ  ── Mᵢⱼ₊₁    ──   ...  = λLᵢⱼ FLᵢⱼ
 │     │        │                        │
 ┕──  ALᵢⱼ  ─  ALᵢⱼ₊₁   ──   ...         ┕──
```
"""
obs2x2FL(AL, M, FL = FLint(AL,M); kwargs...) = obs2x2FL!(AL, M, copy(FL); kwargs...)
function obs2x2FL!(AL, M, FL; kwargs...)
    Ni,Nj = size(AL)
    λL = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        λLs, FL1s, _= eigsolve(X->FLmap(AL[i,:], AL[i,:], M[i,:], X, j), FL[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                FL[i,j] = real(FL1s[1])
                λL[i,j] = real(λLs[1])
            else
                FL[i,j] = real(FL1s[2])
                λL[i,j] = real(λLs[2])
            end
        else
            FL[i,j] = real(FL1s[1])
            λL[i,j] = real(λLs[1])
        end
    end
    return λL, FL
end

"""
    λR, FR = obs2x2FR(AR, M, FR = FRint(AR,M); kwargs...)

This function is designed specifically for 2x2 longitudinally symmetric cell to get correct `FL` environment.
```
   ... ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐          ──┐
            │          │      │            │
   ... ──── Mᵢⱼ₋₁  ── Mᵢⱼ  ──FRᵢⱼ  = λRᵢⱼ──FRᵢⱼ
            │          │      │            │
   ... ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┘          ──┘
```
"""
obs2x2FR(AR, M, FR = FRint(AR,M); kwargs...) = obs2x2FR!(AR, M, copy(FR); kwargs...)
function obs2x2FR!(AR, M, FR; kwargs...)
    Ni,Nj = size(AR)
    λR = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        λRs, FR1s, _= eigsolve(X->FRmap(AR[i,:], AR[i,:], M[i,:], X, j), FR[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λRs) > 1 && norm(abs(λRs[1]) - abs(λRs[2])) < 1e-12
            @show λRs
            if real(λRs[1]) > 0
                FR[i,j] = real(FR1s[1])
                λR[i,j] = real(λRs[1])
            else
                FR[i,j] = real(FR1s[2])
                λR[i,j] = real(λRs[2])
            end
        else
            FR[i,j] = real(FR1s[1])
            λR[i,j] = real(λRs[1])
        end
    end
    return λR, FR
end

function BgFLint(AL, M)
    Ni,Nj = size(AL)
    BgFL = Array{Array{Float64,4},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        D1 = size(AL[i,j],1)
        D2 = size(AL[irr,j],1)
        dL1 = size(M[i,j],1)
        dL2 = size(M[ir,j],1)
        BgFL[i,j] = rand(Float64, D1, dL1, dL2, D2)
    end
    return BgFL
end

"""
    BgFLm = BgFLmap(ALi, ALip, Mi, Mip, BgFLij, J)

ALip means ALᵢ₊₂
Mip means Mᵢ₊₁
```
  ┌──        ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁    ──   ...
  │          │     │        │
  │          │  ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...
BgFLm   =  BgFLᵢⱼ  │        │
  │          │  ─ Mᵢ₊₁ⱼ ── Mᵢ₊₁ⱼ₊₁   ──   ...
  │          │     │        │
  ┕──        ┕──  ALᵢ₊₂ⱼ ─ ALᵢ₊₂ⱼ₊₁  ──   ...
```
"""
function BgFLmap(ALi, ALip, Mi, Mip, BgFLij, J)
    Nj = size(ALi,1)
    BgFLm = copy(BgFLij)
    for j=1:Nj
        jr = J+j-1 - (J+j-1 > Nj)*Nj
        BgFLm = ein"dcba,def,ckge,bjhk,aji -> fghi"(BgFLm,ALi[jr],Mi[jr],Mip[jr],conj(ALip[jr]))
    end
    return BgFLm
end

"""
    λL, BgFL = bigleftenv(AL, M, BgFL = BgFLint(AL,M); kwargs...)

Compute the left environment tensor for MPS A and MPO M, by finding the left fixed point
of AL - M - M - conj(AL) contracted along the physical dimension.
```
   ┌──  ALᵢⱼ  ── ALᵢⱼ₊₁    ──   ...           ┌──
   │     │        │                           │
   │  ─ Mᵢⱼ   ── Mᵢⱼ₊₁     ──   ...           │
 BgFLᵢⱼ  │        │                   = λLᵢⱼ BgFLᵢⱼ
   │  ─ Mᵢ₊₁ⱼ ── Mᵢ₊₁ⱼ₊₁   ──   ...           │
   │     │        │                           │
   ┕──  ALᵢ₊₂ⱼ ─ ALᵢ₊₂ⱼ₊₁  ──   ...           ┕──
```
"""
bigleftenv(AL, M, BgFL = BgFLint(AL,M); kwargs...) = bigleftenv!(AL, M, copy(BgFL); kwargs...)
function bigleftenv!(AL, M, BgFL; kwargs...)
    Ni,Nj = size(AL)
    λL = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        # irr = i + 2 - Ni * (i + 2 > Ni) # modified for 2x2
        λLs, BgFL1s, _= eigsolve(X->BgFLmap(AL[i,:], AL[ir,:], M[i,:], M[ir,:], X, j), BgFL[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λLs) > 1 && norm(abs(λLs[1]) - abs(λLs[2])) < 1e-12
            @show λLs
            if real(λLs[1]) > 0
                BgFL[i,j] = real(BgFL1s[1])
                λL[i,j] = real(λLs[1])
            else
                BgFL[i,j] = real(BgFL1s[2])
                λL[i,j] = real(λLs[2])
            end
        else
            BgFL[i,j] = real(BgFL1s[1])
            λL[i,j] = real(λLs[1])
        end
    end
    return λL, BgFL
end

function BgFRint(AR, M)
    Ni,Nj = size(AR)
    BgFR = Array{Array{Float64,4},2}(undef, Ni, Nj)
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        D1 = size(AR[i,j],3)
        D2 = size(AR[irr,j],3)
        dR1 = size(M[i,j],3)
        dR2 = size(M[ir,j],3)
        BgFR[i,j] = rand(Float64, D1, dR1, dR2, D2)
    end
    return BgFR
end

"""
    FRm = FRmap(ARi, ARip, Mi, FR, J)

ARip means ARᵢ₊₁
```
 ──┐          ...  ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐
   │                    │          │      │
 ──│          ... ──── Mᵢⱼ₋₁  ──  Mᵢⱼ   ──│
  BgFRm   =             │          │     BgFRm
 ──│          ... ──── Mᵢ₊₁ⱼ₋₁ ── Mᵢ₊₁ⱼ ──│
   │                    │          │      │
 ──┘          ...  ─ ARᵢ₊₂ⱼ₋₁ ─── ARᵢ₊₂ⱼ──┘
```
"""
function BgFRmap(ARi, ARip, Mi, Mip, BgFR, J)
    Nj = size(ARi,1)
    BgFRm = copy(BgFR)
    for j=1:Nj
        jr = J-(j-1) + (J-(j-1) < 1)*Nj
        BgFRm = ein"fghi,def,ckge,bjhk,aji -> dcba"(BgFRm,ARi[jr],Mi[jr],Mip[jr],conj(ARip[jr]))
    end
    return BgFRm
end

"""
    λR, BgFR = bigrightenv(AR, M, BgFR = BgFRint(AR,M); kwargs...)

Compute the right environment tensor for MPS A and MPO M, by finding the left fixed point
of AR - M - M - conj(AR) contracted along the physical dimension.
```
     ──┐          ...  ─── ARᵢⱼ₋₁  ── ARᵢⱼ  ──┐
       │                    │          │      │
     ──│          ... ──── Mᵢⱼ₋₁  ──  Mᵢⱼ   ──│
λRᵢⱼ BgFRᵢⱼ   =             │          │     BgFRᵢⱼ
     ──│          ... ──── Mᵢ₊₁ⱼ₋₁ ── Mᵢ₊₁ⱼ ──│
       │                    │          │      │
     ──┘          ...  ─ ARᵢ₊₂ⱼ₋₁ ─── ARᵢ₊₂ⱼ──┘
```
"""
bigrightenv(AR, M, BgFR = BgFRint(AR,M); kwargs...) = bigrightenv!(AR, M, copy(BgFR); kwargs...)
function bigrightenv!(AR, M, BgFR; kwargs...)
    Ni,Nj = size(AR)
    λR = zeros(Ni,Nj)
    for j = 1:Nj,i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        # irr = i + 2 - Ni * (i + 2 > Ni) # modified for 2x2
        λRs, BgFR1s, _= eigsolve(X->BgFRmap(AR[i,:], AR[ir,:], M[i,:], M[ir,:], X, j), BgFR[i,j], 1, :LM; ishermitian = false, kwargs...)
        if length(λRs) > 1 && norm(abs(λRs[1]) - abs(λRs[2])) < 1e-12
            @show λRs
            if real(λRs[1]) > 0
                BgFR[i,j] = real(BgFR1s[1])
                λR[i,j] = real(λRs[1])
            else
                BgFR[i,j] = real(BgFR1s[2])
                λR[i,j] = real(λRs[2])
            end
        else
            BgFR[i,j] = real(BgFR1s[1])
            λR[i,j] = real(λRs[1])
        end
    end
    return λR, BgFR
end


text:
abstract type MCDMResult end
abstract type SCDMResult end

struct TopsisResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    normalizedDecisionMatrix::DataFrame
    normalizedWeightedDecisionMatrix::DataFrame
    bestIndex::Int64
    scores::Array{Float64,1}
end

struct VikorResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    weightedDecisionMatrix::DataFrame
    bestIndex::Int64
    scores::Array{Float64,1}
end

struct ElectreResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    weightedDecisionMatrix::DataFrame
    fitnessTable::Array{Dict,1}
    nonfitnessTable::Array{Dict,1}
    fitnessMatrix::Array{Float64,2}
    nonfitnessMatrix::Array{Float64,2}
    C::Array{Float64,1}
    D::Array{Float64,1}
    bestIndex::Tuple
end

struct MooraResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    weightedDecisionMatrix::DataFrame
    referenceMatrix::DataFrame
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct DematelResult <: MCDMResult
    comparisonMatrix::DataFrame
    NDMatrix::Array{Float64,2}
    relationShipMatrix::Array{Float64,2}
    c::Array{Float64,1}
    r::Array{Float64,1}
    c_plus_r::Array{Float64,1}
    c_minus_r::Array{Float64,1}
    threshold::Float64
    influenceMatrix::Array{Float64,2}
    weights::Array{Float64,1}
end

struct AHPConsistencyResult <: MCDMResult
    comparisonMatrix::DataFrame
    normalizedComparisonMatrix::DataFrame
    consistencyVector::Array{Float64,1}
    priority::Array{Float64,1}
    pc::Array{Float64,1}
    lambda_max::Float64
    CI::Float64
    RI::Float64
    CR::Float64
    isConsistent::Bool
end

struct AHPResult <: MCDMResult
    comparisonMatrixList::Array{DataFrame,1}
    criteriaComparisonMatrix::DataFrame
    criteriaConsistency::AHPConsistencyResult
    decisionMatrix::DataFrame
    scores::Array{Float64,1}
    weights::Array{Float64,1}
    bestIndex::Int64
end

struct NDSResult <: MCDMResult
    ranks::Array{Int64,1}
    bestIndex::Int64
end

struct LaplaceResult <: SCDMResult
    expected_values::Array{Float64,1}
    bestIndex::Int64
end

struct MaximinResult <: SCDMResult
    rowmins::Array{Float64,1}
    bestIndex::Int64
end

struct MaximaxResult <: SCDMResult
    rowmaxs::Array{Float64,1}
    bestIndex::Int64
end

struct MinimaxResult <: SCDMResult
    rowmaxs::Array{Float64,1}
    bestIndex::Int64
end

struct MiniminResult <: SCDMResult
    rowmins::Array{Float64,1}
    bestIndex::Int64
end

struct SavageResult <: SCDMResult
    regretMatrix::DataFrame
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct HurwiczResult <: SCDMResult
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct MLEResult <: SCDMResult
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct ExpectedRegretResult <: SCDMResult
    scores::Array{Float64,1}
    bestIndex::Int64
end

struct GameResult <: MCDMResult
    row_player_probabilities::Array{Float64,1}
    value::Float64
end

struct DataEnvelopResult <: MCDMResult
    efficiencies::Array{Float64,1}
    references::DataFrame
    orderedcases::Array{Symbol,1}
end

struct GreyResult <: MCDMResult
    referenceRow::Array{Float64,1}
    normalizedMat::Matrix
    absoluteValueMat::Matrix
    greyTable::Matrix
    scores::Array{Float64,1}
    ordering::Array{Int64,1}
    bestIndex::Int64
end

struct SawResult <: MCDMResult
    decisionMatrix::DataFrame
    normalizedDecisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct ARASResult <: MCDMResult
    referenceRow::Array{Float64,1}
    extendMat::Array{Float64,2}
    normalizedMat::Array{Float64,2}
    optimalitydegrees::Array{Float64,1}
    scores::Array{Float64,1}
    orderings::Array{Int64,1}
    bestIndex::Int64
end

struct WPMResult <: MCDMResult
    decisionMatrix::DataFrame
    normalizedDecisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct WASPASResult <: MCDMResult
    decisionMatrix::DataFrame
    normalizedDecisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct EDASResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct MARCOSResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct MABACResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct MAIRCAResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct COPRASResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct PrometheeResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct CoCoSoResult <: MCDMResult
    decisionMatrix::DataFrame
    weights::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct CRITICResult <: MCDMResult
    decisionMatrix::DataFrame
    w::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct EntropyResult <: MCDMResult
    decisionMatrix::DataFrame
    w::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

struct CODASResult <: MCDMResult
    decisionMatrix::DataFrame
    w::Array{Float64,1}
    scores::Array{Float64,1}
    ranking::Array{Int64,1}
    bestIndex::Int64
end

ext: jl

max_stars_repo_name: bahadirfyildirim/JMcDM

max_issues_repo_name: bahadirfyildirim/JMcDM

max_forks_repo_name: bahadirfyildirim/JMcDM


text:
## parse.jl : functions to parse input arguments of the Estimator
#
# Parse inputs passed to the Estimator type.
#
# This file is part of MultilevelEstimators.jl - A Julia toolbox for
# Multilevel Monte Carlo Methods (c) Pieterjan Robbe, 2019

## parse! ##
# Make a function parse! for the key key_name.
# This macro has three arguments:
#  - an expression with the key name;
#  - an expression with the default value;
#  - an expression with the checks to perform on this input argument.
macro parse!(key_name, default_value, checks_to_perform)
    @eval begin
        function parse!(index_set::AbstractIndexSet, sample_method::AbstractSampleMethod, options::Dict{Symbol, T} where T, symbol::Val{$key_name})
            key = eltype(symbol)
            if !haskey(options, key)
                options[key] = $default_value
            else
                val = options[key]
                $checks_to_perform
            end
        end
    end
end

parse!(index_set::AbstractIndexSet, sample_method::AbstractSampleMethod, options::Dict{Symbol, T} where T, symbol::Symbol) = parse!(index_set, sample_method, options, Val(symbol))

to_string(key, val) = Estimator, val, string("optional key ", key)

eltype(::Type{<:Val{T}}) where {T} = T

## nb_of_warm_up_samples ##
@parse!(:nb_of_warm_up_samples,
		sample_method isa QMC ? (index_set isa U ? 2 : 1) : 20,
        begin
            check_type(to_string(key, val)..., Signed)
            check_larger_than(to_string(key, val)..., 1)
        end
       )

## nb_of_qoi ##
@parse!(:nb_of_qoi,
        1,
        begin
            check_type(to_string(key, val)..., Signed)
            check_larger_than(to_string(key, val)..., 0)
        end
       )

## continuate ##
@parse!(:continuate,
        true,
        check_type(to_string(key, val)..., Bool)
       )

## nb_of_tols ##
@parse!(:nb_of_tols,
        10,
        begin
            check_type(to_string(key, val)..., Integer)
            check_larger_than(to_string(key, val)..., 0)
        end
       )

## continuation_mul_factor ##
@parse!(:continuation_mul_factor,
        1.3,
        begin
            check_type(to_string(key, val)..., Real)
            check_finite(to_string(key, val)...)
            check_larger_than(to_string(key, val)..., 1)
        end
       )

## min_splitting ##
@parse!(:min_splitting,
        0.5,
        begin
            check_type(to_string(key, val)..., Real)
            check_finite(to_string(key, val)...)
            check_larger_than(to_string(key, val)..., 0)
            check_smaller_than_or_equal_to(to_string(key, val)..., 1)
            haskey(options, :max_splitting) || parse!(index_set, sample_method, options, :max_splitting)
            check_ordered_or_equal(Estimator, val, options[:max_splitting], "optional key min_splitting", "optional key max_splitting")

        end
       )

## max_splitting ##
@parse!(:max_splitting,
        0.99,
        begin
            check_type(to_string(key, val)..., Real)
            check_finite(to_string(key, val)...)
            check_larger_than(to_string(key, val)..., 0)
            check_smaller_than_or_equal_to(to_string(key, val)..., 1)
            haskey(options, :min_splitting) || parse!(index_set, sample_method, options, :min_splitting)
            check_ordered_or_equal(Estimator, options[:min_splitting], val, "optional key min_splitting", "optional key max_splitting")
        end
       )

## folder ##
@parse!(:folder,
        pwd(),
        begin
            check_type(to_string(key, val)..., String)
            haskey(options, :save) || parse!(index_set, sample_method, options, :save)
            if options[:save]
                isdir(val) || throw(ArgumentError(string(val, "is not a directory!")))
                ispath(val) || makepath(val)
            end
        end
       )

## name ##
@parse!(:name,
        get_valid_filename(index_set, sample_method, options),
        begin
            check_type(to_string(key, val)..., String)
            parse!(index_set, sample_method, options, Val(:folder))
            !endswith(val, ".jld2") && occursin(".", val) && throw(ArgumentError("in Estimator, optional key name must not contain a ."))
            val = endswith(val, ".jld2") ? val : string(val, ".jld2")
            options[key] = val
            isfile(joinpath(options[:folder], val)) && @warn string("filename ", val, " exists, will be overwritten!")
        end
       )

function get_valid_filename(index_set, sample_method, options)
    parse!(index_set, sample_method, options, Val(:folder))
    filename = "UntitledEstimator"
    cntr = 0
    if isfile(joinpath(options[:folder], string(filename, ".jld2")))
        cntr += 1
        while isfile(joinpath(options[:folder], string(filename, cntr, ".jld2")))
            cntr = cntr+1
        end
    end
    string(filename, cntr == 0 ? "" : cntr, ".jld2")
end

## save_samples ##
@parse!(:save_samples,
        false,
        check_type(to_string(key, val)..., Bool)
       )

## save ##
@parse!(:save,
        true,
        check_type(to_string(key, val)..., Bool)
       )

## verbose ##
@parse!(:verbose,
        true,
        check_type(to_string(key, val)..., Bool)
       )

## cost_model ##
@parse!(:cost_model,
        EmptyFunction(),
        check_type(to_string(key, val)..., Function)
       )

struct EmptyFunction <: Function end

## robustify_bias_estimate ##
@parse!(:robustify_bias_estimate,
        true,
        check_type(to_string(key, val)..., Bool)
       )

## do_mse_splitting ##
@parse!(:do_mse_splitting,
        true,
        check_type(to_string(key, val)..., Bool)
       )

## max_index_set_param ##
@parse!(:max_index_set_param,
        10 * ndims(index_set),
        begin
            check_type(to_string(key, val)..., Signed)
            check_larger_than(to_string(key, val)..., 0)
            haskey(options, :min_index_set_param) || parse!(index_set, sample_method, options, :min_index_set_param)
            check_ordered_or_equal(Estimator, options[:min_index_set_param], val, "optional key min_index_set_param", "optional key max_index_set_param")
        end
       )

## min_index_set_param ##
@parse!(:min_index_set_param,
        2,
        begin
            check_type(to_string(key, val)..., Signed)
            check_larger_than(to_string(key, val)..., 0)
            haskey(options, :max_index_set_param) || parse!(index_set, sample_method, options, :max_index_set_param)
            check_ordered_or_equal(Estimator, val, options[:max_index_set_param], "optional key min_index_set_param", "optional key max_index_set_param")
        end
       )

## sample_mul_factor ##
@parse!(:sample_mul_factor,
        1.2,
        begin
            check_type(to_string(key, val)..., Real)
            check_finite(to_string(key, val)...)
        end
       )

## nb_of_workers ##
@parse!(:nb_of_workers,
        i -> nworkers(),
        begin
            check_type(to_string(key, val)..., Union{Integer, Function})
            if val isa Integer
                check_larger_than(to_string(key, val)..., 0)
                delete!(options, val)
                options[key] = i -> val
            end
        end
       )

## nb_of_shifts ##
@parse!(:nb_of_shifts,
        i -> 10,
        begin
            check_type(to_string(key, val)..., Union{Integer, Function})
            val isa Function && index_set isa U && throw(ArgumentError("in Estimator, optional key nb_of_shifts cannot be a Function when using unbiased estimation"))
            if val isa Integer
                check_larger_than(to_string(key, val)..., 1)
                delete!(options, val)
                options[key] = i -> val
            end
        end
       )

## point_generator ##
@parse!(:point_generator,
        LatticeRule32(length(options[:distributions])),
        check_type(to_string(key, val)..., AbstractRNG)
       )

## do_regression ##
@parse!(:do_regression,
        true,
        check_type(to_string(key, val)..., Bool)
       )

## max_search_space ##
@parse!(:max_search_space,
        TD(ndims(index_set)),
        begin
            check_type(to_string(key, val)..., AbstractIndexSet)
            check_equal_to(Estimator, ndims(val), "dimensions of max_search_space", ndims(index_set))
        end
       )

## nb_of_uncertainties ##
@parse!(:nb_of_uncertainties,
        i -> length(options[:distributions]),
        check_type(to_string(key, val)..., Function)
       )

## penalization ##
@parse!(:penalization,
        0.5,
        begin
            check_type(to_string(key, val)..., Real)
            check_finite(to_string(key, val)...)
            check_larger_than_or_equal_to(to_string(key, val)..., 0)
            check_smaller_than_or_equal_to(to_string(key, val)..., 1)
        end
       )

## acceptance_rate ##
@parse!(:acceptance_rate,
        1.0,
        begin
            check_type(to_string(key, val)..., Real)
            check_finite(to_string(key, val)...)
            check_larger_than_or_equal_to(to_string(key, val)..., 0)
            check_smaller_than_or_equal_to(to_string(key, val)..., 1)
        end
       )

## qoi_with_max_var
@parse!(:qoi_with_max_var,
        0,
        begin
            check_larger_than(to_string(key, val)..., 0)
            haskey(options, :nb_of_qoi) || parse!(index_set, sample_method, options, :nb_of_qoi)
            check_smaller_than_or_equal_to(to_string(key, val)..., length(options[:nb_of_qoi]))
        end
       )


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CreateOrUpdateTrustedIdProviderParameters <: SwaggerModel
    properties::Any # spec type: Union{ Nothing, CreateOrUpdateTrustedIdProviderProperties } # spec name: properties

    function CreateOrUpdateTrustedIdProviderParameters(;properties=nothing)
        o = new()
        validate_property(CreateOrUpdateTrustedIdProviderParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type CreateOrUpdateTrustedIdProviderParameters

const _property_map_CreateOrUpdateTrustedIdProviderParameters = Dict{Symbol,Symbol}(Symbol("properties")=>Symbol("properties"))
const _property_types_CreateOrUpdateTrustedIdProviderParameters = Dict{Symbol,String}(Symbol("properties")=>"CreateOrUpdateTrustedIdProviderProperties")
Base.propertynames(::Type{ CreateOrUpdateTrustedIdProviderParameters }) = collect(keys(_property_map_CreateOrUpdateTrustedIdProviderParameters))
Swagger.property_type(::Type{ CreateOrUpdateTrustedIdProviderParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateOrUpdateTrustedIdProviderParameters[name]))}
Swagger.field_name(::Type{ CreateOrUpdateTrustedIdProviderParameters }, property_name::Symbol) =  _property_map_CreateOrUpdateTrustedIdProviderParameters[property_name]

function check_required(o::CreateOrUpdateTrustedIdProviderParameters)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

function validate_property(::Type{ CreateOrUpdateTrustedIdProviderParameters }, name::Symbol, val)
end


text:
set_device(dev_name::T, pid::Int) where {T<:AbstractString} = Distributed.remotecall_eval(ImagineWorker, pid, :(_set_device($dev_name)))

#sets device for driver process and all workers
function set_device(dev_name::T) where T<:AbstractString
    global DEFAULT_DEVICE = String(dev_name)
    map(pid->set_device(dev_name, pid), WORKERS)
end

function load_worker_code(id)
    #below was suggested by a julia issue but seems unnecessary now
    #Distributed.remotecall_eval(Main, id, :(using Pkg))
    #Distributed.remotecall_eval(Main, id, :(Pkg.activate(".")))
    Distributed.remotecall_eval(Main, id, :(using ImagineInterface))
    Distributed.remotecall_eval(Main, id, :(using ImagineWorker))
end

function add_workers(n::Int, dev=DEFAULT_DEVICE)
    ids = addprocs(n; dir=joinpath(@__DIR__, ".."))
    append!(FREE_WORKERS, ids)
    append!(WORKERS, ids)
    print("Initializing $n new worker processes...\n")
    load_worker_code(ids[1])
    @sync begin
        for i = 2:length(ids)
	    @async load_worker_code(ids[i])
        end
    end
    for i = 1:length(ids)
        set_device(dev, ids[i])
    end
    return ids
end

function get_worker(; dev=DEFAULT_DEVICE)
    proc_id = -1
    if isempty(FREE_WORKERS)
        print("No free workers found.  Adding more.\n")
        add_workers(1, dev)[1]
        proc_id = pop!(FREE_WORKERS)
    else
        proc_id = pop!(FREE_WORKERS)
        set_device(dev, proc_id)
    end
    push!(USED_WORKERS, proc_id)
    return proc_id
end

function free_workers(proc_ids::Vector{Int})
    for proc_id in proc_ids
        if !in(proc_id, WORKERS)
            error("Worker id $proc_id is not in this module's worker pool")
        end
        if length(WORKERS) > NPERSISTENT_WORKERS
            print("Freeing worker AND removing process\n")
            rmprocs(proc_id)
            deleteat!(WORKERS, findall(x->x==proc_id, WORKERS))
            deleteat!(USED_WORKERS, findall(x->x==proc_id, USED_WORKERS))
        else
            push!(FREE_WORKERS, proc_id)
            deleteat!(USED_WORKERS, findall(x->x==proc_id, USED_WORKERS))
        end
    end
    return FREE_WORKERS
end
free_workers(proc_id::Int) = free_workers([proc_id;])


text:
using Plots

@fastmath function SI4R_fs(N::Int64,bet::Float64,gamm::Float64)
    Psi = (N+1)*(N+2)*(N+3)*(N+4)/24;
    final_size = zeros(N+1,1);
    p_vec = zeros(Psi,1);
    p_vec[2]=1;
    @inbounds for Z5 in 0:N
        w::Int64 = Psi - (N-Z5+1)*(N-Z5+2)*(N-Z5+3)*(N-Z5+4)/24 + 1;
        @inbounds for Z4 in Z5:N
            a5 = 4*gamm*(Z4-Z5);
            @inbounds for Z3 in Z4:N
                a4 = 4*gamm*(Z3-Z4);
                @inbounds for Z2 in Z3:N;
                    a3 = 4*gamm*(Z2-Z3);
                    @inbounds for Z1 in Z2:N
                    a1 = bet*(N-Z1)*(Z1-Z5);
                    a2 = 4*gamm*(Z1-Z2);
                    tot = a1+a2+a3+a4+a5;
                    if Z1-Z5 == 0
                        final_size[Z5+1] = p_vec[w];
                    end
                    if a1 > 0
                       p_vec[w+1] = p_vec[w+1]+ p_vec[w]*a1/tot;
                    end
                    if a2 > 0
                        p_vec[w+N-Z2] = p_vec[w+N-Z2]+ p_vec[w]*a2/tot;
                    end
                    if a3 > 0
                        place3::Int64 = (N-Z3)*(N-Z3+1)/2;
                        p_vec[w+place3] = p_vec[w+place3]+ p_vec[w]*a3/tot;
                    end
                    if a4 > 0
                        place4::Int64 = (N-Z4)*(N-Z4+1)*(N-Z4+2)/6;
                        p_vec[w+place4] = p_vec[w+place4] + p_vec[w]*a4/tot;
                    end
                    if a5 > 0
                        p_vec[w] = p_vec[w]*a5/tot;
                    end
                    w = w + 1;
                end
            end
        end
    end
end
    return final_size
end

final_size = SI4R_fs(N,bet,gamm);

@time final_size = SI4R_fs(N,bet,gamm);

bar(0:N,final_size)

max_stars_repo_name: epimodels/epicookbook

max_issues_repo_name: epimodels/epicookbook

max_forks_repo_name: epimodels/epicookbook


text:
# include("FieldCalc.jl")
# # VecDist(X::Array) = [ X[2,1]-X[1,1], X[2,2]-X[1,2] , X[2,3]-X[1,3]]
# #X is the 2 coordinates to take distance between in format of [x₁,y₁,z₁;x₂,y₂,z₂]

"""
This function takes in an inner and outer radius and makes the ideal "D-shaped" toroid core

It returns a list of coordinates (r,z) where r is the radial location
on the toroid, and z is the height from the center axis. In a cross-section
view this is essentially x and y.

r₁ - inner radius
r₂ - outer radius
dr (kwarg) is the dr between each sampled location

"""
function DCoreGeom(r₁, r₂; dr = 0.0001,NPts = nothing,PlotOn=false,UpsamplePoints = 1e4)
    dZdr(r) = log.(sqrt(r₁ * r₂) ./ r) ./ sqrt.(log.(r ./ r₁) .* log.(r₂ ./ r))
    if ~(NPts===nothing)
        dr = (r₂ - r₁) / UpsamplePoints
    end
    rᵢ = r₁ + dr / 10#starting at a small fraction above zero, to avoid an inf.
    r = [rᵢ]
    zᵢ = 0.0
    z = [zᵢ] # initializing vector
    II = 1 #iterator
    NumIters = floor((r₂ - r₁) / dr)

    while (NumIters > II)
        dz = dZdr(rᵢ)
        zᵢ +=dz*dr
        rᵢ += dr

        push!(z, zᵢ)
        push!(r, rᵢ)
        II += 1
    end
    z .-= z[end] #the boundary condition at the end must be fixed to equal 0

    z[1] = 0 #because of the inf initial dz/dr, the initial boundary must be manually set to zero
    ZFlatHeight = z[2]
    StartPts = Int(ceil(ZFlatHeight/dr))
    for i in 1:StartPts
        insert!(z,1+i,z[1+i]*i/StartPts)
        insert!(r,1+i,r[1])
    end

    z = vcat(z[:],-1 .*reverse(z[2:end-2]))
    r = vcat(r[:],reverse(r[2:end-2]))
    DownsampleIndex = 1:Int(round(UpsamplePoints/NPts)):length(z)
    z = z[DownsampleIndex]
    r = r[DownsampleIndex]
    if PlotOn
        plot(r, z)
    end
    return r, z
end

function DCore_PointPath(r₁, r₂; NPts = 100)
    r, z = DCoreGeom(r₁, r₂; NPts = NPts*10,PlotOn=false)
    CumulativeDist = 0
    for i in 1:(length(r)-1)
        CumulativeDist += sqrt((r[i+1]-r[i])^2+(z[i+1]-z[i])^2)
    end
    TargetDist = CumulativeDist/(NPts+1)
    rᵣ = [r[1]]
    zᵣ = [z[1]]
    CumulativeDist2 = 0
    for i in 1:(length(r)-1)
        CumulativeDist2 += sqrt((r[i+1]-r[i])^2+(z[i+1]-z[i])^2)
        if CumulativeDist2>TargetDist
            push!(rᵣ,r[i])
            push!(zᵣ,z[i])
            CumulativeDist2 = 0
        end
    end

    PointPathZeros = zeros(size(rᵣ))
    PointPath = hcat(rᵣ,zᵣ,PointPathZeros)
    return PointPath
end

function FieldMapDToroid(r₁, r₂; NPts = 100, TestLayers = 20)

    PointPath = DCore_PointPath(r₁, r₂;NPts = NPts)
    FieldMapPointPath(PointPath, TestLayers; WeightRadius = true, InvWeights = true)
end

"""
RectCoords=[ 1   2.5 0
         5   2.5 0
         5  -2.5 0
         1  -2.5 0]

         Dr₁ = 1
         Dr₂ = 5

         Cr₁ = 2
         Cr₂ = 2
         CircCenter = [3,0,0]
         CompareD_Circ_Rect_Toroid(Dr₁, Dr₂,RectCoords,Cr₁, Cr₂,CircCenter,20)

"""
function CompareD_Circ_Rect_Toroid(Dr₁, Dr₂,RectCoords,Cr₁, Cr₂,CircCenter,TestLayers = 20)
    figure(75)
    FieldMapDToroid(Dr₁,Dr₂)

    figure(76)
    R_PointPath = MakeRectPointPath(RectCoords;NElement = 50, PlotOn=false)
    FieldMapPointPath(R_PointPath, 20; WeightRadius = true,InvWeights = true)

     figure(77)
    C_PointPath = MakeEllip(Cr₁,Cr₂;Center = CircCenter,NPts=100)
    C_PointPath = reverse(C_PointPath,dims=1)
FieldMapPointPath(C_PointPath, 20; WeightRadius = true,InvWeights = true)

end


text:
using ActinRingsMC

# System variables
ks = 7.82e-9
kd = 7.82e-9/5
T = 300
delta = 36e-9
Xc = 12e-9
EI = 3.6e-26
lf = 84
Nfil = 6
Nsca = 4

# Best to set lf directly and calculate the nearest Lf
Lf = calc_Lf(lf, delta)

# Calculate the lattice parameters
max_height = calc_max_lattice_height(Nsca, lf)
min_height = calc_min_lattice_height(Nsca, lf)

# This sets how many sites overlap between the filaments in the initial configuration
overlap = 1
start_height = max_height - Nsca*(overlap - 1)
start_radius = calc_radius(delta, start_height)
lattice = Lattice(start_height, max_height, min_height)

filaments = generate_starting_config(lattice, Nfil, Nsca, lf, overlap)

# This is required to make sure internal of lattice matches the starting configuration
update_occupancies!(filaments, lattice)

sysparms = SystemParams(ks, kd, T, delta, Xc, EI, Lf, lf, Nfil, Nsca)
system = System(sysparms, filaments, start_radius)

# Simulation parameters
iters = 2
steps = 1e4
max_bias_diff = 5
write_interval = 1e3
radius_move_freq = 0.5
filebase = "outs/test_run-0"
analytical_biases = true
read_biases = false
biases_filename = ""
restart_iter = 0
binwidth = 1

simparms = SimulationParams(
    iters,
    steps,
    max_bias_diff,
    write_interval,
    radius_move_freq,
    filebase,
    analytical_biases,
    read_biases,
    biases_filename,
    restart_iter,
    binwidth
)

mkpath("outs")

# First run
run_us!(system, lattice, simparms)

# Update simulation parameters for continuation run
filebase = "outs/test_run-1"
read_biases = true
biases_filename = "outs/test_run-0.biases"
restart_iter = 2
simparms = SimulationParams(
    iters,
    steps,
    max_bias_diff,
    write_interval,
    radius_move_freq,
    filebase,
    analytical_biases,
    read_biases,
    biases_filename,
    restart_iter,
    binwidth
)

# Continuation run
run_us!(system, lattice, simparms)

max_stars_repo_name: cumberworth/ActinRingsMC

max_issues_repo_name: cumberworth/ActinRingsMC

max_forks_repo_name: cumberworth/ActinRingsMC


text:
"""
    sparse_actions(pomcp::POMCPPlanner, pomdp::POMDPs.POMDP, h::BeliefNode, num_actions::Int)

Return an iterable object containing no more than `num_actions` actions to be considered at the current node.

Override this if you want to choose specific actions (you can override based on the POMDP type at the node level, or the belief type). If only a limited number of actions are to be considered, this function will be used to generate that set of actions. By default, it simply returns a random sampling of actions from the action space generated by `POMDPs.actions`.

If your problem has a continuous action space, you will want to override this to try a sensible set of action samples.
"""
function sparse_actions(pomcp::POMCPPlanner, pomdp::POMDPs.POMDP, h::BeliefNode, num_actions::Int)
    return sparse_actions(pomcp, pomdp, h.B, num_actions)
end

function sparse_actions{S,A,O}(pomcp::POMCPPlanner, pomdp::POMDPs.POMDP{S,A,O}, b::Any, num_actions::Int)
    if num_actions > 0
        as = POMDPs.actions(pomdp, b)
        if POMDPs.implemented(POMDPs.iterator, Tuple{typeof(as)})
            all_act = collect(POMDPs.iterator(as))
            selected_act = Array(A, min(num_actions, length(all_act)))
            len = length(selected_act)
            for i in 1:len
                j = rand(pomcp.solver.rng, 1:length(all_act))
                selected_act[i] = all_act[j]
                deleteat!(all_act, j)
            end
            return selected_act
        else
            return A[rand(pomcp.solver.rng, as) for j in 1:num_actions]
        end
    else
        return POMDPs.iterator(POMDPs.actions(pomdp, b))
    end
end

"""
Generate a new action when the set of actions is widened.
"""
function next_action(gen::RandomActionGenerator, mdp::POMDPs.POMDP, b, hnode::BeliefNode)
    rand(gen.rng, POMDPs.actions(mdp, b))
end
next_action(f::Function, pomdp::POMDPs.POMDP, b, hnode::BeliefNode) = f(pomdp, s, snode)


text:
using Test, GuessworkQuantumSideInfo, LinearAlgebra, Statistics, Random, UnPack
using GenericLinearAlgebra
using SCS, StableRNGs

using Pajarito, Cbc

TEST_MATLAB = false # requires MATLAB and MATLAB.jl installed
TEST_MISDP = true # requires Pajarito or another MISDP solver
TEST_BB84_MISDP = false # takes ~100 seconds; requires TEST_MISDP
TEST_MOI = false # incompatible with the current version of Pajarito

default_sdp_solver() = SCSSolver(verbose = 0, eps = 1e-6)

function misdp_solver(; verbose = false)
    sdp_solver = default_sdp_solver()
    mip_solver = CbcSolver(loglevel = 0)
    PajaritoSolver(
        cont_solver = sdp_solver,
        mip_solver = mip_solver,
        mip_solver_drives = false,
        use_mip_starts = true,
        solve_relax = false,
        log_level = verbose ? 3 : 0,
    )
end

include("../test_problems.jl")


text:
## test_plot.jl : test plotting commands

@testset "Random field plotting    " begin

    ## 1D ##
    m = Matern(0.75,1,p=1.)
    cov = CovarianceFunction(1,m)
    pts = 0:0.01:1
    grf1 = GaussianRandomField(cov,KarhunenLoeve(500),pts)
    plot(grf1)
    plot_eigenvalues(grf1)
    plot_eigenfunction(grf1, 3)
    plot_eigenfunction(grf1, 25)

    ## 2D ##
    e = Exponential(0.01)
    cov = CovarianceFunction(2,e)
    pts = 0:0.05:1
    grf2 = GaussianRandomField(cov,KarhunenLoeve(10),pts,pts)
    surface(grf2)
    heatmap(grf2)
    contourf(grf2)
    contour(grf2)
    plot_eigenvalues(grf2)
    plot_eigenfunction(grf1, 2)
    plot_eigenfunction(grf1, 10)

    ## Separable ##
    e1 = Exponential(0.1, p=1)
    sp1 = SeparableCovarianceFunction(e1)
    pts = 0:0.01:1
    sg1 = GaussianRandomField(sp1,KarhunenLoeve(10),pts)
    plot_eigenvalues(sg1)
    plot_eigenfunction(sg1, 1)
    plot_eigenfunction(sg1, 4)
    plot_eigenfunction(sg1, 7)

    e2 = Exponential(0.01, p=1)
    sp2 = SeparableCovarianceFunction(e1,e2)
    sg2 = GaussianRandomField(sp2,KarhunenLoeve(10),pts,pts)
    plot_eigenvalues(sg2)
    plot_eigenfunction(sg2, 2)
    plot_eigenfunction(sg2, 5)
    plot_eigenfunction(sg2, 8)

end


text:
function _transmission_pathway(pop::Population,
                               i::Int64, j::Int64)
  return pop.risks[[i; j], :x], pop.risks[[i; j], :y]
end

@recipe function f(tn::TransmissionNetwork,
                   pop::Population,
                   events::Events{T},
                   time::Float64;
                   show_individuals=true::Bool,
                   show_states=true::Bool) where T <: EpidemicModel
  if show_states & !show_individuals
    @warn "Will not show disease states without also showing individuals"
  end
  xguide --> ""
  yguide --> ""
  legend --> :none
  xlims --> extrema(pop.risks[!, :x]) .+ (sum(extrema(pop.risks[!, :x]).*(-1,1)) .* (-0.05, 0.05))
  ylims --> extrema(pop.risks[!, :y]) .+ (sum(extrema(pop.risks[!, :y]).*(-1,1)) .* (-0.05, 0.05))
  aspect_ratio := :equal
  linecolor --> :black
  axis --> nothing
  titlefontcolor --> :black
  framestyle --> :none
  n = pop.individuals
  susceptible = _ids_by_state(events, State_S, time)
  for i = 1:n
    if i ∉ susceptible
      source = findfirst(tn.internal[:, i])
      if source != nothing
        @series begin
          x, y = _transmission_pathway(pop, i, source)
          seriestype := :path
          label := ""
          x, y
        end
      end
    end
  end
  if show_individuals
    if show_states
      @series begin
        pop, events, time
      end
    else
      @series begin
        pop
      end
    end
  end
end

@recipe function f(tn::TransmissionNetwork,
                   pop::Population;
                   show_individuals=true::Bool) where T <: EpidemicModel
  xguide --> ""
  yguide --> ""
  legend --> :none
  xlims --> extrema(pop.risks[!, :x]) .+ (sum(extrema(pop.risks[!, :x]).*(-1,1)) .* (-0.05, 0.05))
  ylims --> extrema(pop.risks[!, :y]) .+ (sum(extrema(pop.risks[!, :y]).*(-1,1)) .* (-0.05, 0.05))
  aspect_ratio := :equal
  axis --> nothing
  titlefontcolor --> :black
  linecolor --> :black
  framestyle --> :none
  n = pop.individuals
  for i = 1:n
    if !tn.external[i] & any(tn.internal[:, i])
      @series begin
        source = findfirst(tn.internal[:, i])
        x, y = _transmission_pathway(pop, i, source)
        seriestype := :path
        label := ""
        x, y
      end
    end
  end
  if show_individuals
    @series begin
      pop
    end
  end
end


text:
export BrickDCSPITFPBaudrateConfig
struct BrickDCSPITFPBaudrateConfig
    enable_dynamic_baudrate::Bool
    minimum_dynamic_baudrate::Integer
end

export BrickDCSPITFPErrorCount
struct BrickDCSPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickDCProtocol1BrickletName
struct BrickDCProtocol1BrickletName
    protocol_version::Integer
    firmware_version::Vector{Integer}
    name::String
end

export BrickDCIdentity
struct BrickDCIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickDC
"""
Drives one brushed DC motor with up to 28V and 5A (peak)
"""
mutable struct BrickDC <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickDC(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.brick_dc")
        deviceInternal = package.BrickDC(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export set_velocity
"""
    $(SIGNATURES)

Sets the velocity of the motor. Whereas -32767 is full speed backward,
0 is stop and 32767 is full speed forward. Depending on the
acceleration (see :func:`Set Acceleration`), the motor is not immediately
brought to the velocity but smoothly accelerated.

The velocity describes the duty cycle of the PWM with which the motor is
controlled, e.g. a velocity of 3277 sets a PWM with a 10% duty cycle.
You can not only control the duty cycle of the PWM but also the frequency,
see :func:`Set PWM Frequency`.
"""
function set_velocity(device::BrickDC, velocity)
    device.deviceInternal.set_velocity(velocity)
end

export get_velocity
"""
    $(SIGNATURES)

Returns the velocity as set by :func:`Set Velocity`.
"""
function get_velocity(device::BrickDC)
    return device.deviceInternal.get_velocity()
end

export get_current_velocity
"""
    $(SIGNATURES)

Returns the *current* velocity of the motor. This value is different
from :func:`Get Velocity` whenever the motor is currently accelerating
to a goal set by :func:`Set Velocity`.
"""
function get_current_velocity(device::BrickDC)
    return device.deviceInternal.get_current_velocity()
end

export set_acceleration
"""
    $(SIGNATURES)

Sets the acceleration of the motor. It is given in *velocity/s*. An
acceleration of 10000 means, that every second the velocity is increased
by 10000 (or about 30% duty cycle).

For example: If the current velocity is 0 and you want to accelerate to a
velocity of 16000 (about 50% duty cycle) in 10 seconds, you should set
an acceleration of 1600.

If acceleration is set to 0, there is no speed ramping, i.e. a new velocity
is immediately given to the motor.
"""
function set_acceleration(device::BrickDC, acceleration)
    device.deviceInternal.set_acceleration(acceleration)
end

export get_acceleration
"""
    $(SIGNATURES)

Returns the acceleration as set by :func:`Set Acceleration`.
"""
function get_acceleration(device::BrickDC)
    return device.deviceInternal.get_acceleration()
end

export set_pwm_frequency
"""
    $(SIGNATURES)

Sets the frequency of the PWM with which the motor is driven.
Often a high frequency
is less noisy and the motor runs smoother. However, with a low frequency
there are less switches and therefore fewer switching losses. Also with
most motors lower frequencies enable higher torque.

If you have no idea what all this means, just ignore this function and use
the default frequency, it will very likely work fine.
"""
function set_pwm_frequency(device::BrickDC, frequency)
    device.deviceInternal.set_pwm_frequency(frequency)
end

export get_pwm_frequency
"""
    $(SIGNATURES)

Returns the PWM frequency as set by :func:`Set PWM Frequency`.
"""
function get_pwm_frequency(device::BrickDC)
    return device.deviceInternal.get_pwm_frequency()
end

export full_brake
"""
    $(SIGNATURES)

Executes an active full brake.

.. warning::
 This function is for emergency purposes,
 where an immediate brake is necessary. Depending on the current velocity and
 the strength of the motor, a full brake can be quite violent.

Call :func:`Set Velocity` with 0 if you just want to stop the motor.
"""
function full_brake(device::BrickDC)
    device.deviceInternal.full_brake()
end

export get_stack_input_voltage
"""
    $(SIGNATURES)

Returns the stack input voltage. The stack input voltage is the
voltage that is supplied via the stack, i.e. it is given by a
Step-Down or Step-Up Power Supply.
"""
function get_stack_input_voltage(device::BrickDC)
    return device.deviceInternal.get_stack_input_voltage()
end

export get_external_input_voltage
"""
    $(SIGNATURES)

Returns the external input voltage. The external input voltage is
given via the black power input connector on the DC Brick.

If there is an external input voltage and a stack input voltage, the motor
will be driven by the external input voltage. If there is only a stack
voltage present, the motor will be driven by this voltage.

.. warning::
 This means, if you have a high stack voltage and a low external voltage,
 the motor will be driven with the low external voltage. If you then remove
 the external connection, it will immediately be driven by the high
 stack voltage.
"""
function get_external_input_voltage(device::BrickDC)
    return device.deviceInternal.get_external_input_voltage()
end

export get_current_consumption
"""
    $(SIGNATURES)

Returns the current consumption of the motor.
"""
function get_current_consumption(device::BrickDC)
    return device.deviceInternal.get_current_consumption()
end

export enable
"""
    $(SIGNATURES)

Enables the driver chip. The driver parameters can be configured (velocity,
acceleration, etc) before it is enabled.
"""
function enable(device::BrickDC)
    device.deviceInternal.enable()
end

export disable
"""
    $(SIGNATURES)

Disables the driver chip. The configurations are kept (velocity,
acceleration, etc) but the motor is not driven until it is enabled again.

.. warning::
 Disabling the driver chip while the motor is still turning can damage the
 driver chip. The motor should be stopped calling :func:`Set Velocity` with 0
 before disabling the motor power. The :func:`Set Velocity` function will **not**
 wait until the motor is actually stopped. You have to explicitly wait for the
 appropriate time after calling the :func:`Set Velocity` function before calling
 the :func:`Disable` function.
"""
function disable(device::BrickDC)
    device.deviceInternal.disable()
end

export is_enabled
"""
    $(SIGNATURES)

Returns *true* if the driver chip is enabled, *false* otherwise.
"""
function is_enabled(device::BrickDC)
    return device.deviceInternal.is_enabled()
end

export set_minimum_voltage
"""
    $(SIGNATURES)

Sets the minimum voltage, below which the :cb:`Under Voltage` callback
is triggered. The minimum possible value that works with the DC Brick is 6V.
You can use this function to detect the discharge of a battery that is used
to drive the motor. If you have a fixed power supply, you likely do not need
this functionality.
"""
function set_minimum_voltage(device::BrickDC, voltage)
    device.deviceInternal.set_minimum_voltage(voltage)
end

export get_minimum_voltage
"""
    $(SIGNATURES)

Returns the minimum voltage as set by :func:`Set Minimum Voltage`
"""
function get_minimum_voltage(device::BrickDC)
    return device.deviceInternal.get_minimum_voltage()
end

export set_drive_mode
"""
    $(SIGNATURES)

Sets the drive mode. Possible modes are:

* 0 = Drive/Brake
* 1 = Drive/Coast

These modes are different kinds of motor controls.

In Drive/Brake mode, the motor is always either driving or braking. There
is no freewheeling. Advantages are: A more linear correlation between
PWM and velocity, more exact accelerations and the possibility to drive
with slower velocities.

In Drive/Coast mode, the motor is always either driving or freewheeling.
Advantages are: Less current consumption and less demands on the motor and
driver chip.
"""
function set_drive_mode(device::BrickDC, mode)
    device.deviceInternal.set_drive_mode(mode)
end

export get_drive_mode
"""
    $(SIGNATURES)

Returns the drive mode, as set by :func:`Set Drive Mode`.
"""
function get_drive_mode(device::BrickDC)
    return device.deviceInternal.get_drive_mode()
end

export set_current_velocity_period
"""
    $(SIGNATURES)

Sets a period with which the :cb:`Current Velocity` callback is triggered.
A period of 0 turns the callback off.
"""
function set_current_velocity_period(device::BrickDC, period)
    device.deviceInternal.set_current_velocity_period(period)
end

export get_current_velocity_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Current Velocity Period`.
"""
function get_current_velocity_period(device::BrickDC)
    return device.deviceInternal.get_current_velocity_period()
end

export set_spitfp_baudrate_config
"""
    $(SIGNATURES)

The SPITF protocol can be used with a dynamic baudrate. If the dynamic baudrate is
enabled, the Brick will try to adapt the baudrate for the communication
between Bricks and Bricklets according to the amount of data that is transferred.

The baudrate will be increased exponentially if lots of data is sent/received and
decreased linearly if little data is sent/received.

This lowers the baudrate in applications where little data is transferred (e.g.
a weather station) and increases the robustness. If there is lots of data to transfer
(e.g. Thermal Imaging Bricklet) it automatically increases the baudrate as needed.

In cases where some data has to transferred as fast as possible every few seconds
(e.g. RS485 Bricklet with a high baudrate but small payload) you may want to turn
the dynamic baudrate off to get the highest possible performance.

The maximum value of the baudrate can be set per port with the function
:func:`Set SPITFP Baudrate`. If the dynamic baudrate is disabled, the baudrate
as set by :func:`Set SPITFP Baudrate` will be used statically.

.. versionadded:: 2.3.5\$nbsp;(Firmware)
"""
function set_spitfp_baudrate_config(device::BrickDC, enable_dynamic_baudrate, minimum_dynamic_baudrate)
    device.deviceInternal.set_spitfp_baudrate_config(enable_dynamic_baudrate, minimum_dynamic_baudrate)
end

export get_spitfp_baudrate_config
"""
    $(SIGNATURES)

Returns the baudrate config, see :func:`Set SPITFP Baudrate Config`.

.. versionadded:: 2.3.5\$nbsp;(Firmware)
"""
function get_spitfp_baudrate_config(device::BrickDC)
    return device.deviceInternal.get_spitfp_baudrate_config()
end

export get_send_timeout_count
"""
    $(SIGNATURES)

Returns the timeout count for the different communication methods.

The methods 0-2 are available for all Bricks, 3-7 only for Master Bricks.

This function is mostly used for debugging during development, in normal operation
the counters should nearly always stay at 0.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function get_send_timeout_count(device::BrickDC, communication_method)
    return device.deviceInternal.get_send_timeout_count(communication_method)
end

export set_spitfp_baudrate
"""
    $(SIGNATURES)

Sets the baudrate for a specific Bricklet port.

If you want to increase the throughput of Bricklets you can increase
the baudrate. If you get a high error count because of high
interference (see :func:`Get SPITFP Error Count`) you can decrease the
baudrate.

If the dynamic baudrate feature is enabled, the baudrate set by this
function corresponds to the maximum baudrate (see :func:`Set SPITFP Baudrate Config`).

Regulatory testing is done with the default baudrate. If CE compatibility
or similar is necessary in your applications we recommend to not change
the baudrate.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function set_spitfp_baudrate(device::BrickDC, bricklet_port, baudrate)
    device.deviceInternal.set_spitfp_baudrate(bricklet_port, baudrate)
end

export get_spitfp_baudrate
"""
    $(SIGNATURES)

Returns the baudrate for a given Bricklet port, see :func:`Set SPITFP Baudrate`.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function get_spitfp_baudrate(device::BrickDC, bricklet_port)
    return device.deviceInternal.get_spitfp_baudrate(bricklet_port)
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Brick side. All
Bricklets have a similar function that returns the errors on the Bricklet side.

.. versionadded:: 2.3.3\$nbsp;(Firmware)
"""
function get_spitfp_error_count(device::BrickDC, bricklet_port)
    return device.deviceInternal.get_spitfp_error_count(bricklet_port)
end

export enable_status_led
"""
    $(SIGNATURES)

Enables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.

.. versionadded:: 2.3.1\$nbsp;(Firmware)
"""
function enable_status_led(device::BrickDC)
    device.deviceInternal.enable_status_led()
end

export disable_status_led
"""
    $(SIGNATURES)

Disables the status LED.

The status LED is the blue LED next to the USB connector. If enabled is is
on and it flickers if data is transfered. If disabled it is always off.

The default state is enabled.

.. versionadded:: 2.3.1\$nbsp;(Firmware)
"""
function disable_status_led(device::BrickDC)
    device.deviceInternal.disable_status_led()
end

export is_status_led_enabled
"""
    $(SIGNATURES)

Returns *true* if the status LED is enabled, *false* otherwise.

.. versionadded:: 2.3.1\$nbsp;(Firmware)
"""
function is_status_led_enabled(device::BrickDC)
    return device.deviceInternal.is_status_led_enabled()
end

export get_protocol1_bricklet_name
"""
    $(SIGNATURES)

Returns the firmware and protocol version and the name of the Bricklet for a
given port.

This functions sole purpose is to allow automatic flashing of v1.x.y Bricklet
plugins.
"""
function get_protocol1_bricklet_name(device::BrickDC, port)
    return device.deviceInternal.get_protocol1_bricklet_name(port)
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has an
accuracy of ±15%. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickDC)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Brick. Calling this function
on a Brick inside of a stack will reset the whole stack.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickDC)
    device.deviceInternal.reset()
end

export write_bricklet_plugin
"""
    $(SIGNATURES)

Writes 32 bytes of firmware to the bricklet attached at the given port.
The bytes are written to the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_bricklet_plugin(device::BrickDC, port, offset, chunk)
    device.deviceInternal.write_bricklet_plugin(port, offset, chunk)
end

export read_bricklet_plugin
"""
    $(SIGNATURES)

Reads 32 bytes of firmware from the bricklet attached at the given port.
The bytes are read starting at the position offset * 32.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function read_bricklet_plugin(device::BrickDC, port, offset)
    return device.deviceInternal.read_bricklet_plugin(port, offset)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Brick is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position is the position in the stack from '0' (bottom) to '8' (top).

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickDC)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickDC, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
A composite action during which a random user (called challenger)
select the worst candidate in the registry, pays a deposit
and call other users to vote if a item should be exluded from the registry.
If the mojority of users agree with the challenger, he receives a bonus an an item is removed.
Otherwise, the challenger loses his deposit and item remains in the registry.
"""
function challenge(registry, agents, deposit, voteFunc, redistributionFunc)
    if (length(registry) >= 10)
        challenger = agents[rand(1:end)]
        evaluations = evaluateCandidateByAgent.(registry, challenger)
        challengedIndex = indmin(evaluations);
        challengedValue = evaluations[challengedIndex]

        # if (min < mean(registry))
            challenger.balance -= deposit

            votingResult = voteFunc(registryMean(registry), challengedValue, agents)
            if (!votingResult[1])
                deleteat!(registry, challengedIndex)
            end
            redistributionFunc(votingResult, challenger, deposit)
        # end
    end
end

"""
A composite action during which a random user (called challenger)
select the worst candidate in the registry, pays a deposit
and call other users to vote if a item should be exluded from the registry.
If the mojority of users agree with the challenger, he receives a bonus an an item is removed.
Otherwise, the challenger loses his deposit and item remains in the registry.
"""
function application(registry, history, agents, voteFunc)
    candidate = getRegistryCandidate()
    push!(history, candidate)
    # println("Candidate: $candidate")
    if voteFunc(registryMean(registry), candidate, agents)[1]
        push!(registry, candidate);
    end
end


text:
################################################################################
function ∇loss( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace; kwargs...)

	predictions = unique([ s.z for branch ∈ branches for s ∈ branch if s.bif ]; atol=2*step(targets.parameter))
	λ = length(targets.targets)-length(predictions)

	if λ≠0
		Φ,∇Φ = measure(F,branches,θ,targets), ∇measure(F,branches,θ,targets)
		return errors(predictions,targets) - λ*log(Φ), ∇errors(F,predictions,θ,targets) - λ*∇Φ/Φ
	else
		return errors(predictions,targets), ∇errors(F,predictions,θ,targets)
	end
end

################################################################################
function ∇errors( F::Function, predictions::AbstractVector{<:BorderedArray}, θ::AbstractVector, targets::StateSpace)
	return mean( p′-> mean( z->(z.p-p′)^2, predictions; type=:geometric )*mean( z-> 2velocity(F,z,θ)/(z.p-p′), predictions; type=:arithmetic ), targets.targets; type=:arithmetic )
end

function ∇measure( F::Function, z::BorderedArray, θ::AbstractVector, targets::StateSpace; newtonOptions=NewtonPar(verbose=false,maxIter=800,tol=1e-6) )
	∂implicit, _, _ = newtonOptions.linsolver( -∂Fz(F,z,θ)', gradient(z->measure(F,z,θ,targets),z) )
	return gradient( θ -> measure(F,z,θ,targets) + F(z,θ)'∂implicit , θ ) + measure(F,z,θ,targets)*∇region(F,z,θ)
end

###########################################################################
function ∇measure( F::Function, branch::Branch, θ::AbstractVector, targets::StateSpace )
	return sum( s -> ∇measure(F,s.z,θ,targets)*s.ds, branch )
end

function ∇measure( F::Function, branches::AbstractVector{<:Branch}, θ::AbstractVector, targets::StateSpace )
	return sum( branch -> ∇measure(F,branch,θ,targets), branches )
end

############################################## gradient term due to changing integration region dz
deformation( F::Function, z::BorderedArray, θ::AbstractVector ) = -∂Fz(F,z,θ)\∂Fθ(F,z,θ) # todo@(gszep) matrix inverse is computational bottleneck
function ∇region( F::Function, z::BorderedArray, θ::AbstractVector )

	∇deformation = reshape( jacobian(z->deformation(F,z,θ),z), length(z),length(θ),length(z) )
	tangent = tangent_field(F,z,θ)
	t = [tangent.u;tangent.p]

	∇region = similar(θ)
	for k ∈ 1:length(θ)
		∇region[k] = t'∇deformation[:,k,:]t
	end

	return ∇region
end

########################################################################### jacobians

# statespace jacobian
∂Fu(F::Function,z::BorderedArray,θ::AbstractVector) = ∂Fz(F,z,θ)[:,Not(end)]

# parameter jacobian
∂Fθ(F::Function,z::BorderedArray,θ::AbstractVector) = jacobian( θ->F(z,θ), θ )

# augmented jacobian
∂Fz(F::Function,z::BorderedArray,θ::AbstractVector) = jacobian( z->F(z,θ), z )

############################################################# autodiff wrappers for BorderedArray
import ForwardDiff: derivative,gradient,jacobian
derivative(f,x::BorderedArray,s::BorderedArray) = derivative( λ -> f(x + λ*s), 0.0 )

function gradient( f, z::BorderedArray )
	return gradient( z -> f( BorderedArray(z[Not(end)],z[end]) ), [z.u; z.p] )
end

function jacobian( f, z::BorderedArray )
	return jacobian( z -> f( BorderedArray(z[Not(end)],z[end]) ), [z.u; z.p] )
end


text:
function hmatvec(A::HMat,v)
    assert(A.width == size(v,1));
    if A.blockType == LOWRANK
        return A.UMat*(A.VMat'*v);
    elseif A.blockType == DENSE
        return A.DMat*v;
    else
        uoffset = 0;
        u = zeros(A.height,size(v,2));
        for i = 1:4
            voffset = 0;
            for j = 1:4
                C = A.childHMat[i,j];
                u[uoffset+(1:C.height),:] += hmatvec(C,v[voffset+(1:C.width),:]);
                voffset += C.width;
                if j == 4
                    uoffset += C.height;
                end
            end
        end
    end
    return u;
end

function hmatTvec(A::HMat,v)
    assert(A.height == size(v,1));
    if A.blockType == LOWRANK
        return A.VMat*(A.UMat'*v);
    elseif A.blockType == DENSE
        return A.DMat'*v;
    else
        uoffset = 0;
        u = zeros(A.width,size(v,2));
        for j = 1:4
            voffset = 0;
            for i = 1:4
                C = A.childHMat[i,j];
                u[uoffset+(1:C.width),:] += hmatTvec(C,v[voffset+(1:C.height),:]);
                voffset += C.height;
                if i == 4
                    uoffset += C.width;
                end
            end
        end
    end
    return u;
end

function hvecmat(v,A::HMat)
    assert(A.height == size(v,2));
    if A.blockType == LOWRANK
        return (v*A.UMat)*A.VMat';
    elseif A.blockType == DENSE
        return v*A.DMat;
    else
        uoffset = 0;
        u = zeros(size(v,1),A.width);
        for j = 1:4
            voffset = 0;
            for i = 1:4
                C = A.childHMat[i,j];
                u[:,uoffset+(1:C.width)] += hvecmat(v[:,voffset+(1:C.height)],C);
                voffset += C.height;
                if i == 4
                    uoffset += C.width;
                end
            end
        end
    end
    return u;
end

function hvecmatT(v,A::HMat)
    assert(A.width == size(v,2));
    if A.blockType == LOWRANK
        return v*A.VMat*A.UMat';
    elseif A.blockType == DENSE
        return v*A.DMat';
    else
        uoffset = 0;
        u = zeros(size(v,1),A.height);
        for i = 1:4
            voffset = 0;
            for j = 1:4
                C = A.childHMat[i,j];
                u[:,uoffset+(1:C.height)] += hvecmatT(v[:,voffset+(1:C.width)],C);
                voffset += C.width;
                if j == 4
                    uoffset += C.height;
                end
            end
        end
    end
    return u;
end

max_stars_repo_name: YingzhouLi/HMat

max_issues_repo_name: YingzhouLi/HMat


text:
__precompile__()

module GetColumn

export getcolumn

using DataFrames
using TimeSeries

include("getcolumn.jl")

end # GetColumn


text:
function Hellothere()
  println("hello")
  x=1
  if x==1
    print("thuis is one")

    else
      println("so you are not one....")
  end
end
function get_fib()
  println("Hello world this is a sample program in julia");
end
Hellothere()

max_stars_repo_name: vbsteja/code

max_issues_repo_name: vbsteja/code

max_forks_repo_name: vbsteja/code


text:
@static if Sys.iswindows()
    using WinReg
end

function findSteamInstallDir()
    if Sys.iswindows()
        try
            # 64bit
            return querykey(WinReg.HKEY_LOCAL_MACHINE, "SOFTWARE\\WOW6432Node\\Valve\\Steam", "InstallPath")

        catch e
            try
                # 32bit
                return querykey(WinReg.HKEY_LOCAL_MACHINE, "SOFTWARE\\Valve\\Steam", "InstallPath")

            catch e
                # User doesn't have steam installed
                return nothing
            end
        end

    elseif Sys.isapple()
        return joinpath(homedir(), "Library", "Application Support", "Steam")

    elseif Sys.islinux()
        linuxSteamDirs = [
            joinpath(homedir(), ".local", "share", "Steam"),
            joinpath(homedir(), ".steam", "steam")
        ]

        for path in linuxSteamDirs
            if isdir(path)
                return path
            end
        end
    end
end

function cleanupPath(path)
    filename = basename(path)

    if lowercase(filename) == "celeste.exe"
        return dirname(path)
    end

    if Sys.isapple() && lowercase(filename) == "celeste.app"
        celesteContentDirs = [
            joinpath(path, "Contents", "Resources"), # Celeste 1.3.3.0 and newer
            joinpath(path, "Contents", "MacOS") # Before Celeste 1.3.3.0
        ]

        for path in celesteContentDirs
            if isdir(path)
                return path
            end
        end
    end

    return path
end

function findCelesteDir()
    steam = findSteamInstallDir()

    if steam !== nothing
        celesteDir = joinpath(steam, "steamapps", "common", "Celeste")

        if Sys.isapple()
            celesteDir = cleanupPath(joinpath(celesteDir, "Celeste.app"))
        end

        if isdir(celesteDir)
            return true, celesteDir
        end
    end

    return false, ""
end

function configureCelesteDir()
    if !haskey(config, "celeste_dir")
        found, celesteDir = findCelesteDir()
        macOS = Sys.isapple()
        target = if macOS "Celeste.app" else "Celeste.exe" end

        if !found || ask_dialog("Looks like you installed Celeste using Steam!\nAhorn depends on the Celeste installation directory to function; would you like Ahorn to use the Steam installation, or would you rather select a different installation of the game for it to use?", "Use Steam installation", "Manually select Celeste dir", window)
            info_dialog("Ahorn depends on the Celeste installation directory to function.\nPlease use the file dialog to select '$target' on your computer.", window)
            filters = if macOS tuple("*.app", "*.exe") else tuple("*.exe") end
            filename = openDialog("Select $target", window, filters)

            if filename == ""
                return false

            elseif lowercase(basename(filename)) != lowercase(target)
                warn_dialog("The file you selected is not named '$target'\nWill attempt to use this as the install directory anyway.", window)
            end

            celesteDir = cleanupPath(filename)
        end

        celesteGraphics = joinpath(celesteDir, "Content", "Graphics")
        celesteAtlases = joinpath(celesteGraphics, "Atlases")

        if !isdir(celesteGraphics) || !isdir(celesteAtlases)
            warn_dialog("The directory you selected does not contain the required paths.\nThe subdirectory Content/Graphics/ of a celeste installation with its contents is required for Ahorn to function.", window)
            exit()
        end

        config["celeste_dir"] = celesteDir
        saveConfig(config, true) # Manually save because it likely hasn't saved this to file
        info_dialog("Ahorn will now extract all files needed. This might take a while!\nPlease close this window to continue.", window)
    end

    return true
end

function needsExtraction(file::String, force::Bool=false)
    return !isfile(file) || force || filesize(file) == 0
end

function needsExtraction(from::String, to::String, force::Bool=false)
    return needsExtraction(to, force) || mtime(from) > mtime(to)
end

function extractGamedata(storage::String, force::Bool=false)
    celesteGraphics = joinpath(config["celeste_dir"], "Content", "Graphics")
    celesteAtlases = joinpath(celesteGraphics, "Atlases")

    storageXML = joinpath(storage, "XML")
    storageSprites = joinpath(storage, "Sprites")

    requiredFiles = Tuple{String, String}[
        (joinpath(celesteGraphics, "ForegroundTiles.xml"), joinpath(storageXML, "ForegroundTiles.xml")),
        (joinpath(celesteGraphics, "BackgroundTiles.xml"), joinpath(storageXML, "BackgroundTiles.xml")),
        (joinpath(celesteGraphics, "AnimatedTiles.xml"), joinpath(storageXML, "AnimatedTiles.xml")),
        (joinpath(celesteGraphics, "Sprites.xml"), joinpath(storageXML, "Sprites.xml")),
        (joinpath(celesteGraphics, "Portraits.xml"), joinpath(storageXML, "Portraits.xml")),

        (joinpath(celesteAtlases, "Gameplay.meta"), joinpath(storageSprites, "Gameplay.meta")),
        (joinpath(celesteAtlases, "Gui.meta"), joinpath(storageSprites, "Gui.meta")),
        (joinpath(celesteAtlases, "Portraits.meta"), joinpath(storageSprites, "Portraits.meta")),
        (joinpath(celesteAtlases, "Misc.meta"), joinpath(storageSprites, "Misc.meta"))
    ]

    requiredAtlases = Tuple{String, String}[
        (joinpath(celesteAtlases, "Gameplay0.data"), joinpath(storageSprites, "Gameplay.png")),
        (joinpath(celesteAtlases, "Gui0.data"), joinpath(storageSprites, "Gui.png"))
    ]

    needsInclude = any(Bool[needsExtraction(from, to, force) for (from, to) in requiredAtlases])

    if needsInclude
        include(Ahorn.abs"extract_sprites_images.jl")
    end

    for (from, to) in requiredAtlases
        if needsExtraction(from, to, force)
             # Making sure the method just loaded is used.
            Base.invokelatest(dumpSprites, from, to)
        end
    end

    for (from, to) in requiredFiles
        if needsExtraction(from, to, force)
            mkpath(dirname(to))

            cp(from, to, force=true)
        end
    end
end

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# CameraList.jl: interface to CameraList objects
"""
  CameraList(::System)

  Spinnaker SDK Camera List object constructor. Returns new Camera List object which
  enumerates available devices.
"""
mutable struct CameraList
  handle::spinCameraList

  function CameraList()
    hcamlist_ref = Ref(spinCameraList(C_NULL))
    spinCameraListCreateEmpty(hcamlist_ref)
    @assert hcamlist_ref[] != C_NULL
    camlist = new(hcamlist_ref[])
    _refresh!(camlist)
    finalizer(_release!, camlist)
    return camlist
  end

end

unsafe_convert(::Type{spinCameraList}, camlist::CameraList) = camlist.handle
unsafe_convert(::Type{Ptr{spinCameraList}}, camlist::CameraList) = pointer_from_objref(camlist)

# Clear list and release handle
function _release!(camlist::CameraList)
  spinCameraListClear(camlist)
  spinCameraListDestroy(camlist)
  camlist.handle = C_NULL
  return nothing
end

# Clear the list and reload enumerated cameras
function _refresh!(camlist::CameraList)
  spinCameraListClear(camlist)
  spinSystemGetCameras(spinsys, camlist)
end

"""
  length(::CameraList) -> Int

  Refresh the CameraList object and return the number of enumerated devices.
"""
function length(camlist::CameraList)
  _refresh!(camlist)
  nc = Ref(Csize_t(0))
  spinCameraListGetSize(camlist, nc)
  return Int(nc[])
end

"""
  show(::IO, ::CameraList)

  Write list of enumerated devices to supplied IO.
"""
function show(io::IO, camlist::CameraList)
  nc = length(camlist)  # This call will refresh the list
  write(io, "CameraList with $(nc[]) enumerated devices:\n")
  write(io, "ID\tSerial No.\tDescription\n")
  for i in 0:nc-1
    cam = camlist[i]
    vendorname = vendor(cam)
    modelname = model(cam)
    serialno = serial(cam)
    write(io, "$i\t$serialno\t$vendorname $modelname")
    _release!(cam) # Make sure this is cleaned up immediately
  end

end

"""
  getindex(::CameraList, ::Int) -> Camera

  Return camera by ID in specified CameraList. Note that IDs are zero-indexed.
"""
function getindex(camlist::CameraList, id::Int)
  nc = length(camlist)
  if (id >= 0) && (id < nc)
    hcam_ref = Ref(spinCamera(C_NULL))
    spinCameraListGet(camlist, id, hcam_ref)
    return Camera(hcam_ref[])
  else
    @error "Invalid camera index"
  end

end


text:
function recursiveCycle(param::MGparam,b::ArrayTypes,x::ArrayTypes,level::Int64)

# println(string("Starting level ", level));

gmresTol = 1e-5;

n = size(b,1);
nrhs = size(b,2);
numCores = param.numCores;
As = param.As;
nlevels = length(As);

if level==nlevels # This actually does not need to happen unless one level only is used (i.e. exact solver).
	x = solveCoarsest(param,b,x);
	return x;
end

Ps = param.Ps;
Rs = param.Rs;
AT = As[level];

const oneType = one(eltype(b));
const zeroType = zero(eltype(b));
r = param.memCycle[level].r;

r[:] = b;
if vecnorm(x)>0.0
   	SpMatMul(-oneType,AT,x,oneType,r,numCores)#  r -= A'*x;
end
D = param.relaxPrecs[level];
MM = identity;
Afun = identity;
if param.relaxType=="Jac-GMRES"
	y = param.memRelax[level].v_prec;
	MM(xx::ArrayTypes) = (SpMatMul(D,xx,y,numCores);return y;);
end

PT = Ps[level];
RT = Rs[level];

npresmth  = param.relaxPre(level);
npostsmth = param.relaxPost(level);

if param.relaxType=="Jac-GMRES"
    Afun = getAfun(AT,param.memRelax[level].Az,numCores)
	if nrhs == 1
		x = FGMRES_relaxation(Afun,r,x,npresmth,MM,gmresTol,false,true,numCores,param.memRelax[level])[1];
		# x = FGMRES(Afun,r,x,npresmth,MM,gmresTol,false,true,numCores,param.memRelax[level])[1];
	else
		for ii = 1:nrhs
			x[:,ii] = FGMRES_relaxation(Afun,r[:,ii],x[:,ii],npresmth,MM,gmresTol,false,true,numCores,param.memRelax[level])[1];
		end
		# x = BlockFGMRES(Afun,r,x,npresmth,MM,gmresTol,false,false,numCores, param.memRelax[level])[1];
		# x = BlockFGMRES_relaxation(AT,r,x,npresmth,MM,gmresTol,false,false,numCores, param.memRelax[level])[1];
	end
elseif param.relaxType == "VankaFaces"
	x = RelaxVankaFacesColor(AT,x,b,r,D,npresmth,numCores,param.Meshes[level],param.transferOperatorType=="SystemsFacesMixedLinear");
else
	x = relax(AT,r,x,b,D,npresmth,numCores);
end

SpMatMul(-oneType,AT,x,zeroType,r,numCores); #  r = -A'*x;
addVectors(oneType,b,r); # r = r + b;

xc = param.memCycle[level+1].x;
xc[:] = 0.0;
bc = param.memCycle[level+1].b;
bc = SpMatMul(RT,r,bc,numCores);
if level==nlevels-1
	# println("solving coarsest");
	xc = solveCoarsest(param,bc,xc);
else
	Ac = As[level+1];
    if param.cycleType == 'K'
		yzK = param.memKcycle[level].v_prec;
		AfunK = getAfun(Ac,param.memKcycle[level].Az,numCores);
		MMG(x) = (yzK[:] = 0.0; recursiveCycle(param,x,yzK,level+1)); # x does not change...
		if nrhs==1
			xc = FGMRES_relaxation(AfunK,bc,xc,2,MMG,gmresTol,false,true,numCores,param.memKcycle[level])[1];
			# xc = FGMRES(AfunK,bc,xc,2,MMG,gmresTol,false,true,numCores,param.memKcycle[level])[1];
		else
			# xc = BlockFGMRES_relaxation(Ac,bc,xc,2,MMG,gmresTol,false,true,numCores,param.memKcycle[level+1])[1];
			xc = BlockFGMRES(AfunK,bc,xc,2,MMG,gmresTol,false,true,numCores,param.memKcycle[level])[1];
		end
    else
		xc = recursiveCycle(param,bc,xc,level+1);# xc changes, bc does not change
		if param.cycleType=='W'
            xc = recursiveCycle(param,bc,xc,level+1);
		elseif param.cycleType=='F'
            param.cycleType='V';
            xc = recursiveCycle(param,bc,xc,level+1);
			param.cycleType='F';
        end
    end
end

# Base.A_mul_B!(oneType,P,xc,oneType,x); # x += P*xc;

x = SpMatMul(oneType,PT,xc,oneType,x,numCores); # x += PT'*xc;

r[:] = b;
SpMatMul(-oneType,AT,x,oneType,r,numCores); #  r -= A'*x;
if param.relaxType=="Jac-GMRES"
	Afun = getAfun(AT,param.memRelax[level].Az,numCores)
	if nrhs == 1
		x = FGMRES_relaxation(Afun,r,x,npostsmth,MM,gmresTol,false,true,numCores,param.memRelax[level])[1];
		# x = FGMRES(Afun,r,x,npostsmth,MM,gmresTol,false,true,numCores,param.memRelax[level])[1];
	else
		for ii=1:nrhs
			x[:,ii] = FGMRES_relaxation(Afun,r[:,ii],x[:,ii],npostsmth,MM,gmresTol,false,true,numCores,param.memRelax[level])[1];
		end
		# x = BlockFGMRES(Afun,r,x,npostsmth,MM,gmresTol,false,false,numCores, param.memRelax[level])[1];
		# x = BlockFGMRES_relaxation(AT,r,x,npostsmth,MM,gmresTol,false,false,numCores, param.memRelax[level])[1];
	end
elseif param.relaxType == "VankaFaces"
	x = RelaxVankaFacesColor(AT,x,b,r,D,npostsmth,numCores,param.Meshes[level],param.transferOperatorType=="SystemsFacesMixedLinear");
else
	x = relax(AT,r,x,b,D,npostsmth,numCores);
end

return x
end

function relax(AT::SparseMatrixCSC,r::ArrayTypes,x::ArrayTypes,b::ArrayTypes,D::SparseMatrixCSC,numit::Int64,numCores::Int64)
# x is actually the output... x and is being changed in the iterations.
# r does not end up accurate becasue we do not update it in the last iteration.
const oneType = one(eltype(r));
const zeroType = zero(eltype(r));
# nr0 = vecnorm(r);
# println(nr0)
for i=1:numit-1
	SpMatMul(oneType,D,r,oneType,x,numCores); # x = x + D'*r
	SpMatMul(-oneType,AT,x,zeroType,r,numCores) # r = -A'*x
	addVectors(oneType,b,r); # r = r + b;
	# println("Reduced: ", vecnorm(r)/nr0)
end
SpMatMul(oneType,D,r,oneType,x,numCores); # x = x + D'*r
return x
end

function solveCoarsest(param::MGparam,b::ArrayTypes,x::ArrayTypes,doTranspose::Int64=0)
if param.coarseSolveType == "MUMPS"
	applyMUMPS(param.LU,b,x,param.doTranspose);
elseif param.coarseSolveType == "BiCGSTAB"
	AT = param.As[end];
	maxIter = 100;
	tol = 0.1;
	out= -2;
	Afun = getAfun(AT,zeros(eltype(b),size(b)),param.numCores);
	D = param.LU;
	y = zeros(eltype(b),size(b));
	M1(xx::ArrayTypes) = (SpMatMul(D,xx,y,param.numCores);return y;);
	if size(b,2)==1
		b = vec(b);
		x, flag,rnorm,iter = KrylovMethods.bicgstb(Afun,b,tol = tol,maxIter = maxIter,M1 = M1,M2 = identity,out=out);
	else
		x, flag,rnorm,iter = KrylovMethods.blockBiCGSTB(Afun,b,tol = tol,maxIter = maxIter,M1 = M1,M2 = identity,out=out);
	end
	# println("Applied ",iter," BiCGSTAB");
elseif param.coarseSolveType == "GMRES"
	AT = param.As[end];
	maxIter = 1;
	tol = 0.01;
	out= -2;
	Afun = getAfun(AT,zeros(eltype(b),size(b)),param.numCores);
	D = param.LU;
	y = zeros(eltype(b),size(b));
	M2(xx::ArrayTypes) = (SpMatMul(D,xx,y,param.numCores);return y;);
	if size(b,2)==1
		b = vec(b);
		# x, flag,rnorm,iter = KrylovMethods.gmres(Afun,b,40,tol = tol,maxIter = maxIter,M = M2,out=out);
		x[:] = 0.0;
		(x,) = FGMRES(Afun,b,x,10,M2,tol,false,false,param.numCores);
	else
		error("Multiple RHS not supported");
	end
else
	x = param.LU\b;
end
return x;
end