text:
using ChainRulesCore
using KrylovKit
using LinearAlgebra
using Random
using Zygote

export num_grad

Zygote.@nograd StopFunction
Zygote.@nograd error
Zygote.@nograd FLint
Zygote.@nograd FRint
Zygote.@nograd BgFLint
Zygote.@nograd BgFRint
Zygote.@nograd leftorth
Zygote.@nograd rightorth
Zygote.@nograd ALCtoAC
Zygote.@nograd LRtoC
Zygote.@nograd initialA
Zygote.@nograd save
Zygote.@nograd load
Zygote.@nograd Random.seed!

@doc raw"
    num_grad(f, K::Real; [δ = 1e-5])

return the numerical gradient of `f` at `K` calculated with
`(f(K+δ/2) - f(K-δ/2))/δ`

# example

```jldoctest; setup = :(using TensorNetworkAD)
julia> TensorNetworkAD.num_grad(x -> x * x, 3) ≈ 6
true
```
"
function num_grad(f, K; δ::Real=1e-5)
    if eltype(K) == ComplexF64
        (f(K + δ / 2) - f(K - δ / 2)) / δ +
            (f(K + δ / 2 * 1.0im) - f(K - δ / 2 * 1.0im)) / δ * 1.0im
    else
        (f(K + δ / 2) - f(K - δ / 2)) / δ
    end
end

@doc raw"
    num_grad(f, K::AbstractArray; [δ = 1e-5])

return the numerical gradient of `f` for each element of `K`.

# example

```jldoctest; setup = :(using TensorNetworkAD, LinearAlgebra)
julia> TensorNetworkAD.num_grad(tr, rand(2,2)) ≈ I
true
```
"
function num_grad(f, a::AbstractArray; δ::Real=1e-5)
    b = Array(copy(a))
    df = map(CartesianIndices(b)) do i
        foo = x -> (ac = copy(b); ac[i] = x; f(_arraytype(a)(ac)))
        num_grad(foo, b[i], δ=δ)
    end
    return _arraytype(a)(df)
end

# patch since it's currently broken otherwise
function ChainRulesCore.rrule(::typeof(Base.typed_hvcat), ::Type{T}, rows::Tuple{Vararg{Int}}, xs::S...) where {T,S}
    y = Base.typed_hvcat(T, rows, xs...)
    function back(ȳ)
        return NoTangent(), NoTangent(), NoTangent(), permutedims(ȳ)...
    end
    return y, back
end

# improves performance compared to default implementation, also avoids errors
# with some complex arrays
function ChainRulesCore.rrule(::typeof(LinearAlgebra.norm), A::AbstractArray)
    n = norm(A)
    function back(Δ)
        return NoTangent(), Δ .* A ./ (n + eps(0f0)), NoTangent()
    end
    return n, back
end

function ChainRulesCore.rrule(::typeof(Base.sqrt), A::AbstractArray)
    As = Base.sqrt(A)
    function back(dAs)
        dA =  As' \ dAs ./2
        return NoTangent(), dA
    end
    return As, back
end

# adjoint for QR factorization
# https://journals.aps.org/prx/abstract/10.1103/PhysRevX.9.031041 eq.(5)
function ChainRulesCore.rrule(::typeof(qrpos), A::AbstractArray{T,2}) where {T}
    Q, R = qrpos(A)
    function back((dQ, dR))
        M = Array(R * dR' - dQ' * Q)
        dA = (UpperTriangular(R + I * 1e-12) \ (dQ + Q * _arraytype(Q)(Hermitian(M, :L)))' )'
        return NoTangent(), _arraytype(Q)(dA)
    end
    return (Q, R), back
end

function ChainRulesCore.rrule(::typeof(lqpos), A::AbstractArray{T,2}) where {T}
    L, Q = lqpos(A)
    function back((dL, dQ))
        M = Array(L' * dL - dQ * Q')
        dA = LowerTriangular(L + I * 1e-12)' \ (dQ + _arraytype(Q)(Hermitian(M, :L)) * Q)
        return NoTangent(), _arraytype(Q)(dA)
    end
    return (L, Q), back
end

"""
    dAMmap(Ai, Aip, Mi, L, R, j, J)

Aip means Aᵢ₊₁
```
               ┌──  Aᵢⱼ  ── ... ── AᵢJ   ──   ...  ──┐
               │     │              │                │
dMᵢJ    =     Lᵢⱼ ─ Mᵢⱼ  ── ... ──     ────── ...  ──Rᵢⱼ
               │     │              │                │
               └──  Aᵢ₊₁ⱼ  ─... ── Aᵢ₊₁J  ──  ...  ──┘

               ┌──  Aᵢⱼ  ── ... ──     ────── ...  ──┐
               │     │              │                │
dAᵢJ    =     Lᵢⱼ ─ Mᵢⱼ  ── ... ── dMᵢJ  ──── ...  ──Rᵢⱼ
               │     │              │                │
               └──  Aᵢ₊₁ⱼ  ─... ── Aᵢ₊₁J  ─── ...  ──┘

               ┌──  Aᵢⱼ  ── ... ── AᵢJ  ────  ...  ──┐        a ────┬──── c
               │     │              │                │        │     b     │
dAᵢ₊₁J   =     Lᵢⱼ ─ Mᵢⱼ  ── ... ── dMᵢJ ────  ...  ──Rᵢⱼ     ├─ d ─┼─ e ─┤
               │     │              │                │        │     g     │
               └──  Aᵢ₊₁ⱼ  ─... ──     ─────  ...  ──┘        f ────┴──── h

```
"""
function dAMmap(Ai, Aip, Mi, L, R, j, J; ifconj = false)
    Nj = size(Ai, 1)
        NL = (J - j + (J - j < 0) * Nj)
    NR = Nj - NL - 1
    L = copy(L)
    R = copy(R)
    for jj = 1:NL
        jr = j + jj - 1 - (j + jj - 1 > Nj) * Nj
        L = ein"((adf,abc),dgeb),fgh -> ceh"(L, Ai[jr], Mi[jr], Aip[jr])
    end
    for jj = 1:NR
        jr = j - jj + (j - jj < 1) * Nj
        R = ein"((ceh,abc),dgeb),fgh -> adf"(R, Ai[jr], Mi[jr], Aip[jr])
    end
    dAiJ = -ein"((adf,fgh),dgeb),ceh -> abc"(L, Aip[J], Mi[J], R)
    dAipJ = -ein"((adf,abc),dgeb),ceh -> fgh"(L, Ai[J], Mi[J], R)
    dMiJ = -ein"(adf,abc),(fgh,ceh) -> dgeb"(L, Ai[J], Aip[J], R)
    return conj(dAiJ), ifconj ? dAipJ : conj(dAipJ), conj(dMiJ)
end

function ChainRulesCore.rrule(::typeof(leftenv), ALu, ALd, M, FL; kwargs...)
    λL, FL = leftenv(ALu, ALd, M, FL)
    Ni, Nj = size(ALu)
    T = eltype(ALu[1,1])
    atype = _arraytype(M[1,1])
    function back((dλL, dFL))
        dALu = fill!(similar(ALu, atype), atype(zeros(T,size(ALu[1,1]))))
        dALd = fill!(similar(ALd, atype), atype(zeros(T,size(ALd[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        for j = 1:Nj, i = 1:Ni
            ir = i + 1 - Ni * (i == Ni)
            jr = j - 1 + Nj * (j == 1)
            dFL[i,j] -= Array(ein"abc,abc ->"(conj(FL[i,j]), dFL[i,j]))[] * FL[i,j]
            ξl, info = linsolve(FR -> FRmap(ALu[i,:], conj(ALd[ir,:]), M[i,:], FR, jr), conj(dFL[i,j]), -λL[i,j], 1; maxiter = 1)
            info.converged == 0 && @warn "ad's linsolve not converge"
            for J = 1:Nj
                dAiJ, dAipJ, dMiJ = dAMmap(ALu[i,:], conj(ALd[ir,:]), M[i,:], FL[i,j], ξl, j, J; ifconj = true)
                dALu[i,J] += dAiJ
                dALd[ir,J] += dAipJ
                dM[i,J] += dMiJ
            end
        end
        return NoTangent(), dALu, dALd, dM, NoTangent()
    end
    return (λL, FL), back
end

function ChainRulesCore.rrule(::typeof(rightenv), ARu, ARd, M, FR; kwargs...)
    λR, FR = rightenv(ARu, ARd, M, FR)
    Ni, Nj = size(ARu)
    T = eltype(ARu[1,1])
    atype = _arraytype(M[1,1])
    function back((dλ, dFR))
        dARu = fill!(similar(ARu, atype), atype(zeros(T,size(ARu[1,1]))))
        dARd = fill!(similar(ARd, atype), atype(zeros(T,size(ARd[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        for j = 1:Nj, i = 1:Ni
            ir = i + 1 - Ni * (i == Ni)
            jr = j - 1 + Nj * (j == 1)
            dFR[i,jr] -= Array(ein"abc,abc ->"(conj(FR[i,jr]), dFR[i,jr]))[] * FR[i,jr]
            ξr, info = linsolve(FL -> FLmap(ARu[i,:], conj(ARd[ir,:]), M[i,:], FL, j), conj(dFR[i,jr]), -λR[i,jr], 1; maxiter = 1)
            info.converged == 0 && @warn "ad's linsolve not converge"
            for J = 1:Nj
                dAiJ, dAipJ, dMiJ = dAMmap(ARu[i,:], conj(ARd[ir,:]), M[i,:], ξr, FR[i,jr], j, J; ifconj = true)
                dARu[i,J] += dAiJ
                dARd[ir,J] += dAipJ
                dM[i,J] += dMiJ
            end
        end
        return NoTangent(), dARu, dARd, dM, NoTangent()
    end
    return (λR, FR), back
end

"""
    ACdmap(ACij, FLj, FRj, Mj, II)

```
.        .         .
.        .         .
.        .         .
│        │         │          a ────┬──── c
FLᵢ₋₁ⱼ ─ Mᵢ₋₁ⱼ ──  FRᵢ₋₁ⱼ     │     b     │
│        │         │          ├─ d ─┼─ e ─┤
FLᵢⱼ ─── Mᵢⱼ ───── FRᵢⱼ       │     g     │
│        │         │          f ────┴──── h
└─────── ACᵢⱼ ─────┘
```
"""
function ACdmap(ACij, FLj, FRj, Mj, II)
    Ni = size(FLj,1)
    ACdm = copy(ACij)
    for i=1:Ni
        ir = II-(i-1) + (II-(i-1) < 1)*Ni
        ACdm = ein"((adf,fgh),dgeb),ceh -> abc"(FLj[ir],ACdm,Mj[ir],FRj[ir])
    end
    return ACdm
end

"""
    ACdFMmap(FLj, Mi, FRj, AC, ACd, i, II)

```
               ┌─────  ACᵢⱼ ─────┐
               │        │        │
              FLᵢⱼ ─── Mᵢⱼ ──── FRᵢⱼ
               │        │        │
               ⋮         ⋮        ⋮
               │        │        │
dMIⱼ    =     FLIⱼ ───     ──── FRIⱼ
               │        │        │
               ⋮         ⋮        ⋮
               │        │        │
               └─────  ACdᵢ₋₁ⱼ ──┘

               ┌─────  ACᵢⱼ ─────┐
               │        │        │
              FLᵢⱼ ─── Mᵢⱼ  ─── FRᵢⱼ
               │        │        │
               ⋮         ⋮        ⋮
               │        │        │
dFLIⱼ   =        ───── MIⱼ ───  FRIⱼ
               │        │        │
               ⋮         ⋮        ⋮
               │        │        │
               └─────  ACdᵢ₋₁ⱼ ──┘

               ┌─────  ACᵢⱼ ─────┐
               │        │        │
              FLᵢⱼ ─── Mᵢⱼ  ─── FRᵢⱼ
               │        │        │
               ⋮         ⋮        ⋮
               │        │        │
dFRIⱼ   =     FLIⱼ ─── MIⱼ ─────       a ────┬──── c
               │        │        │     │     b     │
               ⋮         ⋮        ⋮     ├─ d ─┼─ e ─┤
               │        │        │     │     g     │
               └─────  ACdᵢ₋₁ⱼ ──┘     f ────┴──── h
```
"""
function ACdFMmap(FLj, Mj, FRj, AC, ACd, i, II)
    Ni = size(FLj, 1)
    Nu = (II - i + (II - i < 0) * Ni)
    Nd = Ni - Nu - 1
    AC = copy(AC)
    ACd = copy(ACd)
    for ii = 1:Nu
        ir = i + ii - 1 - (i + ii - 1 > Ni) * Ni
        AC = ein"((adf,abc),dgeb),ceh -> fgh"(FLj[ir], AC, Mj[ir], FRj[ir])
    end
    for ii = 1:Nd
        ir = i - ii + (i - ii < 1) * Ni
        ACd = ein"((adf,fgh),dgeb),ceh -> abc"(FLj[ir], ACd, Mj[ir], FRj[ir])
    end
    dFLIj = -ein"((abc,ceh),dgeb),fgh -> adf"(AC, FRj[II], Mj[II], ACd)
    dMIj = -ein"(adf,abc),(fgh,ceh) -> dgeb"(FLj[II], AC, ACd, FRj[II])
    dFRIj = -ein"((abc,adf),dgeb),fgh -> ceh"(AC, FLj[II], Mj[II], ACd)
    return conj(dFLIj), conj(dMIj), conj(dFRIj)
end

function ChainRulesCore.rrule(::typeof(ACenv), AC, FL, M, FR; kwargs...)
    λAC, AC = ACenv(AC, FL, M, FR)
    Ni, Nj = size(AC)
    T = eltype(AC[1,1])
    atype = _arraytype(M[1,1])
    function back((dλ, dAC))
        dFL = fill!(similar(FL, atype), atype(zeros(T,size(FL[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        dFR = fill!(similar(FR, atype), atype(zeros(T,size(FR[1,1]))))
        for j = 1:Nj, i = 1:Ni
            if dAC[i,j] !== nothing
                ir = i - 1 + Ni * (i == 1)
                dAC[i,j] -= Array(ein"abc,abc ->"(conj(AC[i,j]), dAC[i,j]))[] * AC[i,j]
                ξAC, info = linsolve(ACd -> ACdmap(ACd, FL[:,j], FR[:,j], M[:,j], ir), conj(dAC[i,j]), -λAC[i,j], 1; maxiter = 1)
                info.converged == 0 && @warn "ad's linsolve not converge"
                # errAC = ein"abc,abc ->"(AC[i,j], ξAC)[]
                # abs(errAC) > 1e-1 && throw("AC and ξ aren't orthometric. $(errAC) $(info)")
                # @show info ein"abc,abc ->"(AC[i,j], ξAC)[] ein"abc,abc -> "(AC[i,j], dAC[i,j])[]
                for II = 1:Ni
                    dFLIj, dMIj, dFRIj = ACdFMmap(FL[:,j], M[:,j], FR[:,j], AC[i,j], ξAC, i, II)
                    dFL[II,j] += dFLIj
                    dM[II,j] += dMIj
                    dFR[II,j] += dFRIj
                end
            end
        end
        return NoTangent(), NoTangent(), dFL, dM, dFR
    end
    return (λAC, AC), back
end

"""
    Cdmap(Cij, FLj, FRj, II)

```
.                .
.                .
.                .
│                │
FLᵢ₋₁ⱼ₊₁ ─────  FRᵢ₋₁ⱼ      a ─── b
│                │          │     │
FLᵢⱼ₊₁ ───────  FRᵢⱼ        ├─ c ─┤
│                │          │     │
└────── Cᵢⱼ ─────┘          d ─── e
```
"""
function Cdmap(Cij, FLjp, FRj, II)
    Ni = size(FLjp,1)
    Cdm = copy(Cij)
    for i=1:Ni
        ir = II-(i-1) + (II-(i-1) < 1)*Ni
        Cdm = ein"(acd,de),bce -> ab"(FLjp[ir],Cdm,FRj[ir])
    end
    return Cdm
end

"""
    CdFMmap(FLj, FRj, C, Cd, i, II)

```
               ┌────  Cᵢⱼ ────┐
               │              │
              FLᵢⱼ₊₁───────  FRᵢⱼ
               │              │
               ⋮               ⋮
               │              │
dFLIⱼ₊₁ =        ──────────  FRIⱼ
               │              │
               ⋮               ⋮
               │              │
               └──── Cdᵢⱼ ────┘

               ┌────  Cᵢⱼ ────┐
               │              │
              FLᵢⱼ₊₁ ──────  FRᵢⱼ
               │              │
               ⋮               ⋮
               │              │
dFRIⱼ   =     FLᵢ₊Iⱼ₊₁ ──────           a ─── b
               │              │         │     │
               ⋮               ⋮         ├─ c ─┤
               │              │         │     │
               └──── Cdᵢⱼ ────┘         d ─── e
```
"""
function CdFMmap(FLjp, FRj, C, Cd, i, II)
    Ni = size(FLjp, 1)
    Nu = (II - i + (II - i < 0) * Ni)
    Nd = Ni - Nu - 1
    C = copy(C)
    Cd = copy(Cd)
    for ii = 1:Nu
        ir = i + ii - 1 - (i + ii - 1 > Ni) * Ni
        C = ein"(acd,ab),bce -> de"(FLjp[ir], C, FRj[ir])
    end
    for ii = 1:Nd
        ir = i - ii + (i - ii < 1) * Ni
        Cd = ein"(acd,de),bce -> ab"(FLjp[ir], Cd, FRj[ir])
    end
    dFLIjp = -ein"(ab,bce),de -> acd"(C, FRj[II], Cd)
    dFRIj = -ein"(ab,acd),de -> bce"(C, FLjp[II], Cd)
    return conj(dFLIjp), conj(dFRIj)
end

function ChainRulesCore.rrule(::typeof(Cenv), C, FL, FR; kwargs...)
    λC, C = Cenv(C, FL, FR)
    Ni, Nj = size(C)
    T = eltype(C[1,1])
    atype = _arraytype(FL[1,1])
    function back((dλ, dC))
        dFL = fill!(similar(FL, atype), atype(zeros(T,size(FL[1,1]))))
        dFR = fill!(similar(FR, atype), atype(zeros(T,size(FR[1,1]))))
        for j = 1:Nj, i = 1:Ni
            if dC[i,j] !== nothing
                ir = i - 1 + Ni * (i == 1)
                jr = j + 1 - (j==Nj) * Nj
                dC[i,j] -= Array(ein"ab,ab ->"(conj(C[i,j]), dC[i,j]))[] * C[i,j]
                ξC, info = linsolve(Cd -> Cdmap(Cd, FL[:,jr], FR[:,j], ir), conj(dC[i,j]), -λC[i,j], 1; maxiter = 1)
                info.converged == 0 && @warn "ad's linsolve not converge"
                # errC = ein"ab,ab ->"(C[i,j], ξC)[]
                # abs(errC) > 1e-1 && throw("C and ξ aren't orthometric. $(errC) $(info)")
                # @show info ein"ab,ab ->"(C[i,j], ξC)[] ein"ab,ab -> "(C[i,j], dC[i,j])[]
                for II = 1:Ni
                    dFLIjp, dFRIj = CdFMmap(FL[:,jr], FR[:,j], C[i,j], ξC, i, II)
                    dFL[II,jr] += dFLIjp
                    dFR[II,j] += dFRIj
                end
            end
        end
        return NoTangent(), NoTangent(), dFL, dFR
    end
    return (λC, C), back
end

function ChainRulesCore.rrule(::typeof(obs_FL), ALu, ALd, M, FL; kwargs...)
    λL, FL = obs_FL(ALu, ALd, M, FL)
    Ni, Nj = size(ALu)
    T = eltype(ALu[1,1])
    atype = _arraytype(M[1,1])
    function back((dλL, dFL))
        dALu = fill!(similar(ALu, atype), atype(zeros(T,size(ALu[1,1]))))
        dALd = fill!(similar(ALd, atype), atype(zeros(T,size(ALd[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        for j = 1:Nj, i = 1:Ni
            ir = Ni + 1 - i
            jr = j - 1 + Nj * (j == 1)
            dFL[i,j] -= Array(ein"abc,abc ->"(conj(FL[i,j]), dFL[i,j]))[] * FL[i,j]
            ξl, info = linsolve(FR -> FRmap(ALu[i,:], ALd[ir,:], M[i,:], FR, jr), conj(dFL[i,j]), -λL[i,j], 1; maxiter = 1)
            info.converged == 0 && @warn "ad's linsolve not converge"
            for J = 1:Nj
                dAiJ, dAipJ, dMiJ = dAMmap(ALu[i,:], ALd[ir,:], M[i,:], FL[i,j], ξl, j, J)
                dALu[i,J] += dAiJ
                dALd[ir,J] += dAipJ
                dM[i,J] += dMiJ
            end
        end
        return NoTangent(), dALu, dALd, dM, NoTangent()
    end
    return (λL, FL), back
end

function ChainRulesCore.rrule(::typeof(obs_FR), ARu, ARd, M, FR; kwargs...)
    λR, FR = obs_FR(ARu, ARd, M, FR)
    Ni, Nj = size(ARu)
    T = eltype(ARu[1,1])
    atype = _arraytype(M[1,1])
    function back((dλ, dFR))
        dARu = fill!(similar(ARu, atype), atype(zeros(T,size(ARu[1,1]))))
        dARd = fill!(similar(ARd, atype), atype(zeros(T,size(ARd[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        for j = 1:Nj, i = 1:Ni
            ir = Ni + 1 - i
            jr = j - 1 + Nj * (j == 1)
            dFR[i,jr] -= Array(ein"abc,abc ->"(conj(FR[i,jr]), dFR[i,jr]))[] * FR[i,jr]
            ξr, info = linsolve(FL -> FLmap(ARu[i,:], ARd[ir,:], M[i,:], FL, j), conj(dFR[i,jr]), -λR[i,jr], 1; maxiter = 1)
            info.converged == 0 && @warn "ad's linsolve not converge"
            for J = 1:Nj
                dAiJ, dAipJ, dMiJ = dAMmap(ARu[i,:], ARd[ir,:], M[i,:], ξr, FR[i,jr], j, J)
                dARu[i,J] += dAiJ
                dARd[ir,J] += dAipJ
                dM[i,J] += dMiJ
            end
        end
        return NoTangent(), dARu, dARd, dM, NoTangent()
    end
    return (λR, FR), back
end

"""
    dBgAMmap(Ai, Aip, Mi, Mip, L, R, j, J)
Aip means Aᵢ₊₁
```
               ┌──  Aᵢⱼ  ── ... ── AᵢJ   ──  ...   ──┐
               │     │              │                │
               │─   Mᵢⱼ  ── ...  ──   ────── ...   ──│
dMᵢJ    =     Lᵢⱼ    │              │                Rᵢⱼ
               │─   Mᵢ₊₁ⱼ ──... ──Mᵢ₊₁J ───  ...   ──│
               │     │              │                │
               └──  Aᵢₚⱼ  ─  ...── AᵢₚJ  ──── ...   ──┘
               ┌──  Aᵢⱼ  ── ... ── AᵢJ   ──  ...   ──┐
               │     │              │                │
               │─   Mᵢⱼ  ── ...  ──MᵢJ   ─── ...   ──│
dMᵢₚJ    =    Lᵢⱼ    │              │               Rᵢⱼ
               │─   Mᵢ₊₁ⱼ ──... ──     ───   ...   ──│
               │     │              │                │
               └──  Aᵢₚⱼ  ─  ...── AᵢₚJ  ──── ...   ──┘
               ┌──  Aᵢⱼ  ── ... ──       ──  ...   ──┐
               │     │              │                │
               │─   Mᵢⱼ  ── ... ── MᵢJ ────  ...   ──│
dAᵢJ    =     Lᵢⱼ    │              │                Rᵢⱼ
               │─   Mᵢ₊₁ⱼ ──... ── Mᵢ₊₁J ─── ...   ──│
               │     │              │                │
               └──  Aᵢₚⱼ  ─  ...── AᵢₚJ  ──── ...   ──┘
               ┌──  Aᵢⱼ  ── ... ── AᵢJ   ──  ...   ──┐      a ────┬──── c
               │     │              │                │      │     b     │
               │─   Mᵢⱼ  ── ... ── MᵢJ ────  ...   ──│      ├─ d ─┼─ e ─┤
dAᵢₚJ    =     Lᵢⱼ    │              │               Rᵢⱼ    │     f     │
               │─   Mᵢ₊₁ⱼ ──... ── Mᵢ₊₁J ─── ...   ──│      ├─ g ─┼─ h ─┤
               │     │              │                │      │     j     │
               └──  Aᵢₚⱼ  ─ ... ──      ──── ...   ──┘      i ────┴──── k
```
"""
function dBgAMmap(Ai, Aip, Mi, Mip, L, R, j, J)
    Nj = size(Ai, 1)
        NL = (J - j + (J - j < 0) * Nj)
    NR = Nj - NL - 1
    L = copy(L)
    R = copy(R)
    for jj = 1:NL
        jr = j + jj - 1 - (j + jj - 1 > Nj) * Nj
        L = ein"(((adgi,abc),dfeb),gjhf),ijk -> cehk"(L, Ai[jr], Mi[jr], Mip[jr], Aip[jr])
    end
    for jj = 1:NR
        jr = j - jj + (j - jj < 1) * Nj
        R = ein"(((cehk,abc),dfeb),gjhf),ijk -> adgi"(R, Ai[jr], Mi[jr], Mip[jr], Aip[jr])
    end
    dAiJ = -ein"(((adgi,ijk),gjhf),dfeb), cehk -> abc"(L, Aip[J], Mip[J], Mi[J], R)
    dAipJ = -ein"(((adgi,abc),dfeb),gjhf), cehk -> ijk"(L, Ai[J], Mi[J], Mip[J], R)
    dMiJ = -ein"(adgi,abc),(gjhf,(ijk,cehk)) -> dfeb"(L, Ai[J], Mip[J], Aip[J], R)
    dMipJ = -ein"((adgi,abc),dfeb),(ijk,cehk)-> gjhf"(L, Ai[J], Mi[J], Aip[J], R)
    return conj(dAiJ), conj(dAipJ), conj(dMiJ), conj(dMipJ)
end

function ChainRulesCore.rrule(::typeof(bigleftenv), ALu, ALd, M, BgFL; kwargs...)
    λL, BgFL = bigleftenv(ALu, ALd, M, BgFL)
    Ni, Nj = size(ALu)
    T = eltype(ALu[1,1])
    atype = _arraytype(M[1,1])
    function back((dλL, dBgFL))
        dALu = fill!(similar(ALu, atype), atype(zeros(T,size(ALu[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        dALd = fill!(similar(ALd, atype), atype(zeros(T,size(ALd[1,1]))))
        for j = 1:Nj, i = 1:Ni
            if dBgFL[i,j] !== nothing
                ir = i + 1 - Ni * (i==Ni)
                irr = i + 2 - Ni * (i + 2 > Ni)
                jr = j - 1 + Nj * (j == 1)
                dBgFL[i,j] -= Array(ein"abcd,abcd ->"(conj(BgFL[i,j]), dBgFL[i,j]))[] * BgFL[i,j]
                ξl, info = linsolve(BgFR -> BgFRmap(ALu[i,:], ALd[irr,:], M[i,:], M[ir,:], BgFR, jr), conj(dBgFL[i,j]), -λL[i,j], 1; maxiter = 1)
                info.converged == 0 && @warn "ad's linsolve not converge"
                # errL = ein"abc,cba ->"(FL[i,j], ξl)[]
                # abs(errL) > 1e-1 && throw("FL and ξl aren't orthometric. $(errL) $(info)")
                # @show info ein"abc,cba ->"(FL[i,j], ξl)[] ein"abc,abc -> "(FL[i,j], dFL[i,j])[]
                for J = 1:Nj
                    dAiJ, dAipJ, dMiJ, dMipJ = dBgAMmap(ALu[i,:], ALd[irr,:], M[i,:], M[ir,:], BgFL[i,j], ξl, j, J)
                    dALu[i,J] += dAiJ
                    dALd[i,J] += dAipJ
                    dM[i,J] += dMiJ
                    dM[ir,J] += dMipJ
                end
            end
        end
        return NoTangent(), dALu, dALd, dM, NoTangent()
    end
    return (λL, BgFL), back
end

function ChainRulesCore.rrule(::typeof(bigrightenv), ARu, ARd, M, BgFR; kwargs...)
    λR, BgFR = bigrightenv(ARu, ARd, M, BgFR)
    Ni, Nj = size(ARu)
    T = eltype(ARu[1,1])
    atype = _arraytype(M[1,1])
    function back((dλ, dBgFR))
        dARu = fill!(similar(ARu, atype), atype(zeros(T,size(ARu[1,1]))))
        dM = fill!(similar(M, atype), atype(zeros(T,size(M[1,1]))))
        dARd = fill!(similar(ARd, atype), atype(zeros(T,size(ARd[1,1]))))
        for j = 1:Nj, i = 1:Ni
            ir = i + 1 - Ni * (i==Ni)
            irr = i + 2 - Ni * (i + 2 > Ni)
            jr = j - 1 + Nj * (j == 1)
            if dBgFR[i,jr] !== nothing
                dBgFR[i,jr] -= Array(ein"abcd,abcd ->"(conj(BgFR[i,jr]), dBgFR[i,jr]))[] * BgFR[i,jr]
                ξr, info = linsolve(BgFL -> BgFLmap(ARu[i,:], ARd[irr,:], M[i,:], M[ir,:], BgFL, j), conj(dBgFR[i,jr]), -λR[i,jr], 1; maxiter = 1)
                info.converged == 0 && @warn "ad's linsolve not converge"
                # errR = ein"abc,cba ->"(ξr, FR[i,jr])[]
                # abs(errR) > 1e-1 && throw("FR and ξr aren't orthometric. $(errR) $(info)")
                # @show info ein"abc,cba ->"(ξr, FR[i,jr])[] ein"abc,abc -> "(FR[i,jr], dFR[i,jr])[]
                for J = 1:Nj
                    dAiJ, dAipJ, dMiJ, dMipJ = dBgAMmap(ARu[i,:], ARd[irr,:], M[i,:], M[ir,:], ξr, BgFR[i,jr], j, J)
                    dARu[i,J] += dAiJ
                    dARd[i,J] += dAipJ
                    dM[i,J] += dMiJ
                    dM[ir,J] += dMipJ
                end
            end
        end
        return NoTangent(), dARu, dARd, dM, NoTangent()
    end
    return (λR, BgFR), back
end


text:
export Signature, GitRepo, GitObject, GitAnyObject, GitBlob, GitCommit, GitTag,
       GitTree, GitReference, GitBranch, GitRemote, Sym

typealias GitOffT Int64
typealias GitTimeT Int64

typealias MaybeDict Union(Nothing, Dict{Any,Any})
typealias MaybeString Union(Nothing, String)

# time in a signature

immutable TimeStruct
    time::GitTimeT # time in seconds from epoch
    offset::Cint   # timezone offset in minutes
end
TimeStruct() = TimeStruct(zero(GitTimeT), zero(Cint))

# an action signature (committers, taggers, etc)
immutable SignatureStruct
    name::Ptr{Uint8}  # full name of the author
    email::Ptr{Uint8} # email of the author
    when::TimeStruct  # time when the action happened
end
SignatureStruct() = SignatureStruct(zero(Ptr{Uint8}),
                                    zero(Ptr{Uint8}),
                                    TimeStruct())
immutable BufferStruct
    ptr::Ptr{Uint8}
    asize::Csize_t
    size::Csize_t
end
BufferStruct() = BufferStruct(C_NULL, 0, 0)

Base.bytestring(b::BufferStruct) = bytestring(b.ptr, b.size)

immutable TransferProgressStruct
    total_objects::Cuint
    indexed_objects::Cuint
    received_objects::Cuint
    local_objects::Cuint
    total_deltas::Cuint
    indexed_deltas::Cuint
    received_bytes::Csize_t
end
TransferProgressStruct() = TransferProgressStruct(zero(Cuint),
                                                  zero(Cuint),
                                                  zero(Cuint),
                                                  zero(Cuint),
                                                  zero(Cuint),
                                                  zero(Cuint),
                                                  zero(Csize_t))

# pointers to string data must not be valid references
# for the entire lifetime of the ccall
immutable StrArrayStruct
   strings::Ptr{Ptr{Uint8}}
   count::Csize_t
end
StrArrayStruct() = StrArrayStruct(zero(Ptr{Uint8}), zero(Csize_t))

StrArrayStruct{T<:String}(strs::Vector{T}) = begin
    count = length(strs)
    strings = convert(Ptr{Ptr{Uint8}}, Base.c_malloc(sizeof(Ptr{Uint8}) * count))
    for i=1:count
        len = length(strs[i])
        in_ptr  = convert(Ptr{Uint8}, bytestring(strs[i]))
        out_ptr = convert(Ptr{Uint8}, Base.c_malloc(sizeof(Uint8) * (len + 1)))
        unsafe_copy!(out_ptr, in_ptr, len)
        unsafe_store!(out_ptr, zero(Uint8), len + 1) # NULL byte
        unsafe_store!(strings, out_ptr, i)
    end
    return StrArrayStruct(strings, count)
end
StrArrayStruct{T<:String}(str::T) = StrArrayStruct([str])

Base.convert(::Type{Vector{UTF8String}}, sa::StrArrayStruct) = begin
    arr = Array(UTF8String, sa.count)
    for i=1:sa.count
        arr[i] = utf8(bytestring(unsafe_load(sa.strings, i)))
    end
    return arr
end

free!(sa::StrArrayStruct) = begin
    sa_ptr = [sa]
    ccall((:git_strarray_free, libgit2), Void, (Ptr{StrArrayStruct},), sa_ptr)
    return sa_ptr[1]
end

# git config entry struct
immutable ConfigEntryStruct
    name::Ptr{Uint8}
    value::Ptr{Uint8}
    level::Cint
end

typealias DiffDeltaT Cint

@eval begin
    $(Expr(:type, false, :HeaderStruct,
        Expr(:block,
            [Expr(:(::), symbol("c$i"), :Uint8) for i=1:128]...)))
end

immutable DiffHunkStruct
    old_start::Cint
    old_lines::Cint
    new_start::Cint
    new_lines::Cint
    header_len::Csize_t
    header::HeaderStruct
end

immutable DiffFileStruct
    id::Oid
    path::Ptr{Uint8}
    size::Coff_t
    flags::Uint32
    mode::Uint16
end

immutable DiffDeltaStruct
    status::DiffDeltaT
    flags::Uint32
    similarity::Uint16
    nfiles::Uint16
    old_file::DiffFileStruct
    new_file::DiffFileStruct
end

immutable DiffLineStruct
    origin::Cchar
    old_lineno::Cint
    new_lineno::Cint
    num_lines::Cint
    content_len::Csize_t
    content_offset::Coff_t
    # pointer to diff text, ! this is not NULL byte terminated
    content::Ptr{Uint8}
end

# git diff option struct
immutable DiffOptionsStruct
    version::Cuint
    flags::Uint32

    # options controlling which files are in the diff
    ignore_submodules::Cint
    pathspec::StrArrayStruct
    notify_cb::Ptr{Void}
    notify_payload::Ptr{Void}

    # options controlling how the diff text is generated
    context_lines::Uint16
    interhunk_lines::Uint16
    id_abbrev::Uint16
    max_size::Coff_t
    old_prefix::Ptr{Uint8}
    new_prefix::Ptr{Uint8}
end
DiffOptionsStruct() = DiffOptionsStruct(GitConst.DIFF_OPTIONS_VERSION,
                                        zero(Uint32),
                                        GitConst.SUBMODULE_IGNORE_DEFAULT,
                                        StrArrayStruct(),
                                        zero(Ptr{Void}),
                                        zero(Ptr{Void}),
                                        convert(Uint16, 3),
                                        zero(Uint16),
                                        zero(Uint16),
                                        zero(Coff_t),
                                        zero(Ptr{Uint8}),
                                        zero(Ptr{Uint8}))
immutable RemoteHeadStruct
    islocal::Cint
    id::Oid
    lid::Oid
    name::Ptr{Uint8}
end

immutable MergeTreeOptsStruct
    version::Cuint
    flags::Cint
    rename_threshold::Cuint
    target_limit::Cuint
    metric::Ptr{Void}
    automerge_flags::Cint
end

MergeTreeOptsStruct() = MergeTreeOptsStruct(one(Cuint),
                                            zero(Cint),
                                            zero(Cuint),
                                            zero(Cuint),
                                            zero(Ptr{Void}),
                                            zero(Cint))

immutable DiffFileStruct
    id::Oid
    path::Ptr{Uint8}
    size::Coff_t
    flags::Uint32
    mode::Uint16
end

immutable CheckoutOptionsStruct
    version::Cuint
    checkout_strategy::Cuint
    disable_filters::Cint
    dir_mode::Cuint
    file_mode::Cuint
    file_open_flags::Cint
    notify_flags::Cuint
    notify_cb::Ptr{Void}
    notify_payload::Ptr{Void}
    progress_cb::Ptr{Void}
    progress_payload::Ptr{Void}
    paths::StrArrayStruct
    baseline::Ptr{Void}
    target_directory::Ptr{Uint8}
    ancestor_label::Ptr{Uint8}
    our_label::Ptr{Uint8}
    their_label::Ptr{Uint8}
end
CheckoutOptionsStruct() = CheckoutOptionsStruct(one(Cuint),
                                                GitConst.CHECKOUT_SAFE_CREATE,
                                                zero(Cint),
                                                zero(Cuint),
                                                zero(Cuint),
                                                zero(Cint),
                                                zero(Cuint),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                StrArrayStruct(),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Uint8}),
                                                zero(Ptr{Uint8}),
                                                zero(Ptr{Uint8}),
                                                zero(Ptr{Uint8}))
immutable RemoteCallbacksStruct
    version::Cuint
    sideband_progress::Ptr{Void}
    completion::Ptr{Void}
    credentials::Ptr{Void}
    transfer_progress::Ptr{Void}
    update_tips::Ptr{Void}
    payload::Ptr{Void}
end
RemoteCallbacksStruct() = RemoteCallbacksStruct(one(Cuint),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}),
                                                zero(Ptr{Void}))
immutable CloneOptionsStruct
    version::Cuint
    checkout_opts::CheckoutOptionsStruct
    remote_callbacks::RemoteCallbacksStruct
    bare::Cint
    ignore_cert_errors::Cint
    localclone::Cint # GitCloneLocalT
    remote_name::Ptr{Uint8}
    checkout_branch::Ptr{Uint8}
    signature::Ptr{SignatureStruct}
end

CloneOptionsStruct() = CloneOptionsStruct(one(Cuint),
                                          CheckoutOptionsStruct(),
                                          RemoteCallbacksStruct(),
                                          zero(Cint),
                                          zero(Cint),
                                          zero(Cint),
                                          zero(Ptr{Uint8}),
                                          zero(Ptr{Uint8}),
                                          zero(Ptr{SignatureStruct}))
# --------------
# Git Signature
# --------------

type Signature
    name::UTF8String
    email::UTF8String
    time::Int32
    time_offset::Int32
end

# --------------
# Git Repository
# --------------
type GitRepo
    ptr::Ptr{Void}

    function GitRepo(ptr::Ptr{Void}, own::Bool=true)
        @assert ptr != C_NULL
        r = new(ptr)
        own && finalizer(r, free!)
        return r
    end
end

free!(r::GitRepo) = begin
    if r.ptr != C_NULL
        try
            close(r)
        finally
            ccall((:git_repository_free, libgit2), Void, (Ptr{Void},), r.ptr)
            r.ptr = C_NULL
        end
    end
end

Base.convert(::Type{Ptr{Void}}, r::GitRepo) = r.ptr

# -------------
# Git Objects
# -------------
abstract GitObject

free!(o::GitObject) = begin
    if o.ptr != C_NULL
        ccall((:git_object_free, libgit2), Void, (Ptr{Void},), o.ptr)
        o.ptr = C_NULL
    end
end

type GitAnyObject <: GitObject end

GitAnyObject(ptr::Ptr{Void}) = begin
    @assert ptr != C_NULL
    typ = ccall((:git_object_type, libgit2), Cint, (Ptr{Void},), ptr)
    T = gitobj_const_type(typ)
    return T(ptr)
end

type GitBlob <: GitObject
    ptr::Ptr{Void}

    function GitBlob(ptr::Ptr{Void})
        @assert ptr != C_NULL
        this = new(ptr)
        finalizer(this, free!)
        return this
    end
end

Base.convert(::Type{Ptr{Void}}, b::GitBlob) = b.ptr

type GitCommit <: GitObject
    ptr::Ptr{Void}

    function GitCommit(ptr::Ptr{Void})
        @assert ptr != C_NULL
        this = new(ptr)
        finalizer(this, free!)
        return this
    end
end

Base.convert(::Type{Ptr{Void}}, c::GitCommit) = c.ptr

type GitTag <: GitObject
    ptr::Ptr{Void}

    function GitTag(ptr::Ptr{Void})
        @assert ptr != C_NULL
        this = new(ptr)
        finalizer(this, free!)
        return this
    end
end

Base.convert(::Type{Ptr{Void}}, o::GitTag) = o.ptr

type GitTree <: GitObject
    ptr::Ptr{Void}

    function GitTree(ptr::Ptr{Void})
        @assert ptr != C_NULL
        this = new(ptr)
        finalizer(this, free!)
        return this
    end
end

Base.convert(::Type{Ptr{Void}}, o::GitTree) = o.ptr

git_otype(::Type{GitBlob})      = GitConst.OBJ_BLOB
git_otype(::Type{GitCommit})    = GitConst.OBJ_COMMIT
git_otype(::Type{GitTag})       = GitConst.OBJ_TAG
git_otype(::Type{GitTree})      = GitConst.OBJ_TREE
git_otype(::Type{GitAnyObject}) = GitConst.OBJ_ANY

git_otype{T<:GitObject}(o::T) = git_otype(T)

# ---------------
# Git Reftypes
# ---------------
immutable Sym end
const RefType = Union(Oid, Sym)

type GitReference{T<:RefType}
    ptr::Ptr{Void}
end

free!(r::GitReference) = begin
    if r.ptr != C_NULL
        ccall((:git_reference_free, libgit2), Void, (Ptr{Void},), r.ptr)
        r.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, r::GitReference) = r.ptr

# TODO: git branch should be a subtype of GitReference
type GitBranch
    ptr::Ptr{Void}

    function GitBranch(ptr::Ptr{Void})
        @assert ptr != C_NULL
        this = new(ptr)
        finalizer(this, free!)
        return this
    end
end

Base.convert(::Type{Ptr{Void}}, o::GitBranch) = o.ptr

free!(b::GitBranch) = begin
    if b.ptr != C_NULL
        ccall((:git_reference_free, libgit2), Void, (Ptr{Void},), b.ptr)
        b.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, b::GitBranch) = b.ptr

# ---------------
# Git Remote
# ---------------
type GitRemote
    ptr::Ptr{Void}

    function GitRemote(ptr::Ptr{Void})
        @assert ptr != C_NULL
        r = new(ptr)
        finalizer(r, free!)
        return r
    end
end

free!(r::GitRemote) = begin
    if r.ptr != C_NULL
        ccall((:git_remote_free, libgit2), Void, (Ptr{Void},), r.ptr)
        r.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, r::GitRemote) = r.ptr

ext: jl


text:
# 计算哈密顿矩阵
function fill_H(nhil, states, tol_n, J, Jp, Jz)
    # 哈密顿量, 稀疏矩阵
    # H = spzeros(Float64, nhil, nhil)
    # 哈密顿量，full ED
    H = zeros(Float64, nhil, nhil)
    for i in 1:nhil
        ######################################################################
        # 计算1到tol_n - 2个格点的 S_-(i)S_+(j) + S_-(j)S_+(i) 值
        # 包括最近邻和次近邻
        for x in 1:tol_n - 2
            # 最近邻
            index = fill_up_down(states, i, x, x + 1)
            if index != 0.0
                H[i, index] = H[i, index] + 0.5*J
            end
            # 次近邻
            index = fill_up_down(states, i, x, x + 2)
            if index != 0.0
                H[i, index] = H[i, index] + 0.5*Jp
            end
        end
        # 补充计算tol_n - 1到tol_n格点的 S_-(i)S_+(j) + S_-(j)S_+(i) 值
        # 只有最近邻
        index = fill_up_down(states, i, tol_n - 1, tol_n)
        if index != 0.0
            H[i, index] = H[i, index] + 0.5*J
        end
        ######################################################################
        # 计算含Sz项的值，即对角线上的值（不含随机项
        H[i, i] = H[i, i] + fill_SzSz(states[i], tol_n, Jz)
        ######################################################################
    end
    H
end

# 给出态states[label] 的值: S_-(i)S_+(j) + S_-(j)S_+(i)
# 为0则给出0
# 不为0则给出作用后态的label
function fill_up_down(states, label, i, j)
    the_state = states[label]
    i = i - 1 #从第一号平移i - 1位
    j = j - 1 #从第一号平移j - 1位
    i_state = btest(the_state, i)
    j_state = btest(the_state, j)
    # 态相同直接归零
    if i_state == j_state
        return 0.0
    end
    # 态不同则寻找作用后的态的label
    ## 先给出作用后的态
    if (i_state == 0 && j_state == 1)   # 作用后为i+ j-
        fin_state = the_state + 2^i - 2^j
    elseif (i_state == 1 && j_state == 0)   # 作用后为i- j+
        fin_state = the_state - 2^i + 2^j
    end
    searchsortedfirst(states, fin_state)
end

# 计算态 the_state 所有含Sz的项的值（不含随机项
# 返回计算结果
function fill_SzSz(the_state, tol_n, Jz)
    fin_result = 0.0
    for i in 1:tol_n - 1
        i_state = btest(the_state, i - 1)
        j_state = btest(the_state, i)
        #######################################
        # 计算Jz*Sz(i)*Sz(i + 1)
        if i_state == j_state
            fin_result = fin_result + 0.25*Jz
        else
            fin_result = fin_result - 0.25*Jz
        end
        #######################################
    end
    fin_result
end

# R 模型
# 给出含随机项的值
# 返回计算结果
function fill_Sz(states, tol_n, nhil::Int, W, omega)
    fin_result = zeros(Float64, nhil)
    rand_phi = rand(tol_n)
    for j in 1:nhil
        the_state::Int = states[j]
        for i in 1:tol_n
            i_state::Int = btest(the_state, i - 1)
        # 计算Sz(i)*w*cos(2pi*k*i + phi)
            if i_state == 1
                @inbounds fin_result[j] = fin_result[j] + 0.5*W*cos(omega*i + 2.0*pi*(rand_phi[i] - 0.5))
            else
                @inbounds fin_result[j] = fin_result[j] - 0.5*W*cos(omega*i + 2.0*pi*(rand_phi[i] - 0.5))
            end
        end
    end
    fin_result
end

# 给出含随机项的值,待乘正弦波因子
# 返回计算结果
# function fill_Sz(states, tol_n, nhil::Int)
#     fin_result = zeros(Float64, nhil)
#     for j in 1:nhil
#         the_state = states[j]
#         for i in 1:tol_n
#             i_state = btest(the_state, i - 1)
#         # 计算Sz(i)*w*cos(2pi*k*i + phi)
#             if i_state == 1
#                 fin_result[j] = fin_result[j] + 0.5
#             else
#                 fin_result[j] = fin_result[j] - 0.5
#             end
#         end
#     end
#     fin_result
# end

# QP 模型
# 给出含随机项的值
# 返回计算结果
function fill_Sz(states, tol_n, nhil::Int, W, omega, phi)
    fin_result = zeros(Float64, nhil)
    for j in 1:nhil
        the_state::Int = states[j]
        for i in 1:tol_n
            i_state::Int = btest(the_state, i - 1)
        # 计算Sz(i)*w*cos(2pi*k*i + phi)
            if i_state == 1
                @inbounds fin_result[j] = fin_result[j] + 0.5*W*cos(omega*i + phi)
            else
                @inbounds fin_result[j] = fin_result[j] - 0.5*W*cos(omega*i + phi)
            end
        end
    end
    fin_result
end

# 给出Sz为上数目为up的所有态，共有nhil个
function all_required_states!(states, nhil::Int, tol_n::Int, up::Int)
    num_alls::Int = 2^tol_n
    num_count = 0
    for i::Int in 1:num_alls
        num_up = 0
        for j in 0:tol_n-1
            if btest(i, j) == 1
                num_up = num_up + 1
            end
        end
        if num_up == up
            num_count = num_count + 1
            states[num_count] = i
            # println(bitstring(i))
        end
    end

    # 错误处理
    if num_count != nhil
        throw("dimension of Hilbert space mismatch nhil:$nhil num_count$num_count")
    end
end

# 返回s态第j位结果，上为1，下为0
function btest(s::Int, j::Int)
    u0::UInt = 1
    (s >>> j) & u0
end


text:
v = Array{Float64}(undef, 4)
m = Array{Float64}(undef, 4, 3)
@test @inferred(ArrayInterface.axes(v')) === (StaticInt(1):StaticInt(1),Base.OneTo(4))
@test @inferred(ArrayInterface.axes(m')) === (Base.OneTo(3),Base.OneTo(4))
@test ArrayInterface.axes(v', StaticInt(1)) === StaticInt(1):StaticInt(1)
@test ArrayInterface.axes(v, StaticInt(2)) === StaticInt(1):StaticInt(1)

@testset "LazyAxis" begin
    A = zeros(3,4,5);
    SA = @SArray(zeros(3,4,5))
    lz1 = ArrayInterface.LazyAxis{1}(A)
    slz1 = ArrayInterface.LazyAxis{1}(SA)
    lzc = ArrayInterface.LazyAxis{:}(A)
    slzc = ArrayInterface.LazyAxis{:}(SA)

    @test @inferred(first(lz1)) === @inferred(first(slz1))
    @test @inferred(first(lzc)) === @inferred(first(slzc))
    @test @inferred(last(lz1)) === @inferred(last(slz1))
    @test @inferred(last(lzc)) === @inferred(last(slzc))
    @test @inferred(length(lz1)) === @inferred(length(slz1))
    @test @inferred(length(lzc)) === @inferred(length(slzc))
    @test @inferred(Base.to_shape(lzc)) == length(slzc)
    @test @inferred(Base.checkindex(Bool, lzc, 1)) & @inferred(Base.checkindex(Bool, slzc, 1))
    @test axes(lzc)[1] == Base.axes1(lzc) == axes(Base.Slice(lzc))[1] == Base.axes1(Base.Slice(lzc))

    @test @inferred(getindex(lz1, 2)) == 2
    @test @inferred(getindex(lz1, 1:2)) == 1:2
    @test @inferred(getindex(lz1, 1:1:3)) == 1:1:3

    @test @inferred(ArrayInterface.parent_type(ArrayInterface.LazyAxis{:}(A))) <: Base.OneTo{Int}
    @test @inferred(ArrayInterface.parent_type(ArrayInterface.LazyAxis{4}(SA))) <: ArrayInterface.SOneTo{1}
    @test @inferred(ArrayInterface.parent_type(ArrayInterface.LazyAxis{:}(SA))) <: ArrayInterface.SOneTo{60}
    @test ArrayInterface.can_change_size(ArrayInterface.LazyAxis{1,Vector{Any}})

    Aperm = PermutedDimsArray(A, (3,1,2))
    Aview = @view(Aperm[:,1:2,1])
    @test map(parent, @inferred(ArrayInterface.lazy_axes(Aperm))) === @inferred(ArrayInterface.axes(Aperm))
    @test map(parent, @inferred(ArrayInterface.lazy_axes(Aview))) === @inferred(ArrayInterface.axes(Aview))
    @test map(parent, @inferred(ArrayInterface.lazy_axes(Aview'))) === @inferred(ArrayInterface.axes(Aview'))
    @test map(parent, @inferred(ArrayInterface.lazy_axes((1:2)'))) === @inferred(ArrayInterface.axes((1:2)'))

    @test_throws DimensionMismatch ArrayInterface.LazyAxis{0}(A)
end

@testset "`axes(A, dim)`` with `dim > ndims(A)` (#224)" begin
    m = 2
    n = 3
    B = Array{Float64, 2}(undef, m, n)
    b = view(B, :, 1)

    @test @inferred(ArrayInterface.axes(B, 1)) == 1:m
    @test @inferred(ArrayInterface.axes(B, 2)) == 1:n
    @test @inferred(ArrayInterface.axes(B, 3)) == 1:1

    @test @inferred(ArrayInterface.axes(B, static(1))) == 1:m
    @test @inferred(ArrayInterface.axes(B, static(2))) == 1:n
    @test @inferred(ArrayInterface.axes(B, static(3))) == 1:1

    @test @inferred(ArrayInterface.axes(b, 1)) == 1:m
    @test @inferred(ArrayInterface.axes(b, 2)) == 1:1
    @test @inferred(ArrayInterface.axes(b, 3)) == 1:1

    @test @inferred(ArrayInterface.axes(b, static(1))) == 1:m
    @test @inferred(ArrayInterface.axes(b, static(2))) == 1:1
    @test @inferred(ArrayInterface.axes(b, static(3))) == 1:1
end

@testset "SubArray Adjoint Axis" begin
  N = 4; d = rand(N);

  @test @inferred(ArrayInterface.axes_types(typeof(view(d',:,1:2)))) === Tuple{ArrayInterface.OptionallyStaticUnitRange{StaticInt{1}, StaticInt{1}}, Base.OneTo{Int64}}

end
if isdefined(Base, :ReshapedReinterpretArray)
  @testset "ReshapedReinterpretArray" begin
    a = rand(3, 5)
    ua = reinterpret(reshape, UInt64, a)
    @test ArrayInterface.axes(ua) === ArrayInterface.axes(a)
    @test ArrayInterface.axes(ua, 1) === ArrayInterface.axes(a, 1)
    @test @inferred(ArrayInterface.axes(ua)) isa ArrayInterface.axes_types(ua)
    u8a = reinterpret(reshape, UInt8, a)
    @test @inferred(ArrayInterface.axes(u8a)) isa ArrayInterface.axes_types(u8a)
    @test @inferred(ArrayInterface.axes(u8a, static(1))) isa ArrayInterface.axes_types(u8a, 1)
    @test @inferred(ArrayInterface.axes(u8a, static(2))) isa ArrayInterface.axes_types(u8a, 2)
    fa = reinterpret(reshape, Float64, copy(u8a))
    @inferred(ArrayInterface.axes(fa)) isa ArrayInterface.axes_types(fa)
  end
end


text:
include("../src/eis_series.jl")

function test_eis_series_length()
	print("eis_series.length...")

  @test length(eisenstein_series_poly(4, 5)) == 5
  @test length(eisenstein_series_poly(6, 500)) == 500
  @test length(eisenstein_series_poly(12, 1)) == 1
	println("PASS")
end

function test_eis_series_normalisation()
	print("eis_series.normalisation...")

  e4 = eisenstein_series_qexp(4, 5, QQ, "q", "linear")
  @test coeff(e4, 1) == 1
  @test coeff(e4, 0) == 1//240
  e6 = eisenstein_series_qexp(6, 5, QQ, "q", "constant")
  @test coeff(e6, 0) == 1
  @test coeff(e6, 1) == -504
  e12 = eisenstein_series_qexp(12, 5, ZZ, "q", "integral")
  @test coeff(e12, 0) == 691
  @test coeff(e12, 1) == 65520
	println("PASS")
end

function test_eis_series_coefficients()
	print("eis_series.coefficients...")

  e4 = eisenstein_series_qexp(4, 11, QQ, "q", "constant")
  S = parent(e4)
  q = gen(S)
  @test e4 - (1+240*q+2160*q^2+6720*q^3+17520*q^4+30240*q^5+60480*q^6+82560*q^7+140400*q^8+181680*q^9+272160*q^10+O(q^11)) == 0
  e6 = eisenstein_series_qexp(6, 11, ZZ, "q", "integral")
  S = parent(e6)
  q = gen(S)
  @test e6 - (1-504*q-16632*q^2-122976*q^3-532728*q^4-1575504*q^5-4058208*q^6-8471232*q^7-17047800*q^8-29883672*q^9-51991632*q^10+O(q^11)) == 0
  e12 = eisenstein_series_qexp(12, 11, QQ, "q", "linear")
  S = parent(e12)
  q = gen(S)
  @test e12 - (691//65520+q+2049*q^2+177148*q^3+4196353*q^4+48828126*q^5+362976252*q^6+1977326744*q^7+8594130945*q^8+31381236757*q^9+100048830174*q^10+O(q^11)) == 0
	println("PASS")
end

function test_eis_series()
	test_eis_series_length()
  test_eis_series_normalisation()
  test_eis_series_coefficients()

	println("")
end


text:
@testset "ensembles.jl" begin

    print("\n\n ####### ensembles.jl #######\n\n")

    location = SaveLocation(;prefix=PREFIX)

    pdd = PositionDataDescriptor(:box, 1, 7, 20, [0.1, 0.2], location)
    model = RandomPositionsXXZWithXField(pdd, PowerLaw(6))

    edd = EDDataDescriptor(model, location; suffix="full_ed_threaded")
    ensdd = EnsembleDataDescriptor(edd)

    ensdata = load_or_create(ensdd, location) # should load the ED data since previous test saved it, and save to file without suffix
    save(ensdata)

    ensdata2 = load_ensemble(model, location) # load previously generated data
    ## No idea for good test
    # check access of ensembles ?
    @test ensdata.canonical == ensdata2[:,:,:,:,ENSEMBLE_INDICES[:canonical]]
    @test ensdata.microcanonical == ensdata2[:,:,:,:,ENSEMBLE_INDICES[:microcanonical]]
    @test ensdata.diagonal == ensdata2[:,:,:,:,ENSEMBLE_INDICES[:diagonal]]
    # check attribute forwarding
    @test ensdd.fields == edd.fields
    @test ensdata2.descriptor == ensdata.descriptor
end


text:
## @lc code=start
using LeetCode

function find_ball(grid::Vector{Vector{Int}})::Vector{Int}
    function next_pos(line::Vector{Int}, i::Int)::Int
        i == 0 && return 0
        i > 1 && line[i] == line[i-1] == -1 && return i - 1
        i < n && line[i] == line[i+1] == 1 && return i + 1
        return 0
    end
    n = length(grid[1])
    f(ans, line) = [next_pos(line, j) for j in ans]
    foldl(f, [1:n, grid...]).-1
end

## @lc code=end
## @lc test=start
@testset "1706.where-will-the-ball-fall.jl" begin
    @test find_ball([[1,1,1,-1,-1],[1,1,1,-1,-1],[-1,-1,-1,1,1],[1,1,1,1,-1],[-1,-1,-1,-1,-1]]) == [1,-1,-1,-1,-1]
    @test find_ball([[1,1,1,1,1,1],[-1,-1,-1,-1,-1,-1],[1,1,1,1,1,1],[-1,-1,-1,-1,-1,-1]]) == [0,1,2,3,4,-1]
    @test find_ball([[-1]]) == [-1]
end
## @lc test=end

## @lc info=start
# link: [solution to 1706]()
## @lc info=end


text:
# using DataFrames
using Statistics
using DelimitedFiles
using LinearAlgebra
abstract type EstimationMethod end
struct NonParametric <: EstimationMethod end
struct Parametric <: EstimationMethod end

"""
VAR model can be expressed in state-space form:
```math
y_t = Qz_t
```
```
z_t = Mz_{t-1}+Gu_t
```
where ``z_t`` has ``y_t`` and lagged ``y_t``s
"""
struct VARModel{TA <: AbstractArray, TAm <: AbstractArray, TMm <: AbstractMatrix}
    y::TA
    nlag::Int
    withconst::Bool
    initperiod::Int
    lastperiod::Int
    T::Int
    ns::Int
    resid::TAm
    betahat::TAm
    M::TMm
    Q::TMm
    G::TMm
    seps::TMm
end
mutable struct FactorEstimateStats
    T::Int # number of data periods used for factor estimation
    ns::Int # number of data series
    nobs # total number of observations (=sum_i T_i)
    tss
    ssr
    R2::Vector{Union{Missing, Float64}}
end
struct DFMModel
    data::AbstractArray
    inclcode::Vector{Int}
    T::Int  # number of whole data periods
    ns::Int
    nt_min_factor_estimation::Int
    nt_min_factorloading_estimation::Int
    initperiod::Int
    lastperiod::Int
    nfac_o::Int
    nfac_u::Int
    nfac_t::Int
    tol::Float64
    fes::FactorEstimateStats
    factor::AbstractArray
    lambda::AbstractArray
    uar_coef::AbstractArray
    uar_ser::Vector{Float64}
    n_uarlag::Int
    n_factorlag::Int
    factor_var_model::VARModel
end
"""
Note: `factor` and `factor_var_model.y` are actually using same memory.
"""
function DFMModel(data, inclcode,
    nt_min_factor_estimation::Integer, nt_min_factorloading_estimation::Integer,
    initperiod::Integer, lastperiod::Integer,
    nfac_o::Integer, nfac_u::Integer, tol, n_uarlag::Integer, n_factorlag::Integer)
    size(data, 2) == length(inclcode) || error("length of inclcode must equal to number of data series")
    initperiod < lastperiod || error("initperiod must be smaller than lastperiod")
    ((n_uarlag > 0) && (n_factorlag > 0)) || error("n_uarlag and n_factorlag must be positive")

    T, ns = size(data)
    nfac_t = nfac_o+nfac_u
    fes = FactorEstimateStats(lastperiod - initperiod + 1,
                              count(inclcode.==1),
                              missing, missing, missing,
                              Vector{Union{Missing, Float64}}(undef, count(inclcode.==1)))
    factor = Matrix{Union{Missing, Float64}}(missing, T, nfac_t)
    lambda = Matrix{Float64}(undef, ns, nfac_t)
    uar_coef = Matrix{Float64}(undef, ns, n_uarlag)
    uar_ser = Vector{Float64}(undef, ns)
    factor_var_model = VARModel(factor, n_factorlag, initperiod= initperiod,
                                lastperiod=lastperiod)
    return DFMModel(data, vec(inclcode), T, ns,
                    nt_min_factor_estimation, nt_min_factorloading_estimation,
                    initperiod, lastperiod, nfac_o, nfac_u, nfac_t,
                    tol, fes, factor, lambda, uar_coef, uar_ser,
                    n_uarlag, n_factorlag, factor_var_model)
end

function drop_missing_row(A::AbstractMatrix)
    tmp = .!any(ismissing.(A), dims=2)
    return A[vec(tmp), :], tmp
end

function drop_missing_col(A::AbstractMatrix)
    tmp = drop_missing_row(A')
    return tmp[1]', tmp[2]'
end

function pca_score(X, nfac_u::Integer)
    _, _, V = svd(X)
    score = (X*V)[:, 1:nfac_u]
    return score
end

"""
##### Arguments
- `y::AbstractVector`: length `T` Vector or `TxN` Matrix where `T`
                       is sample size and `N` is the number of
                       regressed variables
- `X::AbstractArray`: `TxK` Matrix where `K` is the number of
                      regressors
##### Outputs
- `b`: OLS estimator of the coefficients
- `e`: residual
"""
function ols(y::AbstractArray, X::AbstractArray)
    b = X\y
    e = y - X*b
    return b, e
end

abstract type OLSSkipRule end
struct Balanced <: OLSSkipRule end
struct Unbalanced <: OLSSkipRule end
"""
do OLS ignoring samples containing missing values
##### Outputs
- `b`: OLS estimator of the coefficients
- `e`: residual
- `numrow`: BitArray inidicating rows used to estimate
"""
function ols_skipmissing(y::AbstractMatrix, X::AbstractArray, ::Balanced)
    N = size(y, 2)
    tmp, numrow = drop_missing_row([y X])
    y_used, x_used = tmp[:, 1:N], tmp[:, N+1:end]
    b, e = ols(y_used, x_used)
    return b, e, vec(numrow)
end
function ols_skipmissing(y::AbstractVector, X::AbstractArray, method::Balanced)
    b, e, numrow = ols_skipmissing(reshape(y, size(y, 1), size(y, 2)), X, method)
    return b, vec(e), numrow
end
"""
##### Arguments
- `y::AbstractMatrix`: `TxN`
- `X::AbstractArray`: `TxK` Matrix or `T` Vector
"""
function ols_skipmissing(y::AbstractMatrix, X::AbstractArray, ::Unbalanced)
    if size(y, 1) != size(X, 1)
        error("Sample size must be same")
    end
    T, N = size(y)
    b = Matrix{Float64}(undef, size(X, 2), N)
    e = Matrix{Union{Missing, Float64}}(missing, T, N)
    numrow = BitArray(undef, T, N)
    for i=1:N
        tmp = ols_skipmissing(y[:, i], X, Balanced())
        b[:, i] = tmp[1]
        e[tmp[3], i] = tmp[2]
        numrow[:, i] = tmp[3]
    end
    return b, e, numrow
end

function lagmat(X::AbstractArray, lags::AbstractVector)
    nc = size(X, 2)
    Xlag = Matrix{Union{Missing, Float64}}(missing, size(X, 1), nc*length(lags))
    for (i, lag) in enumerate(lags)
        Xlag[lag+1:end, nc*(i-1)+1:nc*i] .= X[1:end-lag, :]
    end
    return Xlag
end
lagmat(X::AbstractArray, lag::Integer) = lagmat(X, [lag])

function uar(y::AbstractVector, n_lags::Integer)
    x = lagmat(y, 1:n_lags)
    arcoef, ehat, _ = ols_skipmissing(y, x, Balanced())
    ssr = dot(ehat, ehat)
    ser = sqrt(ssr/(size(x, 1)-size(x, 2)))
    return arcoef, ser
end
"""
estimate factor by iteration using balanced data
"""
function estimate_factor!(m::DFMModel, max_iter::Integer=100000000,
                         computeR2::Bool=true)
    data = m.data
    initperiod, lastperiod, nt_min, nfac_u, nfac_o, tol =
        m.initperiod, m.lastperiod, m.nt_min_factor_estimation,
        m.nfac_u, m.nfac_o, m.tol
    # use part of the data
    est_data = data[:, m.inclcode.==1]
    xdata = est_data[initperiod:lastperiod, :]

    # preprocess data to have unit standard error
    xdata_standardized, xdatastd = standardize_data(xdata)

    m.fes.tss = sum(skipmissing(xdata_standardized.^2))
    m.fes.nobs = length(xdata_standardized[.!ismissing.(xdata_standardized)])

    xbal, _ = drop_missing_col(xdata_standardized)

    # Get initial F_t given Lambda_t using PCA
    f = pca_score(xbal, nfac_u)
    m.fes.ssr = 0
    diff = 1000
    lambda = Array{Union{Missing, Float64}}(undef, m.fes.ns, m.nfac_t)
    for iter = 1:max_iter
        ssr_old = m.fes.ssr
        # given F_t, get Lambda_t
        for i = 1:m.fes.ns
            tmp, usedrows = drop_missing_row([xdata_standardized[:, i] f])
            if count(usedrows) >= nt_min
                lambda[i, :] =
                    ols_skipmissing(xdata_standardized[:, i], f, Balanced())[1]'
            end
        end
        # given Lambda_t, get F_t by regressing X_t on Lambda_t for each t
        tmp = ols_skipmissing(xdata_standardized', lambda[:, nfac_o+1:end], Unbalanced())
        f, ehat = tmp[1]', tmp[2]
        m.fes.ssr = sum(sum(skipmissing(ehat.^2)))
        diff = abs(ssr_old - m.fes.ssr)
        diff >= tol*m.fes.T*m.fes.ns || break
        # println("diff = ", diff)
    end
    m.factor[initperiod:lastperiod,  :] = f
    if computeR2
        for i=1:m.fes.ns
            tmp = drop_missing_row([xdata_standardized[:, i] f])[1]
            if size(tmp, 1) >= nt_min
                _, ehat = ols(tmp[:, 1], tmp[:, 2:end])
                m.fes.R2[i] = compute_r2(tmp[:, 1], ehat)[1]
            end
        end
    end
    return nothing
end

function estimate_factor_loading!(m::DFMModel)
    data, initperiod, lastperiod, fac, nt_min, nfac_t, n_uarlag =
        m.data, m.initperiod, m.lastperiod, m.factor,
        m.nt_min_factorloading_estimation, m.nfac_t, m.n_uarlag
    n_series = size(data, 2)
    r2_mat = Vector{Float64}(undef, n_series)
    for is = 1:n_series
        tmp, numrow = drop_missing_row([data[initperiod:lastperiod, is] fac[initperiod:lastperiod, :]])
        if count(numrow) >= nt_min # if available sample size is large enough
            b, uhat = ols(tmp[:, 1], [tmp[:, 2:end] ones(count(numrow))])
            y_used = data[initperiod:lastperiod, is][vec(numrow), :]
            m.lambda[is, :] .= b[1:end-1]
            r2_mat[is], _, _ = compute_r2(y_used, uhat)
            if r2_mat[is] < 0.9999
                arcoef, ser = uar(uhat, n_uarlag)
            else
                arcoef, ser = zeros(n_uarlag, 1), 0.0
            end
        end
        m.uar_coef[is, :] = arcoef'
        m.uar_ser[is] = ser
    end
    return nothing
end
function VARModel(y::AbstractArray, nlag::Integer=1;
                  withconst::Bool=true,
                  initperiod::Integer=1, lastperiod::Integer=size(y, 1))
    T, ns = size(y, 1), size(y, 2)
    resid = Array{Union{Missing, Float64}}(missing, size(y))
    betahat = Matrix{Union{Missing, Float64}}(missing, ns*nlag+withconst, ns)
    M = Matrix{Union{Missing, Float64}}(missing, ns*nlag, ns*nlag)
    Q = Matrix{Union{Missing, Float64}}(missing, ns, ns*nlag)
    G = Matrix{Union{Missing, Float64}}(missing, ns*nlag, ns)
    seps = Matrix{Union{Missing, Float64}}(missing, ns, ns)
    return VARModel(y, nlag, withconst, initperiod, lastperiod, T, ns, resid, betahat, M, Q, G, seps)
end
function estimate_var!(varm::VARModel, compute_matrices::Bool=true)
    initperiod, lastperiod = varm.initperiod, varm.lastperiod
    withconst, nlag = varm.withconst, varm.nlag
    resid, seps = varm.resid, varm.seps

    y_restricted = varm.y[initperiod:lastperiod, :]

    # regressors
    withconst || (x = lagmat(y_restricted, 1:nlag))
    !withconst || (x = [ones(lastperiod-initperiod+1) lagmat(y_restricted, 1:nlag)])

    # do OLS ignoring the samples containing NaN
    betahat, ehat, numrows = ols_skipmissing(y_restricted, x, Balanced())
    varm.betahat .= betahat

    T_used = count(numrows) # used sample size
    K = size(x, 2) # number of regressors

    ndf = T_used - K # degree of freedom (T-K)
    seps .= ehat'*ehat/ndf # covariance matrix of error term
    resid[initperiod- 1 .+ findall(numrows), :] .= ehat

    !compute_matrices || fill_matrices!(varm, betahat)
    return nothing
end
function fill_matrices!(varm::VARModel, betahat::Array)
    ns, nlag = varm.ns, varm.nlag
    M, Q, G = varm.M, varm.Q, varm.G

    b = betahat[2:end, :]' # now, each row corresponds to each equation

    M .= zeros(ns*nlag, ns*nlag)
    M[1:ns, :] .= b # coefficients of VAR
    M[ns+1:end, 1:end-ns] .= Matrix{Float64}(I, ns*nlag-ns, ns*nlag-ns)　# lag part

    Q .= zeros(ns, ns*nlag)
    Q[1:ns, 1:ns] .= Matrix{Float64}(I, ns, ns)
    G .= zeros(ns*nlag, ns)
    G[1:ns, 1:ns] .= ((cholesky(varm.seps)).U)'
    return nothing
end

function standardize_data(data::AbstractArray)
    datamean = [mean(collect(skipmissing(data[:, i]))) for i = 1:size(data, 2)]'
    # # make correction (which I don't understand why being needed)
    tmp = size(data, 1) .- sum(ismissing.(data), dims = 1)
    tmp = (tmp.-1)./tmp
    datastd = [std(collect(skipmissing(data[:, i]))) for i = 1:size(data, 2)]'.*(tmp.^.5)
    data_standardized = (data .- datamean)./datastd
    return data_standardized, datastd
end

"""
estimate DFM Model non-parametrically
"""
function estimate!(m::DFMModel, ::NonParametric=NonParametric())

    # estimate factor using balanced data
    estimate_factor!(m)

    # estimate factor loading using full sample
    estimate_factor_loading!(m)

    # estimate the equation of factor evolution
    estimate_var!(m.factor_var_model)

    return nothing
end
compute_series(dfmm::DFMModel, is::Integer) = dfmm.factor*dfmm.lambda[is, :]
detrended_year_growth(y::AbstractVector) = vec(sum(lagmat(y, 0:3), dims=2))

find_row_number(date::Tuple{Int, Int}, dates) =
    findall([date == dataset.calds[i] for i=1:length(dataset.calds)])[1]

function compute_r2(y::AbstractArray, e::AbstractVector)
    ssr = dot(e, e)
    tss = dot(y.-mean(y), y.-mean(y))
    return 1-(ssr/tss), ssr, tss
end
function compute_bw_weight(bw_para::Integer)
    bw_weight = zeros(2bw_para+1)
    trend = Vector{Int64}(undef, 2bw_para+1)
    for i=0:100
        trend[bw_para+1+i] = i
        trend[bw_para+1-i] = -i
        bw_weight[bw_para+1+i] = 15/16*(1-(trend[bw_para+1+i]/bw_para)^2)^2
        bw_weight[bw_para+1-i] = 15/16*(1-(trend[bw_para+1-i]/bw_para)^2)^2
    end
    bw_weight = bw_weight./sum(bw_weight)
    return bw_weight, trend
end
function gain(h::AbstractVector, w::Real)
    # Calculate Gain of filter h at Frequency w
    z = exp(-w*im)
    h1 = h[1]
    z1 = 1
    for i = 2:length(h)
        z1 = z1*z
        h1 = h1 + h[i]*z1
    end
    g2=h1*h1'
    fgain = sqrt(g2)
    return fgain
end

function bai_ng_criterion(m::DFMModel)
    fes = m.fes
    nbar = fes.nobs/fes.T # average observation per period
    g = log(min(nbar, fes.T))*(nbar+fes.T)/fes.nobs
    bn_icp = log(fes.ssr/fes.nobs)+ m.nfac_t*g
    return bn_icp
end
"""
- `nfac_max::Integer`: maximum number of factors
"""
struct FactorNumberEstimateStats
    bn_icp
    ssr_static
    R2_static
    aw_icp
    ssr_dynamic
    R2_dynamic
    tss::Float64
    nobs::Int
    T::Int
end
"""
- `m::DFMModel`: `DFMModel` specifying the model except number of unobservable
                 factors.
"""
function estimate_factor_numbers(m::DFMModel, nfacs::Union{Real, AbstractVector})
    max_nfac = maximum(nfacs)
    bn_icp = Vector{Union{Missing, Float64}}(undef, max_nfac)
    ssr_static = Vector{Float64}(undef, max_nfac)
    R2_static = Matrix{Union{Missing, Float64}}(undef, m.fes.ns, max_nfac)
    aw_icp = Matrix{Union{Missing, Float64}}(undef, max_nfac, max_nfac)
    ssr_dynamic = Matrix{Union{Missing, Float64}}(undef, max_nfac, max_nfac)
    R2_dynamic = Array{Union{Missing, Float64}}(undef, m.fes.ns, max_nfac, max_nfac)

    global tss, nobs, T
    for (i, nfac) = enumerate(1:max_nfac)
        dfmm = DFMModel(m.data, m.inclcode,
                m.nt_min_factor_estimation, m.nt_min_factorloading_estimation,
                m.initperiod, m.lastperiod, m.nfac_o, nfac, m.tol, m.n_uarlag, m.n_factorlag)
        estimate_factor!(dfmm)
        bn_icp[i] = bai_ng_criterion(dfmm)
        ssr_static[i] = dfmm.fes.ssr
        R2_static[:, i] = dfmm.fes.R2
        aw_icp[1:nfac, i], ssr_dynamic[1:nfac, i], R2_dynamic[:, 1:nfac, i] =
            amengual_watson_test(dfmm, 4)
        global tss = dfmm.fes.tss
        global nobs = dfmm.fes.nobs
        global T = dfmm.fes.T
    end
    return FactorNumberEstimateStats(bn_icp, ssr_static, R2_static,
                                     aw_icp, ssr_dynamic, R2_dynamic,
                                     tss, nobs, T)
end

function amengual_watson_test(m::DFMModel, nper::Integer)
    factor = m.factor
    T, ns, nfac_static = m.T, m.fes.ns, m.nfac_t
    nvar_lag = m.factor_var_model.nlag
    est_data = m.data[:, m.inclcode.==1]

    # Construct lags of factors and residuals for est_data
    x = [ones(T) lagmat(factor, 1:nvar_lag)]
    est_data_res = Array{Union{Missing, Float64}}(undef, T, ns)
    for is = 1:ns
        tmp, usedrows = drop_missing_row([est_data[:, is] x])
        y = tmp[:, 1]
        z = tmp[:, 2:end]
        ndf = size(z, 1)-size(z, 2)
        if ndf >= m.nt_min_factor_estimation  # Minimum degrees of freedom for series
            b, e = ols(y, z)
            est_data_res[findall(vec(usedrows)), is] = e
        end
    end

    # Carry out calculations for number of dynamic factors
    ssr = Array{Float64}(undef, nfac_static)
    r2  = Array{Union{Missing, Float64}}(undef, ns, nfac_static)
    aw  = Array{Float64}(undef, nfac_static)
    for nfac = 1:nfac_static
        dfmm = DFMModel(est_data_res, ones(count(m.inclcode.==1)),
                m.nt_min_factor_estimation, m.nt_min_factorloading_estimation,
                m.initperiod+4, m.lastperiod, 0, nfac, m.tol, m.n_uarlag, m.n_factorlag)
        estimate_factor!(dfmm)
        aw[nfac] = bai_ng_criterion(dfmm)
        ssr[nfac] = dfmm.fes.ssr
        r2[:, nfac] = dfmm.fes.R2
    end
    return aw, ssr, r2
end

"""
compute impulse response function of `VARModel` with one standard error shock
```math
y_t = Qz_t
```
```
z_t = Mz_{t-1}+Gu_t
```
##### Arguments
- `varm::VARModel`
- `shock_ids`: `<:AbstractVector` or `<:Real` or `:all`.
- `T::Integer`:  horizon of IRF
##### Return
- `irfs::Array` IRFs. It is 2D if `shock_ids<:Real` and 3D otherwise.
                First dimension is each series, second dimension is time,
                third dimension is type of shock (if exist).
"""
function impulse_response(varm::VARModel, shock_ids::AbstractVector, T::Integer)
    irfs = Array{Float64, 3}(undef, size(varm.Q, 1), T, length(shock_ids))
    for (i, shock_id) in enumerate(shock_ids)
        compute_irf_single_shock!(irfs, varm, i, shock_id, T)
    end
    return irfs
end
function compute_irf_single_shock!(irfs::AbstractArray, varm::VARModel,
                       i::Integer, shock_id::Integer, T::Integer)
    x = varm.G[:, shock_id]
    for t = 1:T
        irfs[:, t, i] = varm.Q * x
        x .= varm.M * x
    end
    return nothing
end
function impulse_response(varm::VARModel, shock_id::Real, T::Integer)
    irfs = Matrix{Float64}(undef, size(varm.Q, 1), T)
    compute_irf_single_shock!(irfs, varm, x, 1, shock_id, T)
    return irfs
end
impulse_response(varm::VARModel, shock_id::Symbol, T::Integer) =
    impulse_response(varm, Val(shock_id), T)
impulse_response(varm::VARModel, shock_ids::Val{:all}, T::Integer) =
    impulse_response(varm, 1:size(varm.G, 2), T)


text:
############################################################
# judiJacobian ###############################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiJacobian, judiJacobianException, subsample

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

abstract type judiAbstractJacobian{DDT<:Number,RDT<:Number} <: joAbstractLinearOperator{DDT,RDT} end

# Type for linear operator representing  Pr*A(m)^-1*Ps,
# i.e. it includes source and receiver projections
struct judiJacobian{DDT<:Number,RDT<:Number} <: judiAbstractJacobian{DDT,RDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    model::Model
    source::judiVector
    recGeometry::Geometry
    options::Options
    fop::Function              # forward
    fop_T::Union{Function, Nothing}  # transpose
end

mutable struct judiJacobianException <: Exception
    msg :: String
end

############################################################
## Constructor
"""
    judiJacobian(F,q)

Create a linearized modeling operator from the non-linear modeling operator `F` and \\
the source `q`. `F` is a full modeling operator including source/receiver projections.

Examples
========

1) `F` is a modeling operator without source/receiver projections:

    J = judiJacobian(Pr*F*Ps',q)

2) `F` is the combined operator `Pr*F*Ps'`:

    J = judiJacobian(F,q)

"""
function judiJacobian(F::judiPDEfull, source::judiVector; DDT::DataType=Float32, RDT::DataType=DDT, options=nothing)
# JOLI wrapper for nonlinear forward modeling
    compareGeometry(F.srcGeometry, source.geometry) == true || judiJacobianException("Source geometry mismatch")
    (DDT == Float32 && RDT == Float32) || throw(judiJacobianException("Domain and range types not supported"))
    m = n_samples(F.recGeometry, F.info)
    n = F.info.n

    isnothing(options) && (options = F.options)
    return J = judiJacobian{Float32,Float32}("linearized wave equation", m, n,
        F.info, F.model, source, F.recGeometry, options, bornop, adjbornop)
end

judiJacobian(J::judiJacobian{DDT,RDT}; name=J.name, m=J.m, n=J.n, info=J.info, model=J.model, source=J.source,
     geom=J.recGeometry, opt=J.options, fop=J.fop, fop_T=J.fop_T) where {DDT, RDT} =
            judiJacobian{DDT,RDT}(name, m, n, info, model, source, geom, opt, fop, fop_T)
############################################################
## overloaded Base functions

# conj(judiJacobian), real so return itself
conj(A::judiAbstractJacobian{DDT,RDT}) where {DDT,RDT} = A
# transpose(judiJacobian)
transpose(A::judiAbstractJacobian{DDT,RDT}) where {DDT,RDT} = judiJacobian(A; m=A.n, n=A.m, fop=A.fop_T, fop_T=A.fop)
# adjoint(judiJacobian)
adjoint(A::judiAbstractJacobian{DDT,RDT}) where {DDT,RDT}  = judiJacobian(A; m=A.n, n=A.m, fop=A.fop_T, fop_T=A.fop)

############################################################
## overloaded Base *(...judiJacobian...)

# *(judiJacobian,vec)
function *(A::judiAbstractJacobian{ADDT,ARDT}, v::AbstractVector{Float32}) where {ADDT,ARDT}
    A.n == size(v,1) || throw(judiJacobianException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    V = A.fop(A, v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiJacobian,vec):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

*(A::judiAbstractJacobian{ADDT,ARDT}, v::AbstractMatrix{vDT}) where {ADDT,ARDT,vDT} = *(A, vec(v))
*(A::judiAbstractJacobian{ADDT,ARDT}, v::AbstractVector{Float64}) where {ADDT,ARDT} = *(A, jo_convert(Float32, v, false))

# *(judiJacobian,judiVector)
function *(A::judiAbstractJacobian{ADDT,ARDT}, v::judiVector{vDT, AT}) where {ADDT,ARDT,vDT, AT}
    A.n == size(v,1) || throw(judiJacobianException("Shape mismatch: A:$(size(A)), v: $(size(v))"))
    jo_check_type_match(ADDT,vDT,join(["DDT for *(judiJacobian,judiVector):",A.name,typeof(A),vDT]," / "))
    compareGeometry(A.recGeometry,v.geometry) == true || throw(judiJacobianException("Geometry mismatch"))
    V = A.fop(A, v)
    jo_check_type_match(ARDT,eltype(V),join(["RDT from *(judiJacobian,judiVector):",A.name,typeof(A),eltype(V)]," / "))
    return V
end

# *(num,judiJacobian)
*(a::Number,A::judiAbstractJacobian{ADDT,ARDT}) where {ADDT,ARDT} =  judiJacobian(A; source=a*A.source)
# -(judiJacobian)
-(A::judiAbstractJacobian{DDT,RDT}) where {DDT,RDT} = judiJacobian(A; source=-A.source)

############################################################
## overloaded Bases +(...judiJacobian...), -(...judiJacobian...)

# +(judiJacobian,num)
function +(A::judiAbstractJacobian{ADDT,ARDT}, b::Number) where {ADDT,ARDT}
    return judiJacobian(A; fop=v1=>A.fop(A, v1)+joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
                   fop_T=v2 -> A.fop_T(A, v2)+joConstants(A.n,A.m,b;DDT=ADDT,RDT=ARDT)*v2)
end

# -(judiJacobian,num)
function -(A::judiAbstractJacobian{ADDT,ARDT},b::Number) where {ADDT,ARDT}
    return judiJacobian(A; fop=v1 -> A.fop(A, v1)-joConstants(A.m,A.n,b;DDT=ADDT,RDT=ARDT)*v1,
                   fop_T=v2 -> A.fop_T(A, v2)-joConstants(A.n,A.m,b;DDT=ADDT,RDT=ARDT)*v2)
end

############################################################
## Additional overloaded functions

# Subsample Jacobian
function subsample(J::judiJacobian{ADDT,ARDT}, srcnum) where {ADDT,ARDT}
    recGeometry = subsample(J.recGeometry, srcnum)
    nsrc = typeof(srcnum) <: Int ? 1 : length(srcnum)
    info = Info(J.info.n, nsrc, J.info.nt[srcnum])
    m = n_samples(recGeometry, info)
    return judiJacobian(J; m=m, info=info, source=J.source[srcnum],geom=recGeometry, opt=subsample(J.options, srcnum))
end

getindex(J::judiAbstractJacobian, a) = subsample(J, a)

############################################################
## Forward/adjoint function to avoid unecessary extra declaration

function bornop(J::judiJacobian, v)
    srcnum = 1:J.info.nsrc
    return time_modeling(J.model, J.source.geometry, J.source.data, J.recGeometry, nothing, v, srcnum, 'J', 1, J.options)
end

function adjbornop(J::judiJacobian, w)
    srcnum = 1:J.info.nsrc
    return time_modeling(J.model, J.source.geometry, J.source.data, J.recGeometry,
                         process_input_data(w, J.recGeometry, J.info), nothing, srcnum, 'J', -1, J.options)
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using ManifoldsBase
using Test

struct ErrorTestManifold <: AbstractManifold{ℝ} end

function ManifoldsBase.check_point(::ErrorTestManifold, x)
    if any(u -> u < 0, x)
        return DomainError(x, "<0")
    end
    return nothing
end
function ManifoldsBase.check_vector(M::ErrorTestManifold, x, v)
    mpe = check_point(M, x)
    mpe === nothing || return mpe
    if any(u -> u < 0, v)
        return DomainError(v, "<0")
    end
    return nothing
end

@testset "Domain errors" begin
    M = ErrorTestManifold()
    @test isa(check_point(M, [-1, 1]), DomainError)
    @test check_point(M, [1, 1]) === nothing
    @test !is_point(M, [-1, 1])
    @test is_point(M, [1, 1])
    @test_throws DomainError is_point(M, [-1, 1], true)

    @test isa(check_vector(M, [1, 1], [-1, 1]), DomainError)
    @test check_vector(M, [1, 1], [1, 1]) === nothing
    @test !is_vector(M, [1, 1], [-1, 1])
    @test is_vector(M, [1, 1], [1, 1])
    @test_throws DomainError is_vector(M, [1, 1], [-1, 1], true)
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ 3df9d1c0-164c-11eb-335f-edb4366d2625
begin
	using Statistics
	using Images
	using FFTW
	using Plots
	using DSP
	using ImageFiltering
	using PlutoUI
	using OffsetArrays
end

# ╔═╡ 41d37ef0-164b-11eb-0940-636023edbe0d
md"""
## Image Convolutions
"""

# ╔═╡ f8db5200-1680-11eb-1d10-cbc2b31fc024
function shrink_image(image, ratio=5)
	(height, width) = size(image)
	new_height = height ÷ ratio - 1
	new_width = width ÷ ratio - 1
	list = [
			mean(image[
				ratio * i:ratio * (i + 1),
				ratio * j:ratio * (j + 1),
			])
			for j in 1:new_width
			for i in 1:new_height
	]
	reshape(list, new_height, new_width)
end

# ╔═╡ 35aa24a0-1680-11eb-065b-7973a0dab607
begin
	url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/03/TheOreoCat.jpeg/900px-TheOreoCat.jpeg"
	download(url, "cat_in_a_hat.jpg")
	large_image = load("cat_in_a_hat.jpg")
	image = shrink_image(large_image, 7)
end

# ╔═╡ 2a905470-1687-11eb-30c3-5f635c6c0aeb
kernel = Kernel.gaussian((1, 1))

# ╔═╡ c2515c00-1687-11eb-3cf2-1385aa2ce3d6
function show_colored_kernel(kernel)
	to_rgb(x) = RGB(max(-x, 0), max(x,0), 0)
	to_rgb.(kernel) / maximum(abs.(kernel))
end

# ╔═╡ 0c4eb462-1688-11eb-0f56-cf09d1f53f0a
show_colored_kernel(kernel)

# ╔═╡ 9c21a2b0-16af-11eb-3fae-6d8927359b72
function clamp_at_boundary(M, i, j)
	return M[
		clamp(i, 1, size(M, 1)),
		clamp(j, 1, size(M, 2)),
	]
end

# ╔═╡ 307c6610-16b6-11eb-159c-f18f672027df
begin
	I = [1 2 3; 8 4 9]
	size(I, 2)
end

# ╔═╡ 1592aef0-1688-11eb-2d74-35ce0cb1ddc3
function convolve(M, kernel, M_index_function=clamp_at_boundary)
	height = size(kernel, 1)
	width = size(kernel, 2)

	half_height = height ÷ 2
	half_width = width ÷ 2

	new_image = similar(M)

	# (i, j) loop over the original image
	@inbounds for i in 1:size(M, 1)
		for j in 1:size(M, 2)
			# (k, l) loop over the neighbouring pixels
			new_image[i, j] = sum([
						kernel[k, l] * M_index_function(M, i - k, j - l)
						for k in -half_height:-half_height + height - 1
						for l in -half_width:-half_width + width - 1
					])
		end
	end
	return new_image
end

# ╔═╡ 9650d1b0-16b6-11eb-061b-d94361398d7b
begin
	K = OffsetArray(gaussian((3,3), 0.25), -1:1, -1:1)
	U = rand(1.0:100.0, 6, 6);
	convolve(U, K)
end

# ╔═╡ 765a986a-173c-11eb-29e6-cb155437d4ea
convolve(image, Kernel.gaussian((3, 3)))

# ╔═╡ 7767a37e-173c-11eb-39a8-0d1bd7be9077
convolve(image, Kernel.gaussian((10, 10)))

# ╔═╡ 7785d7e0-173c-11eb-11b4-ebc2a9c3815a
sharpen_kernel = centered([
	-0.5 -1.0 -0.5
	-1.0  7.0 -1.0
	-0.5 -1.0 -0.5
])

# ╔═╡ cea7f490-1746-11eb-2a66-8f75edaec775
edge_detection_kernel_horizontal = Kernel.sobel()[1]

# ╔═╡ b9733322-1747-11eb-1da7-3d4b992d60a2
show_colored_kernel(edge_detection_kernel_horizontal)

# ╔═╡ ed2204f0-1747-11eb-1662-9d8199f940a1
edge_detection_kernel_vertical = Kernel.sobel()[2]

# ╔═╡ 0018a2bc-1748-11eb-090f-8189d8958e3b
show_colored_kernel(edge_detection_kernel_vertical)

# ╔═╡ 25644f9e-1748-11eb-1785-b199708755e1
sum(edge_detection_kernel_vertical)

# ╔═╡ 64f95e42-1748-11eb-1048-7f168caea4b4
edge_enhanced_vertical = 3 * Gray.(abs.(convolve(image, edge_detection_kernel_vertical)))

# ╔═╡ bb0d4492-1748-11eb-05c6-0973e38a9f28
edge_enhanced_horizontal = 3 * Gray.(abs.(convolve(image, edge_detection_kernel_horizontal)))

# ╔═╡ c1823bbe-1745-11eb-10ad-a56948b9eaab
[image convolve(image, sharpen_kernel) convolve(convolve(image, sharpen_kernel), edge_detection_kernel_horizontal) convolve(convolve(image, sharpen_kernel), edge_detection_kernel_vertical) edge_enhanced_vertical edge_enhanced_horizontal]

# ╔═╡ 16d61752-1746-11eb-305e-81831983418b
sum(sharpen_kernel)

# ╔═╡ 0dad691e-174b-11eb-0970-4981859d5f82
md"""
## Trying Fourier Transforms
"""

# ╔═╡ 77d0acd4-173c-11eb-28ea-5d1566fbb940
begin
	function rgb_to_float(color)
    	return mean([color.r, color.g, color.b])
	end

	function fourier_spectrum_magnitudes(img)
		grey_values = rgb_to_float.(img)
		spectrum = fftshift(fft(grey_values))
		return abs.(spectrum)
	end

	function plot_1d_fourier_spectrum(img, dims=1)
		spectrum = fourier_spectrum_magnitudes(img)
		plot(centered(mean(spectrum, dims=1)[1:end]))
	end
end

# ╔═╡ ea23db4a-1745-11eb-07c4-5ffeafc03183
plot_1d_fourier_spectrum(image)

# ╔═╡ ea460db4-1745-11eb-1af5-1d206f5b617b
begin
	herd_zebras_url = "https://i.pinimg.com/originals/3c/66/74/3c6674c2c869cccdd741379fe593294d.jpg"
	download(herd_zebras_url, "herd_zebras.jpg")
	large_zebras = load("herd_zebras.jpg")
	shrink_zebras = shrink_image(large_zebras, 7)
end

# ╔═╡ 91ab8350-174c-11eb-2d29-01e1ef5116c9
plot_1d_fourier_spectrum(shrink_zebras)

# ╔═╡ ea635bb2-1745-11eb-2631-d56be9f6e88c
begin
	gauss_kernel = Kernel.gaussian((2, 2))
	conv_image = convolve(shrink_zebras, gauss_kernel)
end

# ╔═╡ ea80a71c-1745-11eb-3a0d-773e478efcc6
plot_1d_fourier_spectrum(conv_image)

# ╔═╡ ea9dec78-1745-11eb-2473-bbdce2c4c96a
begin
	function heatmap_2d_fourier_spectrum(img)
		heatmap(log.(fourier_spectrum_magnitudes(img)))
	end

	function heatmap_2d_fourier_spectrum(img)
		heatmap(log.(fourier_spectrum_magnitudes(img)))
	end
end

# ╔═╡ 77eb908a-173c-11eb-1323-8192758a3bf8
heatmap_2d_fourier_spectrum(shrink_zebras)

# ╔═╡ 780650a0-173c-11eb-00db-f7c0ec17cf52
heatmap_2d_fourier_spectrum(conv_image)

# ╔═╡ 78a80050-173c-11eb-2fe8-6b70bb6a2507

# ╔═╡ Cell order:
# ╟─41d37ef0-164b-11eb-0940-636023edbe0d
# ╠═3df9d1c0-164c-11eb-335f-edb4366d2625
# ╠═f8db5200-1680-11eb-1d10-cbc2b31fc024
# ╠═35aa24a0-1680-11eb-065b-7973a0dab607
# ╠═2a905470-1687-11eb-30c3-5f635c6c0aeb
# ╠═c2515c00-1687-11eb-3cf2-1385aa2ce3d6
# ╠═0c4eb462-1688-11eb-0f56-cf09d1f53f0a
# ╠═9c21a2b0-16af-11eb-3fae-6d8927359b72
# ╠═307c6610-16b6-11eb-159c-f18f672027df
# ╠═1592aef0-1688-11eb-2d74-35ce0cb1ddc3
# ╠═9650d1b0-16b6-11eb-061b-d94361398d7b
# ╠═765a986a-173c-11eb-29e6-cb155437d4ea
# ╠═7767a37e-173c-11eb-39a8-0d1bd7be9077
# ╠═7785d7e0-173c-11eb-11b4-ebc2a9c3815a
# ╠═cea7f490-1746-11eb-2a66-8f75edaec775
# ╠═b9733322-1747-11eb-1da7-3d4b992d60a2
# ╠═ed2204f0-1747-11eb-1662-9d8199f940a1
# ╠═0018a2bc-1748-11eb-090f-8189d8958e3b
# ╠═25644f9e-1748-11eb-1785-b199708755e1
# ╠═64f95e42-1748-11eb-1048-7f168caea4b4
# ╠═bb0d4492-1748-11eb-05c6-0973e38a9f28
# ╠═c1823bbe-1745-11eb-10ad-a56948b9eaab
# ╠═16d61752-1746-11eb-305e-81831983418b
# ╟─0dad691e-174b-11eb-0970-4981859d5f82
# ╠═77d0acd4-173c-11eb-28ea-5d1566fbb940
# ╠═ea23db4a-1745-11eb-07c4-5ffeafc03183
# ╠═ea460db4-1745-11eb-1af5-1d206f5b617b
# ╠═91ab8350-174c-11eb-2d29-01e1ef5116c9
# ╠═ea635bb2-1745-11eb-2631-d56be9f6e88c
# ╠═ea80a71c-1745-11eb-3a0d-773e478efcc6
# ╠═ea9dec78-1745-11eb-2473-bbdce2c4c96a
# ╠═77eb908a-173c-11eb-1323-8192758a3bf8
# ╠═780650a0-173c-11eb-00db-f7c0ec17cf52
# ╠═78a80050-173c-11eb-2fe8-6b70bb6a2507

max_stars_repo_name: martin-fabbri/colab-notebooks

max_issues_repo_name: martin-fabbri/colab-notebooks

max_forks_repo_name: martin-fabbri/colab-notebooks