text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

struct SMPSModel{N, T <: AbstractFloat, M <: AbstractMatrix}
    raw::RawSMPS{T}
    stages::NTuple{N, SMPSStage{T,M}}

    function SMPSModel(raw::RawSMPS{T}, stages::NTuple{N, SMPSStage{T,M}}) where {N, T <: AbstractFloat, M <: AbstractMatrix}
        return new{N,T,M}(raw, stages)
    end
end

function SMPSModel(raw::RawSMPS)
    N = length(raw.tim.stages)
    constructor = sparsity(raw.cor) >= 0.5 ? SparseLPData : LPData
    stages = ntuple(Val(N)) do i
        parse_stage(raw, i, N, constructor)
    end
    return SMPSModel(raw, stages)
end

function parse_stage(raw::RawSMPS, stage::Integer, N::Integer, constructor)
    1 <= stage <= N || error("$(stage.id) not in range 1 to $N.")
    # Get required data
    tim = raw.tim
    cor = raw.cor
    sto = raw.sto
    if stage == 1
        # Get row and column of the given stage
        row = tim.row_delims[stage]
        col = tim.col_delims[stage]
        # Compute the rowrange up to the next stage
        next_row = tim.row_delims[stage + 1]
        rowrange = collect(keys(cor.rows))[cor.rows[row][1]:cor.rows[next_row][1]-1]
        # Compute the colrange up to the next stage
        next_col = tim.col_delims[stage + 1]
        colrange = collect(keys(cor.vars))[cor.vars[col]:cor.vars[next_col]-1]
        # Extract first-stage model
        model = constructor(cor; colrange = colrange, rowrange = rowrange, include_constant = true)
        # Return stage
        return SMPSStage(model)
    else
        # Get row and column of the current stage
        row = tim.row_delims[stage]
        col = tim.col_delims[stage]
        # and the column of the previous stage
        prev_col = tim.col_delims[stage - 1]
        # Compute the rowrange up to the next stage
        rowrange = if stage < N
            next_row = tim.row_delims[stage + 1]
            rowrange = collect(keys(cor.rows))[cor.rows[row][1]:cor.rows[next_row][1]-1]
        else
            rowrange = collect(keys(cor.rows))[cor.rows[row][1]:end]
        end
        # Compute the colrange from the previos stage up to the next stage
        colrange = if stage < N
            next_col = tim.col_delims[stage + 1]
            colrange = collect(keys(cor.vars))[cor.vars[prev_col]:cor.vars[next_col]-1]
        else
            colrange = collect(keys(cor.vars))[cor.vars[prev_col]:end]
        end
        # Extract second-stage model
        model = constructor(cor; colrange = colrange, rowrange = rowrange)
        # Compute technology and recourse ranges
        delim      = findfirst(x -> x == col, colrange)
        technology = 1:delim-1
        recourse   = delim:length(colrange)
        # Extract uncertainty perturbation
        uncertainty = uncertainty_template(sto, cor, model, stage - 1)
        # Return stage
        return SMPSStage(stage, model, uncertainty, technology, recourse)
    end
end

function stochastic_model(smps::SMPSModel{2})
    # Create stochastic model
    sm = @stochastic_model begin
        @stage 1 begin
            @parameters begin
                data = smps.stages[1]
            end
            # Unpack first-stage SMPS data
            n    = data.model.indexmap.n
            m₂   = data.model.indexmap.m₂
            c₁   = data.model.c₁
            c₂   = data.model.c₂
            A    = data.model.A
            b    = data.model.b
            C, d = canonical(data.model.C,
                             data.model.d₁,
                             data.model.d₂)
            lb   = data.model.lb
            ub   = data.model.ub
            bin  = data.model.is_binary
            int  = data.model.is_integer
            # Define all first-stage variables as decisions
            if any(isfinite.(lb))
                if any(isfinite.(ub))
                    # Double bounds
                    @decision(model, lb[i] <= x[i in 1:n] <= ub[i])
                else
                    # Lower bounds
                    @decision(model, x[i in 1:n] >= lb[i])
                end
            elseif any(isfinite.(ub))
                # Upper bounds
                @decision(model, x[i in 1:n] <= ub[i])
            else
                # Free
                @decision(model, x[i in 1:n])
            end
            # Add any binary or integer restrictions
            for i in 1:n
                if bin[i]
                    @constraint(model, x[i] in MOI.ZeroOne())
                elseif int[i]
                    @constraint(model, x[i] in MOI.Integer())
                end
            end
            # Define objective and constraints
            if any(abs.(c₁) .> sqrt(eps())) || abs(c₂) > sqrt(eps())
                @objective(model, Min, dot(c₁, x) + c₂)
            end
            if length(A) > 0 && length(b) > 0
                @constraint(model, A * x .== b)
            end
            if length(C) > 0 && length(d) > 0
                @constraint(model, C * x .<= d)
            end
        end
        @stage 2 begin
            @parameters begin
                data = smps.stages[2]
            end
            @uncertain Δq ΔT ΔW Δh ΔC Δd₁ Δd₂ from SMPSScenario
            # Unpack first-stage SMPS data
            n    = length(data.recourse)
            m₂   = data.model.indexmap.m₂
            q    = data.model.c₁[data.recourse]
            Δq   = Δq[data.recourse]
            T    = data.model.A[:, data.technology]
            W    = data.model.A[:, data.recourse]
            h    = data.model.b
            C̃, h̃ = canonical(data.model.C + ΔC,
                             data.model.d₁ + Δd₁,
                             data.model.d₂ + Δd₂)
            T̃    = C̃[:, data.technology]
            W̃    = C̃[:, data.recourse]
            lb   = data.model.lb[data.recourse]
            ub   = data.model.ub[data.recourse]
            bin  = data.model.is_binary[data.recourse]
            int  = data.model.is_integer[data.recourse]
            # Define all recourse variables
            if any(isfinite.(lb))
                if any(isfinite.(ub))
                    # Double bounds
                    @recourse(model, lb[i] <= y[i in 1:n] <= ub[i])
                else
                    # Lower bounds
                    @recourse(model, y[i in 1:n] >= lb[i])
                end
            elseif any(isfinite.(ub))
                # Upper bounds
                @recourse(model, y[i in 1:n] <= ub[i])
            else
                # Free
                @recourse(model, y[i in 1:n])
            end
            # Add any binary or integer restrictions
            for i in 1:n
                if bin[i]
                    @constraint(model, y[i] in MOI.ZeroOne())
                elseif int[i]
                    @constraint(model, y[i] in MOI.Integer())
                end
            end
            # Define objective and constraints
            @objective(model, Min, dot((q + Δq), y))
            if (length(T) > 0 || length(W) > 0) && length(h) > 0
                @constraint(model, (T + ΔT) * x + (W + ΔW) * y .== (h + Δh))
            end
            if (length(T̃) > 0 || length(W̃) > 0) && length(h̃) > 0
                @constraint(model, T̃ * x + W̃ * y .<= h̃)
            end
        end
    end
    return sm
end


text:
abstract type AbstractFermionfields_2D{NC} <: AbstractFermionfields{NC,2}
end

function Base.setindex!(x::T,v,i1,i2,i5,i6)  where T <: AbstractFermionfields_2D
    @inbounds x.f[i1,i2 + x.NDW,i5 + x.NDW,i6] = v
end

function Base.getindex(x::T,i1,i2,i5,i6) where T <: AbstractFermionfields_2D
    #=
    i2new = i2 .+ x.NDW
    i3new = i3 .+ x.NDW
    i4new = i4 .+ x.NDW
    i5new = i5 .+ x.NDW
    @inbounds return x.f[i1,i2new,i3new,i4new,i5new,i6]
    =#
    @inbounds return x.f[i1,i2 .+ x.NDW,i5 .+ x.NDW,i6]
end

@inline function get_latticeindex_fermion(i,NC,NX,NT)
    #i =(((((ig-1)*NT+it-1)*NZ+iz-1)*NY+iy-1)*NX+ix-1)*NC+ic
    ic = (i-1) % NC + 1
    ii = div(i-ic,NC)
    #ii = (((ig-1)*NT+it-1)*NZ+iz-1)*NY+iy-1)*NX+ix-1
    ix = ii % NX + 1
    ii = div(ii-(ix-1),NX)
    #ii = ((ig-1)*NT+it-1)*NZ+iz-1)*NY+iy-1
    #ii = (ig-1)*NT+it-1
    it = ii % NT + 1
    ig = div(ii-(it-1),NT) + 1
    return ic,ix,it,ig
end

Base.length(x::T) where T <: AbstractFermionfields_2D = x.NC*x.NX*x.NT*x.NG

function Base.size(x::AbstractFermionfields_2D{NC})  where NC
    return (x.NC,x.NX,x.NT,x.NG)
end

function Base.iterate(x::T,state = 1) where T <: AbstractFermionfields_2D
    if state > length(x)
        return nothing
    end

    return (x[state],state+1)
end

function Base.setindex!(x::T,v,i)  where T <: AbstractFermionfields_2D
    ic,ix,it,ig  = get_latticeindex_fermion(i,x.NC,x.NX,x.NT)
    @inbounds x[ic,ix,it,ig] = v
end

function Base.getindex(x::T,i) where T <: AbstractFermionfields_2D
    ic,ix,it,ig  = get_latticeindex_fermion(i,x.NC,x.NX,x.NT)
    @inbounds return x[ic,ix,it,ig]
end

function Base.getindex(F::Adjoint_fermionfields{T},i1,i2,i5,i6) where T <: Abstractfermion  #F'
    @inbounds return conj(F.parent[i1,i2,i5,i6])
end

function Base.setindex!(F::Adjoint_fermionfields{T},v,i1,i2,i5,i6,μ)  where T <: Abstractfermion
    error("type $(typeof(F)) has no setindex method. This type is read only.")
end

function clear_fermion!(a::Vector{<: AbstractFermionfields_2D{NC}}) where NC
    for μ=1:4
        clear_fermion!(a[μ])
    end
end

function clear_fermion!(a::AbstractFermionfields_2D{NC}) where NC
    n1,n2,n5,n6 = size(a.f)
    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            a.f[i1,i2,i5,i6]= 0
                        end
                    end
                #end
            #end
        end
    end
end

function substitute_fermion!(a::AbstractFermionfields_2D{NC},b::AbstractFermionfields_2D{NC}) where NC
    n1,n2,n5,n6 = size(a.f)
    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            a.f[i1,i2,i5,i6]= b.f[i1,i2,i5,i6]
                        end
                    end
                #end
            #end
        end
    end
end

function substitute_fermion!(a::AbstractFermionfields_2D{NC},b::Abstractfermion) where NC
    NX = a.NX
    #NY = a.NY
    ##NZ = a.NZ
    NT = a.NT
    NG = a.NG
    @inbounds for i6=1:NG
        for i5=1:NT
            #for i4=1:NZ
                #for i3=1:NY
                    for i2=1:NX
                        @simd for i1=1:NC
                            a[i1,i2,i5,i6]= b[i1,i2,i5,i6]
                        end
                    end
                #end
            #end
        end
    end
    set_wing_fermion!(a)
end

struct Shifted_fermionfields_2D{NC,T} <: Shifted_fermionfields{NC,2}
    parent::T
    #parent::T
    shift::NTuple{2,Int8}
    NC::Int64

    #function Shifted_Gaugefields(U::T,shift,Dim) where {T <: AbstractGaugefields}
    function Shifted_fermionfields_2D(F::AbstractFermionfields_2D{NC},shift) where NC
        return new{NC,typeof(F)}(F,shift,NC)
    end
end

function Base.size(x::Shifted_fermionfields_2D)
    return size(x.parent)
end
using InteractiveUtils

#=
function shift_fermion(U::AbstractFermionfields_2D{NC},ν::T) where {T <: Integer,NC}
    return shift_fermion(U,Val(ν))
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{1}) where {T <: Integer,NC}
    shift = (1,0,0,0)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{2}) where {T <: Integer,NC}
    shift = (0,1,0,0)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{3}) where {T <: Integer,NC}
    shift = (0,0,1,0)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{4}) where {T <: Integer,NC}
    shift = (0,0,0,1)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{-1}) where {T <: Integer,NC}
    shift = (-1,0,0,0)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{-2}) where {T <: Integer,NC}
    shift = (0,-1,0,0)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{-3}) where {T <: Integer,NC}
    shift = (0,0,-1,0)
    return Shifted_fermionfields_2D(U,shift)
end

function shift_fermion(U::AbstractFermionfields_2D{NC},::Val{-4}) where {T <: Integer,NC}
    shift = (0,0,0,-1)
    return Shifted_fermionfields_2D(U,shift)
end
=#

        #lattice shift
function shift_fermion(F::AbstractFermionfields_2D{NC},ν::T) where {T <: Integer,NC}
    if ν == 1
        shift = (1,0)
    elseif ν == 2
        shift = (0,1)
    elseif ν == -1
            shift = (-1,0)
    elseif ν == -2
            shift = (0,-1)
    else
        error("ν = $ν")
    end

    return Shifted_fermionfields_2D(F,shift)
end

function shift_fermion(F::TF,shift::NTuple{Dim,T}) where {Dim,T <: Integer,TF <: AbstractFermionfields_2D}
    return Shifted_fermionfields_2D(F,shift)
end

function Base.setindex!(F::T,v,i1,i2,i5,i6)  where T <: Shifted_fermionfields_2D
    error("type $(typeof(F)) has no setindex method. This type is read only.")
end

function Base.getindex(F::T,i1,i2,i5,i6)  where T <: Shifted_fermionfields_2D
    @inbounds return F.parent[i1,i2.+ F.shift[1],i5.+ F.shift[2],i6]
end

function Base.getindex(F::T,i1::N,i2::N,i5::N,i6::N)  where {T <: Shifted_fermionfields_2D,N <: Integer}
    @inbounds return F.parent[i1,i2 + F.shift[1],i5 + F.shift[2],i6]
end

#=

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(1,0,0,0)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2.+ 1, i5,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(-1,0,0,0)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2 .- 1, i5,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(0,1,0,0)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2,i3 .+ 1,i4,i5,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(0,-1,0,0)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2,i3 .- 1,i4,i5,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(0,0,1,0)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2,i3 ,i4 .+ 1,i5,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(0,0,-1,0)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2,i3,i4 .- 1 ,i5,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(0,0,0,1)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2,i3 ,i4 ,i5 .+ 1,i6]
end

function Base.getindex(F::Shifted_fermionfields_2D{NC,T,Val{(0,0,0,-1)}},i1,i2, i5,i6)  where {NC,T,shift}
    @inbounds return F.parent[i1,i2,i3,i4 ,i5 .- 1,i6]
end
=#

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{NC},A::T,x::T3) where {NC,T<:Abstractfields,T3 <:Abstractfermion}
    #@assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    @assert NC != 3 "NC should not be 3"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for k1=1:NC
                            y[k1,ix,it,ialpha] = 0
                            @simd for k2=1:NC
                                y[k1,ix,it,ialpha] += A[k1,k2,ix,it]*x[k2,ix,it,ialpha]
                            end
                        end
                    end
                #end
            #end
        end
    end
end

@inline function updatefunc!(y,A,x,ix, it,ialpha)
    #@code_llvm  @inbounds x[1,ix, it,ialpha]
    #@code_typed x[1,ix, it,ialpha]
    #error("dd")
    #@code_lowered x[1,ix, it,ialpha]
    x1 = x[1,ix, it,ialpha]
    x2 = x[2,ix, it,ialpha]
    x3 = x[3,ix, it,ialpha]
    y[1,ix, it,ialpha] = A[1,1,ix, it]*x1 +
                                A[1,2,ix, it]*x2+
                                A[1,3,ix, it]*x3
    y[2,ix, it,ialpha] = A[2,1,ix, it]*x1+
                                A[2,2,ix, it]*x2 +
                                A[2,3,ix, it]*x3
    y[3,ix, it,ialpha] = A[3,1,ix, it]*x1+
                                A[3,2,ix, it]*x2 +
                                A[3,3,ix, it]*x3
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{3},A::T,x::T3) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 3 == x.NC "dimension mismatch! NC in y is 3 but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            #println("it = ",it, " ialpha = $ialpha")
            ##for iz=1:NZ
            #    #for iy=1:NY
                    for ix=1:NX
                        #updatefunc!(y,A,x,ix, it,ialpha)
                        #error("oo")
                        # #=
                        x1 = x[1,ix, it,ialpha]
                        x2 = x[2,ix, it,ialpha]
                        x3 = x[3,ix, it,ialpha]

                        y[1,ix, it,ialpha] = A[1,1,ix, it]*x1 +
                                                    A[1,2,ix, it]*x2+
                                                    A[1,3,ix, it]*x3
                        y[2,ix, it,ialpha] = A[2,1,ix, it]*x1+
                                                    A[2,2,ix, it]*x2 +
                                                    A[2,3,ix, it]*x3
                        y[3,ix, it,ialpha] = A[3,1,ix, it]*x1+
                                                    A[3,2,ix, it]*x2 +
                                                    A[3,3,ix, it]*x3
                        # =#
                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{3},A::T,x::T3,iseven::Bool) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 3 == x.NC "dimension mismatch! NC in y is 3 but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            #println("it = ",it, " ialpha = $ialpha")
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse((ix + iy + iz + it) % 2 == 0,true,false)
                        if evenodd == iseven
                        #updatefunc!(y,A,x,ix, it,ialpha)
                        #error("oo")
                        # #=
                            x1 = x[1,ix, it,ialpha]
                            x2 = x[2,ix, it,ialpha]
                            x3 = x[3,ix, it,ialpha]
                            y[1,ix, it,ialpha] = A[1,1,ix, it]*x1 +
                                                        A[1,2,ix, it]*x2+
                                                        A[1,3,ix, it]*x3
                            y[2,ix, it,ialpha] = A[2,1,ix, it]*x1+
                                                        A[2,2,ix, it]*x2 +
                                                        A[2,3,ix, it]*x3
                            y[3,ix, it,ialpha] = A[3,1,ix, it]*x1+
                                                        A[3,2,ix, it]*x2 +
                                                        A[3,3,ix, it]*x3
                        end
                        # =#
                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{2},A::T,x::T3) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 2 == x.NC "dimension mismatch! NC in y is 2 but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        x1 = x[1,ix, it,ialpha]
                        x2 = x[2,ix, it,ialpha]
                        y[1,ix, it,ialpha] = A[1,1,ix, it]*x1 +
                                                    A[1,2,ix, it]*x2
                        y[2,ix, it,ialpha] = A[2,1,ix, it]*x1+
                                                    A[2,2,ix, it]*x2

                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{2},A::T,x::T3,iseven::Bool) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 2 == x.NC "dimension mismatch! NC in y is 2 but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse((ix + it) % 2 == 0,true,false)
                        if evenodd == iseven
                            x1 = x[1,ix, it,ialpha]
                            x2 = x[2,ix, it,ialpha]
                            y[1,ix, it,ialpha] = A[1,1,ix, it]*x1 +
                                                        A[1,2,ix, it]*x2
                            y[2,ix, it,ialpha] = A[2,1,ix, it]*x1+
                                                        A[2,2,ix, it]*x2
                        end

                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{NC},x::T3,A::T) where {NC,T<:Abstractfields,T3 <:Abstractfermion}
    #@assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for k1=1:NC
                            y[k1,ix, it,ialpha] = 0
                            @simd for k2=1:NC
                                y[k1,ix, it,ialpha] += x[k1,ix, it,ialpha]*A[k1,k2,ix, it]
                            end
                        end
                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{NC},x::T3,A::T,iseven::Bool) where {NC,T<:Abstractfields,T3 <:Abstractfermion}
    #@assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        evenodd = ifelse((ix + it) % 2 == 0,true,false)
                        if evenodd == iseven

                            for k1=1:NC
                                y[k1,ix, it,ialpha] = 0
                                @simd for k2=1:NC
                                    y[k1,ix, it,ialpha] += x[k1,ix, it,ialpha]*A[k1,k2,ix, it]
                                end
                            end
                        end
                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{3},x::T3,A::T) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 3 == x.NC "dimension mismatch! NC in y is 3 but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        x1 = x[1,ix, it,ialpha]
                        x2 = x[2,ix, it,ialpha]
                        x3 = x[3,ix, it,ialpha]
                        y[1,ix, it,ialpha] = x1*A[1,1,ix, it] +
                                                    x2*A[2,1,ix, it]+
                                                    x3*A[3,1,ix, it]
                        y[2,ix, it,ialpha] = x1*A[1,2,ix, it]+
                                                    x2*A[2,2,ix, it] +
                                                    x3*A[3,2,ix, it]
                        y[3,ix, it,ialpha] = x1*A[1,3,ix, it]+
                                                    x2*A[2,3,ix, it] +
                                                    x3*A[3,3,ix, it]
                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{2},x::T3,A::T) where {T<:Abstractfields,T3 <:Abstractfermion}
    #@assert 2 == x.NC "dimension mismatch! NC in y is 2 but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        x1 = x[1,ix, it,ialpha]
                        x2 = x[2,ix, it,ialpha]
                        y[1,ix, it,ialpha] = x1*A[1,1,ix, it] +
                                                    x2*A[2,1,ix, it]
                        y[2,ix, it,ialpha] = x1*A[1,2,ix, it]+
                                                    x2*A[2,2,ix, it]

                    end
                #end
            #end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{NC},A::T,x::T3) where {NC,T<:Number,T3 <:Abstractfermion}
    @assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG

    @inbounds for ialpha=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for k1=1:NC
                            y[k1,ix, it,ialpha] = A*x[k1,ix, it,ialpha]
                        end
                    end
                #end
            #end
        end
    end
end

"""
mul!(u,x,y) -> u_{ab} = x_a*y_b
"""
function LinearAlgebra.mul!(u::T1,x::AbstractFermionfields_2D{NC},y::AbstractFermionfields_2D{NC}) where {T1 <: AbstractGaugefields,NC}
    NX = x.NX
    #NY = x.NY
    ##NZ = x.NZ
    NT = x.NT
    NG = x.NG
    clear_U!(u)

    for ik=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for ib=1:NC
                            @simd for ia=1:NC
                                u[ia,ib,ix, it] += x[ia,ix, it,ik]*y[ib,ix, it,ik]
                            end
                        end
                    end
                #end
            #end
        end
    end
    set_wing_U!(u)
end

function LinearAlgebra.mul!(u::T1,x::Abstractfermion,y::Adjoint_fermionfields{<: AbstractFermionfields_2D{NC}}) where {T1 <: AbstractGaugefields,NC}
    _,NX,NT,NG = size(y)
    clear_U!(u)

    for ik=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for ib=1:NC
                            @simd for ia=1:NC
                                u[ia,ib,ix, it] += x[ia,ix, it,ik]*y[ib,ix, it,ik]
                            end
                        end
                    end
                #end
            #end
        end
    end
    set_wing_U!(u)
end

function LinearAlgebra.mul!(u::T1,x::Adjoint_fermionfields{<: Shifted_fermionfields_2D{NC,T}},y::Abstractfermion) where {T1 <: AbstractGaugefields,NC,T}
    _,NX,NT,NG = size(x)
    clear_U!(u)

    for ik=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for ib=1:NC
                            @simd for ia=1:NC
                                u[ia,ib,ix, it] += x[ia,ix, it,ik]*y[ib,ix, it,ik]
                            end
                        end
                    end
                #end
            #end
        end
    end
    set_wing_U!(u)
end

function LinearAlgebra.mul!(u::T1,x::Adjoint_fermionfields{<: AbstractFermionfields_2D{NC}},y::Abstractfermion) where {T1 <: AbstractGaugefields,NC}
    _,NX,NT,NG = size(x)
    clear_U!(u)

    for ik=1:NG
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for ib=1:NC
                            @simd for ia=1:NC
                                u[ia,ib,ix, it] += x[ia,ix, it,ik]*y[ib,ix, it,ik]
                            end
                        end
                    end
                #end
            #end
        end
    end
    set_wing_U!(u)
end

function cross!(u::T1,x::Adjoint_fermionfields{<: Shifted_fermionfields_2D{NC,T}},y::Abstractfermion) where {T1 <: AbstractGaugefields,NC,T}
    mul!(u,y,x)
end

function cross!(u::T1,x::AbstractFermionfields_2D{NC},y::Abstractfermion) where {T1 <: AbstractGaugefields,NC}
    mul!(u,y,x)
end

function cross!(u::T1,x::Abstractfermion,y::Adjoint_fermionfields{<: AbstractFermionfields_2D{NC}}) where {T1 <: AbstractGaugefields,NC}
    mul!(u,y,x)
end

function cross!(u::T1,x::Adjoint_fermionfields{<: AbstractFermionfields_2D{NC}},y::Abstractfermion) where {T1 <: AbstractGaugefields,NC}
    mul!(u,y,x)
end

"""
mul!(y,A,x,α,β) -> α*A*x+β*y -> y
"""
function LinearAlgebra.mul!(y::AbstractFermionfields_2D{NC},A::T,x::T3,α::TA,β::TB) where {NC,T<:Number,T3 <:Abstractfermion,TA <: Number,TB <: Number}
    @assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG
    if A == one(A)
        @inbounds for ialpha=1:NG
            for it=1:NT
                ##for iz=1:NZ
                    ##for iy=1:NY
                        for ix=1:NX
                            for k1=1:NC
                                y[k1,ix, it,ialpha] = α*x[k1,ix, it,ialpha] + β*y[k1,ix, it,ialpha]
                            end
                        end
                    #end
                #end
            end
        end
    else
        @inbounds for ialpha=1:NG
            for it=1:NT
                ##for iz=1:NZ
                    ##for iy=1:NY
                        for ix=1:NX
                            for k1=1:NC
                                y[k1,ix, it,ialpha] = A*α*x[k1,ix, it,ialpha] + β*y[k1,ix, it,ialpha]
                            end
                        end
                    #end
                #end
            end
        end
    end
end

function LinearAlgebra.mul!(y::AbstractFermionfields_2D{NC},A::T,x::T3,α::TA,β::TB,iseven::Bool) where {NC,T<:Number,T3 <:Abstractfermion,TA <: Number,TB <: Number}
    @assert NC == x.NC "dimension mismatch! NC in y is $NC but NC in x is $(x.NC)"
    NX = y.NX
    #NY = y.NY
    ##NZ = y.NZ
    NT = y.NT
    NG = y.NG
    if A == one(A)
        @inbounds for ialpha=1:NG
            for it=1:NT
                ##for iz=1:NZ
                    ##for iy=1:NY
                        for ix=1:NX
                            evenodd = ifelse((ix + it ) % 2 == 0,true,false)
                            if evenodd == iseven
                                for k1=1:NC
                                    y[k1,ix, it,ialpha] = α*x[k1,ix, it,ialpha] + β*y[k1,ix, it,ialpha]
                                end
                            end
                        end
                    #end
                #end
            end
        end
    else
        @inbounds for ialpha=1:NG
            for it=1:NT
                ##for iz=1:NZ
                    ##for iy=1:NY
                        for ix=1:NX
                            evenodd = ifelse((ix + it ) % 2 == 0,true,false)
                            if evenodd == iseven
                                for k1=1:NC
                                    y[k1,ix, it,ialpha] = A*α*x[k1,ix, it,ialpha] + β*y[k1,ix, it,ialpha]
                                end
                            end
                        end
                    #end
                #end
            end
        end
    end
end

#Overwrite Y with X*a + Y*b, where a and b are scalars. Return Y.
function LinearAlgebra.axpby!(a::Number, X::T, b::Number, Y::AbstractFermionfields_2D{NC}) where {NC,T <: AbstractFermionfields_2D}
    n1,n2,n5,n6 = size(Y.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            Y.f[i1,i2,i5,i6] = a*X.f[i1,i2,i5,i6] + b*Y.f[i1,i2,i5,i6]
                        end
                    end
                #end
            #end
        end
    end
    return Y
end

function LinearAlgebra.axpy!(a::Number, X::T, Y::AbstractFermionfields_2D{NC}) where {NC,T <: AbstractFermionfields_2D}
    LinearAlgebra.axpby!(a,X,1,Y)
    return Y
end

function Base.:*(a::Number,x::AbstractFermionfields_2D{NC}) where {NC}
    y = similar(x)
    n1,n2,n5,n6 = size(y.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            y.f[i1,i2,i5,i6] = a*x.f[i1,i2,i5,i6]
                        end
                    end
                #end
            #end
        end
    end
    return y
end

function LinearAlgebra.rmul!(x::AbstractFermionfields_2D{NC},a::Number) where {NC}
    n1,n2,n5,n6 = size(x.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            x.f[i1,i2,i5,i6] = a*x.f[i1,i2,i5,i6]
                        end
                    end
                #end
            #end
        end
    end
    return x
end

function add_fermion!(c::AbstractFermionfields_2D{NC},α::Number,a::T1,β::Number,b::T2) where {NC,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b
    n1,n2,n5,n6 = size(c.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            #println(a.f[i1,i2, i5,i6],"\t",b.f[i1,i2, i5,i6] )
                            c.f[i1,i2,i5,i6] += α*a.f[i1,i2, i5,i6] + β*b.f[i1,i2, i5,i6]
                        end
                    end
                #end
            #end
        end
    end
    return
end

function add_fermion!(c::AbstractFermionfields_2D{NC},α::Number,a::T1) where {NC,T1 <: Abstractfermion}#c += alpha*a
    n1,n2,n5,n6 = size(c.f)

    @inbounds for i6=1:n6
        for i5=1:n5
            #for i4=1:n4
                #for i3=1:n3
                    for i2=1:n2
                        @simd for i1=1:NC
                            #println(a.f[i1,i2, i5,i6],"\t",b.f[i1,i2, i5,i6] )
                            c.f[i1,i2, i5,i6] += α*a.f[i1,i2, i5,i6]
                        end
                    end
                #end
            #end
        end
    end
    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::AbstractFermionfields_2D{NC}) where NC
    NX = x.NX
    #NY = x.NY
    ##NZ = x.NZ
    NT = x.NT
    n6 = size(x.f)[end]
    σ = sqrt(1/2)

    for ialpha = 1:n6
        for it=1:NT
            ##for iz=1:NZ
                ##for iy=1:NY
                    for ix=1:NX
                        for ic=1:NC
                            x[ic,ix, it,ialpha] = σ*randn()+im*σ*randn()
                        end
                    end
                #end
            #end
        end
    end
    set_wing_fermion!(x)
    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::AbstractFermionfields_2D{NC},randomfunc,σ) where NC

    NX = x.NX
    #NY = x.NY
    ##NZ = x.NZ
    NT = x.NT
    n6 = size(x.f)[end]
    #σ = sqrt(1/2)

    for mu = 1:n6
        for ic=1:NC
            for it=1:NT
                ##for iz=1:NZ
                    ##for iy=1:NY
                        for ix=1:NX
                            v1 = sqrt(-log(randomfunc()+1e-10))
                            v2 = 2pi*randomfunc()

                            xr = v1*cos(v2)
                            xi = v1 * sin(v2)

                            x[ic,ix, it,mu] = σ*xr + σ*im*xi
                        end
                    #end
                #end
            end
        end
    end

    set_wing_fermion!(x)

    return
end

function gauss_distribution_fermion!(x::AbstractFermionfields_2D{NC},randomfunc) where NC
    σ = 1
    gauss_distribution_fermion!(x,randomfunc,σ)
end

function Z2_distribution_fermion!(x::AbstractFermionfields_2D{NC})  where NC
    NX = x.NX
    #NY = x.NY
    ##NZ = x.NZ
    NT = x.NT
    n6 = size(x.f)[4]
    #σ = sqrt(1/2)

    for mu = 1:n6
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for ic=1:NC
                            x[ic,ix, it,mu] = rand([-1,1])
                        end
                    end
                #end
            #end
        end
    end

    set_wing_fermion!(x)

    return
end

"""
c-------------------------------------------------c
c     Random number function Z4  Noise
c     https://arxiv.org/pdf/1611.01193.pdf
c-------------------------------------------------c
"""
function Z4_distribution_fermi!(x::AbstractFermionfields_2D{NC})  where NC
    NX = x.NX
    #NY = x.NY
    #NZ = x.NZ
    NT = x.NT
    n6 = size(x.f)[4]
    θ = 0.0
    N::Int32 = 4
    Ninv = Float64(1/N)
    for ialpha = 1:n6
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @inbounds @simd for ic=1:NC
                            θ = Float64(rand(0:N-1))*π*Ninv # r \in [0,π/4,2π/4,3π/4]
                            x[ic,ix,iy,iz,it,ialpha] = cos(θ)+im*sin(θ)
                        end
                    end
                #end
            #end
        end
    end

    set_wing_fermion!(x)

    return
end

function uniform_distribution_fermion!(x::AbstractFermionfields_2D{NC})  where NC
    NX = x.NX
    #NY = x.NY
    #NZ = x.NZ
    NT = x.NT
    n6 = size(x.f)[6]
    #σ = sqrt(1/2)

    for mu = 1:n6
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for ic=1:NC
                            x[ic,ix, it,mu] = rand()*2 - 1 #each element should be in (-1,1)
                        end
                    end
                #end
            #end
        end
    end

    set_wing_fermion!(x)

    return
end

function fermion2vector!(vector,x::AbstractFermionfields_2D{NC}) where NC
    vector .= 0
    NX = x.NX
    #NY = x.NY
    #NZ = x.NZ
    NT = x.NT
    n6 = size(x.f)[6]
    #σ = sqrt(1/2)

    count = 0
    for mu = 1:n6
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        for ic=1:NC
                            count += 1
                            vector[count] = x[ic,ix, it,mu]
                        end
                    end
                #end
            #end
        end
    end

end

function fermions2vectors!(vector,x::Vector{<: AbstractFermionfields_2D{NC}}) where NC
    M = length(x)
    n,m = size(vector)
    @assert M == m

    vector .= 0
    NX = x[1].NX
    #NY = x[1].NY
    #NZ = x[1].NZ
    NT = x[1].NT
    n6 = size(x[1].f)[6]
    #σ = sqrt(1/2)

    for im = 1:M
        count = 0
        for mu = 1:n6
            for it=1:NT
                #for iz=1:NZ
                    #for iy=1:NY
                        for ix=1:NX
                            for ic=1:NC
                                count += 1
                                vector[count,im] = x[im][ic,ix, it,mu]
                            end
                        end
                    #end
                #end
            end
        end
    end

end

function LinearAlgebra.dot(a::AbstractFermionfields_2D{NC},b::AbstractFermionfields_2D{NC}) where {NC}
    NT = a.NT
    #NZ = a.NZ
    #NY = a.NY
    NX = a.NX
    NG = a.NG

    c = 0.0im
    @inbounds for α=1:NG
        for it=1:NT
            #for iz=1:NZ
                #for iy=1:NY
                    for ix=1:NX
                        @simd for ic=1:NC
                            c+= conj(a[ic,ix, it,α])*b[ic,ix, it,α]
                        end
                    end
                #end
            #end
        end
    end
    return c
end

ext: jl


text:
# converts BitVector to an integer.
function bitarr_to_int(arr, s = 0)
    v = 1
    for i in view(arr,length(arr):-1:1)
        s += v*i
        v <<= 1
    end
    s
end

# PART 1
function power_consumption(input)
  input_number_size = size(input, 2)
  num = falses(input_number_size)

  for i in 1:input_number_size
    col = input[:, i]
    bit = count(i->i==1, col) > length(col) / 2
    num[i] = bit
  end

  bitarr_to_int(num) * bitarr_to_int(.!num)
end

# PART 2
function life_support_rating(input)
  oxygen = filter_numbers(input, (col, numbers) -> count(i->i==1, col) >= size(numbers, 1) / 2)
  co2 = filter_numbers(input, (col, numbers) -> count(i->i==1, col) < size(numbers, 1) / 2)

  oxygen * co2
end

function filter_numbers(input, common_bit_fn)
  numbers = input

  for i in 1:size(input, 2)
    col = numbers[:, i]
    common_bit = common_bit_fn(col, numbers)
    numbers = numbers[(numbers[:, i] .== common_bit), :]
    if size(numbers, 1) == 1; break; end
  end

  bitarr_to_int(numbers[1, :])
end


text:
"""

    nestedhierarchymodel(S::Int64, L::Int64)

Return `UnipartiteNetwork` where resources are assigned to consumers according
to the nested hierarchy model for `S` species and `L`.

> Cattin, M.-F. et al. (2004) ‘Phylogenetic constraints and adaptation explain
> food-web structure’, Nature, 427(6977), pp. 835–839. doi: 10.1038/nature02327.

See also: `nichemodel`, `cascademodel`, `mpnmodel`
"""
function nestedhierarchymodel(S::Int64, L::Int64)

    # Evaluate input
    L >= S*S && throw(ArgumentError("Number of links L cannot be larger than the richness squared"))
    L <= 0 && throw(ArgumentError("Number of links L must be positive"))

    # Initiate matrix
    A = UnipartiteNetwork(zeros(Bool, (S, S)))

    Co = (L/(S*S))

    # Assign values to all species and sort in an ascending order
    e = sort(rand(S))

    # Beta parameter (after Allesina et al. 2008)
    # Classic niche: β = 1.0/(2.0*C)-1.0
    β = (S - 1.0)/(2.0*Co*S) - 1.0

    # Evaluate input
    β == 0 && throw(ArgumentError("β value is equal to zero! Try decreasing number of links"))

    # Random variable from the Beta distribution
    X = rand(Beta(1.0, β), 1)

    # 'r' values are assigned to all species
    r = e .* X

    # Number of prey for predator i, li rounded to the nearest integer
    expectedlinks = ((r ./ sum(r)) .* L)
    observedlinks = map(x -> round(Int, x), expectedlinks)

    totallinks = sum(observedlinks)

    # User message: number of links to be assigned
    # println("Nested hierarchy model with ", totallinks, " links to assign")

    # Links which may happen to exceed S are set to S-1
    for link in 1:S

        if observedlinks[link] >= S

            observedlinks[link] = S-1

        end

    end

    # Ensure that the first species has no links
    observedlinks[1] = 0

    # Species pool
    species_pool = species(A)

    # Create dictionaries with all the species indices and names

    # Strings into indices
    sp_pos = Dict{String, Int64}()

    for (n, f) in enumerate(species_pool)

        sp_pos[f] = n

    end

    # Indices into strings
    sp_str = Dict{Int64, String}()

    for (n, f) in enumerate(species_pool)

        sp_str[n] = f

    end

    # Species already predated
    predated = []

    # Assign links to all consumers except for the first one (basal)
    for consumer in 2:S

        # STAGE 1: Assign species with smaller niche value

        # Assign species randomly unless one has other predators
        eligible = sample(1:(consumer-1), (consumer-1), replace = false)

        linkstoassign = observedlinks[consumer]

        # Hold the last resource assagned in at this stage
        last_resource = 0

        # Loop through a random resources
        for resource in eligible

            # Update last resource holder
            last_resource = resource

            # Break the loop in case of reaching desired number of links
            linkstoassign == 0 && break

            # Assign link
            A[consumer, resource] = true
            linkstoassign = linkstoassign - 1

            # Update the predatored species array
            push!(predated, sp_str[resource])

            # If there are other predators break this
            (sum(A[:,resource]) > 1) && break

        end

        # Go to another consumer in case there are no more links
        linkstoassign == 0 && continue

        # STAGE 2: if the previos resource had a consumer assign resource from th eunion

        # Make a union of all consumers feeding on last_resource
        shared_consumers = A[:, species_pool[last_resource]]
        resource_union = Set{String}()
        [union!(resource_union, A[cons, :]) for cons in shared_consumers]

        # Remove species already utilized by current consumer
        resource_union = setdiff(resource_union, A[species_pool[consumer], :])

        # Shuffle the resource_union
        resource_union = sample(collect(resource_union), length(resource_union), replace = false)

        for resource in resource_union

            linkstoassign == 0 && break

            A[consumer, sp_pos[resource]] = true
            linkstoassign = linkstoassign - 1

            # Update the no_predator array
            push!(predated, resource)

        end

        # Go to another consumer in case there are no more links
        linkstoassign == 0 && continue

        # STAGE 3: randomly assign resources without predators

        # Remove already utilized species
        not_predated = setdiff(species_pool, predated)

        # Shuffle the current set
        not_predated = sample(collect(not_predated), length(not_predated), replace = false)

        # Assign links
        for resource in not_predated

            linkstoassign == 0 && break

            A[consumer, sp_pos[resource]] = true
            linkstoassign = linkstoassign - 1

        end

        # Go to another consumer in case there are no more links
        linkstoassign == 0 && continue

        # STAGE 4: Randomly assign species from the rest of the resources

        # Rest of the unassigned species!
        diet = collect(A[species(A)[consumer], :])

        # Get rid of specis already assigned to the current consumer
        not_consumed = setdiff(species_pool, diet)

        # Shuffle the set
        not_consumed = sample(collect(not_consumed), length(not_consumed), replace = false)

        # Assign links
        for resource in not_consumed

            linkstoassign == 0 && break

            A[consumer, sp_pos[resource]] = true
            linkstoassign = linkstoassign - 1

        end

    end

    return(A)

end

"""

    nestedhierarchymodel(S::Int64, Co::Float64)

Connectance can be provided instead of number of links.

See also: `nichemodel`, `cascademodel`, `mpnmodel`
"""
function nestedhierarchymodel(S::Int64, Co::Float64)

    # Calculate number of links L from connectednes
    L = round(Int, (Co*(S*S)))

    return nestedhierarchymodel(S, L)

end

"""

    nestedhierarchymodel(N::T) {T <: UnipartiteNetwork}

Applied to empirical `UnipartiteNetwork` return its randomized version.
"""
function nestedhierarchymodel(N::T) where {T <: UnipartiteNetwork}

    return nestedhierarchymodel(richness(N), connectance(N))

end

"""

    nestedhierarchymodel(parameters::Tuple)

Parameters tuple can also be provided in the form (Species::Int64, Co::Float64)
or (Species::Int64, Int::Int64).
"""
function nestedhierarchymodel(parameters::Tuple)

    return nestedhierarchymodel(parameters[1], parameters[2])

end


text:
function ncsize!(S, present, g, i, j)
	!present[i] && return 0
	S[i] != 0 && error("the graph is not acyclic")
	S[i] = 1
	for k in neighbors(g,i)
		if k != j && present[k]
			S[i] += ncsize!(S, present, g, k, i)
		end
	end
	return S[i]
end

"""
Optimally break tree until the giant component size is less than threshold.
"""
function tree_break!(present, g, threshold)
    S = zeros(Int, nv(g))
    H = MutableBinaryMaxHeap{Tuple{Int64,Int64}}()
    for i in vertices(g)
        if present[i] && S[i] == 0
            push!(H, (ncsize!(S, present, g, i, nothing), i))
        end
    end

    attack_nodes = Int[]
    scomps = Int[]
    while !isempty(H)
        scomp, i = pop!(H)
        sender = nothing
        while true
            sizes = [(S[k],k) for k in neighbors(g,i) if k != sender && present[k]]
            isempty(sizes) && break
            M, largest = maximum(sizes)
            if M <= scomp/2
                for k in neighbors(g,i)
                    if S[k] > 1 && present[k]
                        push!(H, (S[k],k))
                    end
                end
                present[i] = false
                push!(attack_nodes, i)
                push!(scomps, scomp)
                break
            end
            S[i] = 1
            for k in neighbors(g,i)
                if k != largest && present[k]
                    S[i] += S[k]
                end
            end
            sender, i = i, largest
        end
		scomp <= threshold && break
    end
    return attack_nodes, scomps
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function simulate_outage(;init_time_step, diesel_kw, fuel_available, b, m, diesel_min_turndown, batt_kwh, batt_kw,
                    batt_roundtrip_efficiency, n_timesteps, n_steps_per_hour, batt_soc_kwh, crit_load)
    """
    Determine how long the critical load can be met with gas generator and energy storage.
    :param init_time_step: Int, initial time step
    :param diesel_kw: float, generator capacity
    :param fuel_available: float, gallons
    :param b: float, diesel fuel burn rate intercept coefficient (y = m*x + b)  [gal/hr]
    :param m: float, diesel fuel burn rate slope (y = m*x + b)  [gal/kWh]
    :param diesel_min_turndown:
    :param batt_kwh: float, battery capacity
    :param batt_kw: float, battery inverter capacity (AC rating)
    :param batt_roundtrip_efficiency:
    :param batt_soc_kwh: float, battery state of charge in kWh
    :param n_timesteps: Int, number of time steps in a year
    :param n_steps_per_hour: Int, number of time steps per hour
    :param crit_load: list of float, load after DER (PV, Wind, ...)
    :return: float, number of hours that the critical load can be met using load following
    """
    for i in 0:n_timesteps-1
        t = (init_time_step - 1 + i) % n_timesteps + 1  # for wrapping around end of year
        load_kw = crit_load[t]

        if load_kw < 0  # load is met
            if batt_soc_kwh < batt_kwh  # charge battery if there's room in the battery
                batt_soc_kwh += minimum([
                    batt_kwh - batt_soc_kwh,     # room available
                    batt_kw / n_steps_per_hour * batt_roundtrip_efficiency,  # inverter capacity
                    -load_kw / n_steps_per_hour * batt_roundtrip_efficiency,  # excess energy
                ])
            end

        else  # check if we can meet load with generator then storage
            fuel_needed = (m * maximum([load_kw, diesel_min_turndown * diesel_kw]) + b) / n_steps_per_hour
            # (gal/kWh * kW + gal/hr) * hr = gal
            if load_kw <= diesel_kw && fuel_needed <= fuel_available  # diesel can meet load
                fuel_available -= fuel_needed
                if load_kw < diesel_min_turndown * diesel_kw  # extra generation goes to battery
                    if batt_soc_kwh < batt_kwh  # charge battery if there's room in the battery
                        batt_soc_kwh += minimum([
                            batt_kwh - batt_soc_kwh,     # room available
                            batt_kw / n_steps_per_hour * batt_roundtrip_efficiency,  # inverter capacity
                            (diesel_min_turndown * diesel_kw - load_kw) / n_steps_per_hour * batt_roundtrip_efficiency  # excess energy
                        ])
                    end
                end
                load_kw = 0

            else  # diesel can meet part or no load
                if fuel_needed > fuel_available && load_kw <= diesel_kw  # tank is limiting factor
                    load_kw -= maximum([0, (fuel_available * n_steps_per_hour - b) / m])  # (gal/hr - gal/hr) * kWh/gal = kW
                    fuel_available = 0

                elseif fuel_needed <= fuel_available && load_kw > diesel_kw  # diesel capacity is limiting factor
                    load_kw -= diesel_kw
                    # run diesel gen at max output
                    fuel_available = maximum([0, fuel_available - (diesel_kw * m + b) / n_steps_per_hour])
                                                                # (kW * gal/kWh + gal/hr) * hr = gal
                else  # fuel_needed > fuel_available && load_kw > diesel_kw  # limited by fuel and diesel capacity
                    # run diesel at full capacity and drain tank
                    load_kw -= minimum([diesel_kw, maximum([0, (fuel_available * n_steps_per_hour - b) / m])])
                    fuel_available = 0
                end

                if minimum([batt_kw, batt_soc_kwh * n_steps_per_hour]) >= load_kw  # battery can carry balance
                    # prevent battery charge from going negative
                    batt_soc_kwh = maximum([0, batt_soc_kwh - load_kw / n_steps_per_hour])
                    load_kw = 0
                end
            end
        end

        if round(load_kw, digits=5) > 0  # failed to meet load in this time step
            return i / n_steps_per_hour
        end
    end

    return n_timesteps / n_steps_per_hour  # met the critical load for all time steps
end

"""
    simulate_outages(;batt_kwh=0, batt_kw=0, pv_kw_ac_hourly=[], init_soc=0, critical_loads_kw=[], wind_kw_ac_hourly=[],
        batt_roundtrip_efficiency=0.829, diesel_kw=0, fuel_available=0, b=0, m=0, diesel_min_turndown=0.3,
    )

Time series simulation of outages starting at every time step of the year. Used to calculate how many time steps the
critical load can be met in every outage, which in turn is used to determine probabilities of meeting the critical load.

# Arguments
- batt_kwh: float, battery storage capacity
- batt_kw: float, battery inverter capacity
- pv_kw_ac_hourly: list of floats, AC production of PV system
- init_soc: list of floats between 0 and 1 inclusive, initial state-of-charge
- critical_loads_kw: list of floats
- wind_kw_ac_hourly: list of floats, AC production of wind turbine
- batt_roundtrip_efficiency: roundtrip battery efficiency
- diesel_kw: float, diesel generator capacity
- fuel_available: float, gallons of diesel fuel available
- b: float, diesel fuel burn rate intercept coefficient (y = m*x + b*rated_capacity)  [gal/kwh/kw]
- m: float, diesel fuel burn rate slope (y = m*x + b*rated_capacity)  [gal/kWh]
- diesel_min_turndown: minimum generator turndown in fraction of generator capacity (0 to 1)

Returns a dict
```
{
    "resilience_by_timestep": r,
    "resilience_hours_min": r_min,
    "resilience_hours_max": r_max,
    "resilience_hours_avg": r_avg,
    "outage_durations": x_vals,
    "probs_of_surviving": y_vals,
    "probs_of_surviving_by_month": y_vals_group_month,
    "probs_of_surviving_by_hour_of_the_day": y_vals_group_hour,
}
```
"""
function simulate_outages(;batt_kwh=0, batt_kw=0, pv_kw_ac_hourly=[], init_soc=[], critical_loads_kw=[], wind_kw_ac_hourly=[],
                     batt_roundtrip_efficiency=0.829, diesel_kw=0, fuel_available=0, b=0, m=0, diesel_min_turndown=0.3,
                     )
    n_timesteps = length(critical_loads_kw)
    n_steps_per_hour = Int(n_timesteps / 8760)
    r = repeat([0], n_timesteps)

    if batt_kw == 0 || batt_kwh == 0
        init_soc = repeat([0], n_timesteps)  # default is 0

        if (isempty(pv_kw_ac_hourly) || (sum(pv_kw_ac_hourly) == 0)) && diesel_kw == 0
            # no pv, generator, nor battery --> no resilience
            return Dict(
                "resilience_by_timestep" => r,
                "resilience_hours_min" => 0,
                "resilience_hours_max" => 0,
                "resilience_hours_avg" => 0,
                "outage_durations" => Int[],
                "probs_of_surviving" => Float64[],
            )
        end
    end

    if isempty(pv_kw_ac_hourly)
        pv_kw_ac_hourly = repeat([0], n_timesteps)
    end
    if isempty(wind_kw_ac_hourly)
        wind_kw_ac_hourly = repeat([0], n_timesteps)
    end
    load_minus_der = [ld - pv - wd for (pv, wd, ld) in zip(pv_kw_ac_hourly, wind_kw_ac_hourly, critical_loads_kw)]
    """
    Simulation starts here
    """
    # outer loop: do simulation starting at each time step

    for time_step in 1:n_timesteps
        r[time_step] = simulate_outage(;
            init_time_step = time_step,
            diesel_kw = diesel_kw,
            fuel_available = fuel_available,
            b = b, m = m,
            diesel_min_turndown = diesel_min_turndown,
            batt_kwh = batt_kwh,
            batt_kw = batt_kw,
            batt_roundtrip_efficiency = batt_roundtrip_efficiency,
            n_timesteps = n_timesteps,
            n_steps_per_hour = n_steps_per_hour,
            batt_soc_kwh = init_soc[time_step] * batt_kwh,
            crit_load = load_minus_der
        )
    end
    results = process_results(r, n_timesteps)
    return results
end

function process_results(r, n_timesteps)

    r_min = minimum(r)
    r_max = maximum(r)
    r_avg = round((float(sum(r)) / float(length(r))), digits=2)

    x_vals = collect(range(1, stop=Int(floor(r_max)+1)))
    y_vals = Array{Float64, 1}()

    for hrs in x_vals
        push!(y_vals, round(sum([h >= hrs ? 1 : 0 for h in r]) / n_timesteps,
                            digits=4))
    end
    return Dict(
        "resilience_by_timestep" => r,
        "resilience_hours_min" => r_min,
        "resilience_hours_max" => r_max,
        "resilience_hours_avg" => r_avg,
        "outage_durations" => x_vals,
        "probs_of_surviving" => y_vals,
    )
end

function simulate_outages(d::Dict, p::REoptInputs; microgrid_only::Bool=false)
    batt_roundtrip_efficiency = p.s.storage.charge_efficiency[:elec] *
                                p.s.storage.discharge_efficiency[:elec]

    # TODO handle generic PV names
    pv_kw_ac_hourly = zeros(length(p.time_steps))
    if "PV" in keys(d)
        pv_kw_ac_hourly = (
            get(d["PV"], "year_one_to_battery_series_kw", zeros(length(p.time_steps)))
          + get(d["PV"], "year_one_curtailed_production_series_kw", zeros(length(p.time_steps)))
          + get(d["PV"], "year_one_to_load_series_kw", zeros(length(p.time_steps)))
          + get(d["PV"], "year_one_to_grid_series_kw", zeros(length(p.time_steps)))
        )
    end
    if microgrid_only && !Bool(get(d, "PV_upgraded", false))
        pv_kw_ac_hourly = zeros(length(p.time_steps))
    end

    batt_kwh = 0
    batt_kw = 0
    init_soc = []
    if "Storage" in keys(d)
        batt_kwh = get(d["Storage"], "size_kwh", 0)
        batt_kw = get(d["Storage"], "size_kw", 0)
        init_soc = get(d["Storage"], "year_one_soc_series_pct", [])
    end
    if microgrid_only && !Bool(get(d, "storage_upgraded", false))
        batt_kwh = 0
        batt_kw = 0
        init_soc = []
    end

    diesel_kw = 0
    if "Generator" in keys(d)
        diesel_kw = get(d["Generator"], "size_kw", 0)
    end
    if microgrid_only
        diesel_kw = get(d, "Generator_mg_kw", 0)
    end

    simulate_outages(;
        batt_kwh = batt_kwh,
        batt_kw = batt_kw,
        pv_kw_ac_hourly = pv_kw_ac_hourly,
        init_soc = init_soc,
        critical_loads_kw = p.s.electric_load.critical_loads_kw,
        wind_kw_ac_hourly = [],
        batt_roundtrip_efficiency = batt_roundtrip_efficiency,
        diesel_kw = diesel_kw,
        fuel_available = p.s.generator.fuel_avail_gal,
        b = p.s.generator.fuel_intercept_gal_per_hr,
        m = p.s.generator.fuel_slope_gal_per_kwh,
        diesel_min_turndown = p.s.generator.min_turn_down_pct
    )
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
# For shared libraries in plugins
# This file only serves as a place to keep track of this constant
# Ahorn will not use it for anything besides loading/reloading plugin libraries

const loadedLibraries = String[]

max_stars_repo_name: CelestialCartographers/Ahorn

max_issues_repo_name: CelestialCartographers/Ahorn

max_forks_repo_name: CelestialCartographers/Ahorn


text:
module CollectiveSpins

export Spin, SpinCollection, CavityMode, CavitySpinCollection,
        geometry, GreenTensor, G_ij, Gamma_ij, Omega_ij, Omega_k_chain, Delta_k_chain, Omega_k_2D, Gamma_k_2D

include("system.jl")
include("timeevolution_base.jl")
include("geometry.jl")
include("interaction.jl")
include("effective_interaction.jl")
include("effective_interaction_simple.jl")
include("effective_interaction_rotated.jl")
include("quantum.jl")
include("reducedspin.jl")
include("independent.jl")
include("meanfield.jl")
include("mpc.jl")
include("io.jl")
include("collective_modes.jl")

using .system
using .quantum
using .meanfield
using .mpc
using .reducedspin
using .interaction
using .collective_modes

end # module


text:
function _gen_bs_docstring()
    call = WhenAt(
            Date(2020, 5, 2),
            Either(
                    Both(
                            Unit(Stock(:PETR4)),
                            Give(Scale(Konst(20.0), Unit(SpotCurrency(Currencies.BRL))))
                        ),
                    Worthless()
                )
        )

    put = WhenAt(
            Date(2020, 5, 2),
            Either(
                    Both(
                            Scale(Konst(20.0), Unit(SpotCurrency(Currencies.BRL))),
                            Give(Unit(Stock(:PETR4)))
                        ),
                    Worthless()
                )
        )

    return """
    BlackScholesModel

Pricing model for european call or put.

# Constructor:

```
BlackScholesModel(
    functional_currency::Currency,
    provider::MarketData.AbstractMarketDataProvider,
    underlying::Stock
)
```

# European Call contract example

```
$call
```

For a Call contract, the leg from the `Both` contract
which is different from `Unit{Stock}`
is recognized as the option *strike*.

# European Put contract example

```
$put
```

For a Put contract, the leg from the `Both` contract
which is not a `Give`
is recognized as the option *strike*.

# *strike* valuation

The *strike* will be priced using the `StaticHedgingModel`
with the same functional currency as the `BlackScholesModel`.
"""
end

@doc _gen_bs_docstring() BlackScholesModel


text:
global PROC
#global GPU
#global CPU

#if ~isdefined(:GPU)
#    GPU=false;
#    CPU=true
#end

if ~isdefined(:PROC)
    PROC="CPU"
end

#macro gpu(ex); if PROC=="GPU"; esc(ex); end; end
#macro cpu(ex); if ~GPU; esc(ex); end; end

macro gpu(ex)
    global PROC
    if PROC=="GPU" || PROC=="GPU32"
         esc(ex)
    end
end

macro cpu(ex)
    global PROC
    if PROC=="CPU"
         esc(ex)
    end
end

#export @gpuex, @cpuex

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module SubsetSelectionCIO

using MathProgBase, JuMP, Gurobi, CPLEX, Compat, LinearAlgebra
using DataFrames, CSV

import Compat.String

include("inner_op.jl")

export oa_formulation

getthreads() = haskey(ENV, "SLURM_JOB_CPUS_PER_NODE") ? parse(Int, ENV["SLURM_JOB_CPUS_PER_NODE"]) : 0

###########################
# FUNCTION oa_formulation
###########################
"""Computes the minimum regression error with Ridge regularization subject an explicit
cardinality constraint using cutting-planes.

w^* := arg min  ∑_i ℓ(y_i, x_i^T w) +1/(2γ) ||w||^2
           st.  ||w||_0 = k

INPUTS
  ℓ           - LossFunction to use
  Y           - Vector of outputs. For classification, use ±1 labels
  X           - Array of inputs
  k           - Sparsity parameter
  γ           - ℓ2-regularization parameter
  indices0    - (optional) Initial solution
  ΔT_max      - (optional) Maximum running time in seconds for the MIP solver. Default is 60s
  gap         - (optional) MIP solver accuracy

OUTPUT
  indices     - Indicates which features are used as regressors
  w           - Regression coefficients
  Δt          - Computational time (in seconds)
  status      - Solver status at termination
  Gap         - Optimality gap at termination
  cutCount    - Number of cuts needed in the cutting-plane algorithm
  """
function oa_formulation(ℓ::LossFunction, Y, X, k::Int, γ;
          indices0=findall(x-> x<k/size(X,2), rand(size(X,2))), ΔT_max=60, verbose=false, Gap=0e-3, solver::Symbol=:Gurobi)

  n,p = size(X)

  miop = (solver == :Gurobi) ?    Model(solver=GurobiSolver(MIPGap=Gap, TimeLimit=ΔT_max,
                                    OutputFlag=1*verbose, LazyConstraints=1, Threads=getthreads())) :
                                  Model(solver=CplexSolver(CPX_PARAM_EPGAP=Gap, CPX_PARAM_TILIM=ΔT_max,
                                    CPX_PARAM_SCRIND=1*verbose))
  s0 = zeros(p); s0[indices0].=1

  c0, ∇c0 = inner_op(ℓ, Y, X, s0, γ)

  # Optimization variables
  @variable(miop, s[j=1:p], Bin, start=s0[j])
  @variable(miop, t>=0, start=c0)

  # Objective
  @objective(miop, Min, t)

  # Constraints
  @constraint(miop, sum(s)<=k)

  cutCount=1; bestObj=c0; bestSolution=s0[:];
  @constraint(miop, t>= c0 + dot(∇c0, s-s0))

  #Info array
  bbdata = DataFrame(time=Float64[], node=Int[], obj=Float64[], bestbound=Float64[])

  # Outer approximation method for Convex Integer Optimization (CIO)
  function outer_approximation(cb)
    cutCount += 1
    c, ∇c = inner_op(ℓ, Y, X, getvalue(s), γ)
    if c<bestObj
      bestObj = c; bestSolution=getvalue(s)[:]
    end

    node      = MathProgBase.cbgetexplorednodes(cb)
    obj       = bestObj
    bestbound = MathProgBase.cbgetbestbound(cb)
    push!(bbdata, [time()-t0,node,obj,bestbound])

    @lazyconstraint(cb, t>=c + dot(∇c, s-getvalue(s)))
  end
  addlazycallback(miop, outer_approximation)

  t0=time()
  status = solve(miop)
  Δt = getsolvetime(miop)

  if status != :Optimal
    Gap = 1 - JuMP.getobjbound(miop) /  getobjectivevalue(miop)
  end
  if status == :Optimal
    bestSolution = getvalue(s)[:]
  end
  # Find selected regressors and run a standard linear regression with Tikhonov regularization
  indices = findall(s->s>0.5, bestSolution)
  w = SubsetSelection.recover_primal(ℓ, Y, X[:, indices], γ)

  return indices, w, Δt, status, Gap, cutCount
end

end # module


text:
using ReservoirComputing

#model parameters
const res_size = 30
const radius = 1.2
const activation = tanh
const degree = 6
const sparsity = 0.5
const sigma = 0.1
const beta = 0.0
const alpha = 1.0
const nla_type = NLADefault()
const in_size = 3
const extended_states = false
const h_steps = 2
const max_value = 0.8

const weight = 0.8
const fb_weight = 0.05
const jumps = 2

W_in = init_input_layer(res_size, in_size, sigma)

train_len = 50
predict_len = 12
data = ones(Float64, in_size, 100)
train = data[:, 1:1+train_len-1]
test = data[:, train_len:train_len+predict_len-1]

#test input layers functions
input_layer = [init_reservoir_givendeg, init_reservoir_givensp, DLR, DLRB, SCR, CRJ]

for t in input_layer

    if t == init_reservoir_givendeg
        W = t(res_size, radius, degree)
    elseif t == init_reservoir_givensp
        W = t(res_size, radius, sparsity)
    elseif t == DLRB
        W = t(res_size, weight, fb_weight)
    elseif t == CRJ
        W = t(res_size, weight, fb_weight, jumps)
    else
        W = t(res_size, weight)
    end
    esn = ESN(W,
        train,
        W_in,
        activation = activation,
        alpha = alpha,
        nla_type = nla_type,
        extended_states = extended_states)

    @test size(W, 1) == res_size
end

W = pseudoSVD(res_size, max_value, sparsity)
@test size(W, 1) == res_size


text:
using MPSProcess
using Test

@testset "MPSProcess.jl" begin
    # Write your tests here.
end


text:
"""
 Tests for PL1WM
"""

import Statistics # for mean values of trajectories
import LinearAlgebra # for the normn
using StochasticDiffEq
using Test
using Random
using DiffEqDevTools
#using DiffEqGPU

function prob_func(prob, i, repeat)
    remake(prob,seed=seeds[i])
end

"""
 Test Scalar SDEs (oop)
"""

@info "Scalar oop noise"

# PC exercise 14.2.2
numtraj = Int(1e4)
u₀ = 0.1
f(u,p,t) = p[1]*u
g(u,p,t) = p[2]*u
dts = 1 .//2 .^(6:-1:3)
tspan = (0.0,1.0)
p = [3//2,1//100]

h1(z) = z
#analytical_sol(t) = E(f(X(t))) =

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(f,g,u₀,tspan,p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end]),false),
        prob_func = prob_func
        )

sim = test_convergence(dts,ensemble_prob,PL1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.34 # order is 2.34
println("PL1WM:", sim.𝒪est[:weak_final])

"""
 Test Scalar SDEs (iip)
"""

@info "Scalar iip noise"

u₀ = [0.1]
f1!(du,u,p,t) =  (du[1] = p[1]*u[1])
g1!(du,u,p,t) =  (du[1] = p[2]*u[1])
dts = 1 .//2 .^(6:-1:3)
tspan = (0.0,1.0)
p = [3//2,1//100]

prob = SDEProblem(f1!,g1!,u₀,tspan,p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h1(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(1e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,PL1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.34 # order is 2.34
println("PL1WM:", sim.𝒪est[:weak_final])

"""
 Test non-commutative noise SDEs (iip)
"""

@info "Non-commutative noise"

u₀ = [1/8,1/8,1/1,1/8]
function f2!(du,u,p,t)
  du[1] = 243//154*u[1] - 27//77*u[2] + 23//154*u[3] - 65//154*u[4]
  du[2] = 27//77*u[1] - 243//154*u[2] + 65//154*u[3] - 23//154*u[4]
  du[3] = 5//154*u[1] - 61//154*u[2] + 162//77*u[3] - 36//77*u[4]
  du[4] = 61//154*u[1] - 5//154*u[2] + 36//77*u[3] - 162//77*u[4]
end
function g2!(du,u,p,t)
  du[1,1] = 1//9*sqrt(u[2]^2+u[3]^2+2//23)*1//13
  du[1,2] = 1//8*sqrt(u[4]^2+u[1]^2+1//11)*1//14
  du[1,3] = p[1]*1//12*sqrt(u[1]^2+u[2]^2+1//9)*1//6
  du[1,4] = p[1]*1//14*sqrt(u[3]^2+u[4]^2+3//29)*1//8
  du[2,1] = 1//9*sqrt(u[2]^2+u[3]^2+2//23)*1//14
  du[2,2] = 1//8*sqrt(u[4]^2+u[1]^2+1//11)*1//16
  du[2,3] = p[1]*1//12*sqrt(u[1]^2+u[2]^2+1//9)*1//5
  du[2,4] = p[1]*1//14*sqrt(u[3]^2+u[4]^2+3//29)*1//9
  du[3,1] = 1//9*sqrt(u[2]^2+u[3]^2+2//23)*1//13
  du[3,2] = 1//8*sqrt(u[4]^2+u[1]^2+1//11)*1//16
  du[3,3] = p[1]*1//12*sqrt(u[1]^2+u[2]^2+1//9)*1//5
  du[3,4] = p[1]*1//14*sqrt(u[3]^2+u[4]^2+3//29)*1//8
  du[4,1] = 1//9*sqrt(u[2]^2+u[3]^2+2//23)*1//15
  du[4,2] = 1//8*sqrt(u[4]^2+u[1]^2+1//11)*1//12
  du[4,3] = p[1]*1//12*sqrt(u[1]^2+u[2]^2+1//9)*1//6
  du[4,4] = p[1]*1//14*sqrt(u[3]^2+u[4]^2+3//29)*1//9
end
dts = 1 .//2 .^(6:-1:1)
tspan = (0.0,1.0)
p = [1]
h2(z) = z
# solution: E(X^i) = 1/8 exp(2*T), for i=1,2,4; E(X^3) = exp(2*T)

prob = SDEProblem(f2!,g2!,u₀,tspan,p,noise_rate_prototype=zeros(4,4))
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h2(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(1e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,PL1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀[1]*exp(2*1.0)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.33 # order is 1.6748033428458136
println("PL1WM:", sim.𝒪est[:weak_final])

"""
 Test Diagonal noise SDEs (iip), SIAM Journal on Numerical Analysis, 47 (2009), pp. 1713–1738
"""

@info "Diagonal noise"

u₀ = [0.1,0.1]
function f3!(du,u,p,t)
  du[1] = 3//2*u[1]
  du[2] = 3//2*u[2]
end
function g3!(du,u,p,t)
  du[1] = 1//10*u[1]
  du[2] = 1//10*u[2]
end
dts = 1 .//2 .^(3:-1:0)
tspan = (0.0,1.0)

h3(z) = z^2 # == 1//10**exp(3//2*t) if h3(z) = z and  == 1//100**exp(301//100*t) if h3(z) = z^2 )

prob = SDEProblem(f3!,g3!,u₀,tspan)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (h3(sol[end][1]),false),
        prob_func = prob_func
        )

numtraj = Int(5e4)
seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

sim = test_convergence(dts,ensemble_prob,PL1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=1//100*exp(301//100)
    )
@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("PL1WM:", sim.𝒪est[:weak_final])

"""
 Test Additive noise SDEs Kloeden & Platen Exercise 14.4.1
"""

numtraj = Int(1e3)
u₀ = 0.1
fadd(u,p,t) = p[1]*u
gadd(u,p,t) = p[2]
dts = 1 .//2 .^(4:-1:0)
tspan = (0.0,1.0)
p = [2//2,1//100]

seed = 100
Random.seed!(seed)
seeds = rand(UInt, numtraj)

prob = SDEProblem(fadd,gadd,u₀,tspan,p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (sol[end],false),
        prob_func = prob_func
        )

sim = test_convergence(dts,ensemble_prob,PL1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
    )

@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("PL1WM:", sim.𝒪est[:weak_final])

sim1 = test_convergence(dts,ensemble_prob,PL1WMA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
    )

@test abs(sim1.𝒪est[:weak_final]-2) < 0.3
println("PL1WMA:", sim1.𝒪est[:weak_final])

@test minimum(sim.solutions .≈ sim1.solutions)

#inplace

u₀ = [0.1]
fadd!(du,u,p,t) = @.(du = p[1]*u)
gadd!(du,u,p,t) = @.(du = p[2])

prob = SDEProblem(fadd!,gadd!,u₀,tspan,p)
ensemble_prob = EnsembleProblem(prob;
        output_func = (sol,i) -> (sol[end][1],false),
        prob_func = prob_func
        )

sim = test_convergence(dts,ensemble_prob,PL1WM(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
)

@test abs(sim.𝒪est[:weak_final]-2) < 0.3
println("PL1WM:", sim.𝒪est[:weak_final])

sim1 = test_convergence(dts,ensemble_prob,PL1WMA(),
    save_everystep=false,trajectories=numtraj,save_start=false,adaptive=false,
    weak_timeseries_errors=false,weak_dense_errors=false,
    expected_value=u₀.*exp(1.0*(p[1]))
    )

@test abs(sim1.𝒪est[:weak_final]-2) < 0.3
println("PL1WMA:", sim1.𝒪est[:weak_final])

@test minimum(sim.solutions .≈ sim1.solutions)


text:
# Copyright 2021 The FirstOrderLp Authors
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

function generate_primal_dual_hybrid_gradient_params(;
  l_inf_ruiz_iterations = 0,
  l2_norm_rescaling = false,
  pock_chambolle_alpha = nothing,
  iteration_limit = 200,
  primal_importance = 1.0,
  scale_invariant_initial_primal_weight = true,
  verbosity = 0,
  record_iteration_stats = true,
  restart_scheme = FirstOrderLp.NO_RESTARTS,
  restart_frequency_if_fixed = 100,
  artificial_restart_threshold = 0.5,
  use_weighted_average = true,
  sufficient_reduction_for_restart = 0.1,
  necessary_reduction_for_restart = 0.8,
  primal_weight_update_smoothing = 0.5,
  termination_evaluation_frequency = 5,
  use_approximate_localized_duality_gap = false,
  restart_to_current_metric = FirstOrderLp.GAP_OVER_DISTANCE_SQUARED,
  step_size_policy = "adaptive",
)
  if step_size_policy == "malitsky-pock"
    step_size_policy_params = FirstOrderLp.MalitskyPockStepsizeParameters(
      0.7,     # contraction_factor
      0.99,    # breaking_factor
      1.0,     # interpolation_coefficient
    )
  elseif step_size_policy == "constant"
    step_size_policy_params = FirstOrderLp.ConstantStepsizeParams()
  else
    step_size_policy_params = FirstOrderLp.AdaptiveStepsizeParams(
      0.3,     # reduction_exponent
      0.6,     # growth_exponent
    )
  end
  restart_params = FirstOrderLp.construct_restart_parameters(
    restart_scheme,
    restart_to_current_metric,
    restart_frequency_if_fixed,
    artificial_restart_threshold,
    sufficient_reduction_for_restart,
    necessary_reduction_for_restart,
    primal_weight_update_smoothing,
    use_approximate_localized_duality_gap,
  )
  parameters = FirstOrderLp.PdhgParameters(
    l_inf_ruiz_iterations,
    l2_norm_rescaling,
    pock_chambolle_alpha,
    primal_importance,
    scale_invariant_initial_primal_weight,
    verbosity,
    record_iteration_stats,
    termination_evaluation_frequency,
    terminate_on_iteration_limit(iteration_limit),
    restart_params,
    step_size_policy_params,
  )
  return parameters
end

@testset "Primal-dual hybrid gradient" begin
  @testset "Low precision" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 300,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-4
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-4
  end
  @testset "Terminate with optimal solution" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 1000,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_lp()
    parameters.termination_criteria.eps_optimal_absolute = 1e-8
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.termination_reason == FirstOrderLp.TERMINATION_REASON_OPTIMAL
  end
  @testset "Test Verbosity" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 300,
      primal_importance = 1.0,
      verbosity = 10,
    )
    problem = example_lp()
    old_stdout = stdout # We don't want to see this output on the screen.
    stdout_rd, stdout_wr = redirect_stdout()
    try
      output = FirstOrderLp.optimize(parameters, problem)
      @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-4
      @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-4
    finally
      redirect_stdout(old_stdout)
    end
  end
  @testset "Test Fixed Frequency Restart" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 500,
      primal_importance = 1.0,
      verbosity = 0,
      # This test breaks if you set restart_scheme=NO_RESTARTS.
      restart_scheme = FirstOrderLp.FIXED_FREQUENCY,
      restart_frequency_if_fixed = 30,
    )
    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9
  end
  @testset "Test Adaptive Restart Heuristic" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 600,
      primal_importance = 1.0,
      verbosity = 0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
    )

    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9

    @test any(
      it -> it.restart_used == FirstOrderLp.RESTART_CHOICE_RESTART_TO_AVERAGE,
      output.iteration_stats,
    )
  end

  @testset "Test adaptive_step_size=false, no smoothing" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 700,
      primal_importance = 1.0,
      verbosity = 0,
      primal_weight_update_smoothing = 0.0, # this test breaks if smoothing=0.5
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
      step_size_policy = "constant",
    )

    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9

    @test any(
      it -> it.restart_used == FirstOrderLp.RESTART_CHOICE_RESTART_TO_AVERAGE,
      output.iteration_stats,
    )
    step_size = output.iteration_stats[1].step_size
    for i in 2:length(output.iteration_stats)
      @test output.iteration_stats[i].step_size == step_size
    end
  end

  @testset "Test restart_to_current_metric = no_restart_to_current" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 600,
      primal_importance = 1.0,
      verbosity = 0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
      restart_to_current_metric = FirstOrderLp.NO_RESTART_TO_CURRENT,
    )

    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9

    @test any(
      it -> it.restart_used == FirstOrderLp.RESTART_CHOICE_RESTART_TO_AVERAGE,
      output.iteration_stats,
    )
  end

  @testset "Test restart_to_current_metric = gap_over_distance" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 600,
      primal_importance = 1.0,
      verbosity = 0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
      restart_to_current_metric = FirstOrderLp.GAP_OVER_DISTANCE,
    )

    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9

    @test any(
      it -> it.restart_used == FirstOrderLp.RESTART_CHOICE_RESTART_TO_AVERAGE,
      output.iteration_stats,
    )
  end

  @testset "Test Adaptive Restart" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 200,
      primal_importance = 1.0,
      verbosity = 0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
    )

    problem = example_lp()
    problem.objective_vector = [0.0, 0.0, 0.0, 0.0]
    parameters.termination_criteria.eps_optimal_absolute = 1e-8
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.termination_reason == FirstOrderLp.TERMINATION_REASON_OPTIMAL
  end

  @testset "use_approximate_localized_duality_gap = true" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 300,
      primal_importance = 1.0,
      verbosity = 0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
      use_approximate_localized_duality_gap = true,
    )

    problem = example_lp()
    problem.objective_vector = [0.0, 0.0, 0.0, 0.0]
    parameters.termination_criteria.eps_optimal_absolute = 1e-8
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.termination_reason == FirstOrderLp.TERMINATION_REASON_OPTIMAL
  end

  @testset "Malitsky and Pock linesearch, no smoothing" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 700,
      primal_importance = 1.0,
      primal_weight_update_smoothing = 0.0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
      step_size_policy = "malitsky-pock",
    )

    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9

  end

  @testset "Malitsky and Pock linesearch, smoothing" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 700,
      primal_importance = 1.0,
      restart_scheme = FirstOrderLp.ADAPTIVE_NORMALIZED,
      step_size_policy = "malitsky-pock",
    )

    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9

  end

  @testset "Quadratic Programming 1" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 200,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_qp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [0.2; 0.8] atol = 1.0e-4
    @test output.dual_solution ≈ [0.2] atol = 1.0e-4
  end
  @testset "Quadratic Programming 2" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 200,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_qp2()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [0.25; 0.0] atol = 1.0e-4
    @test output.dual_solution ≈ [0.0] atol = 1.0e-4
  end
  @testset "Test Preprocessing" begin
    @testset "l2 norm rescaling" begin
      parameters = generate_primal_dual_hybrid_gradient_params(
        l2_norm_rescaling = true,
        iteration_limit = 200,
        primal_importance = 1.0,
        verbosity = 0,
      )
      problem = example_qp2()
      output = FirstOrderLp.optimize(parameters, problem)
      @test output.primal_solution ≈ [0.25; 0.0] atol = 1.0e-4
      @test output.dual_solution ≈ [0.0] atol = 1.0e-4
    end
    @testset "ruiz" begin
      parameters = generate_primal_dual_hybrid_gradient_params(
        l_inf_ruiz_iterations = 10,
        iteration_limit = 200,
        primal_importance = 1.0,
        verbosity = 0,
      )
      problem = example_qp2()
      output = FirstOrderLp.optimize(parameters, problem)
      @test output.primal_solution ≈ [0.25; 0.0] atol = 1.0e-4
      @test output.dual_solution ≈ [0.0] atol = 1.0e-4
    end
    @testset "Pock-Chambolle rescaling" begin
      # TODO: Investigate why so many iterations are required.
      parameters = generate_primal_dual_hybrid_gradient_params(
        pock_chambolle_alpha = 1.0,
        iteration_limit = 3000,
        primal_importance = 1.0,
        verbosity = 0,
      )
      problem = example_lp()
      output = FirstOrderLp.optimize(parameters, problem)
      @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-4
      @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-4
    end
  end
  @testset "High precision" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 800,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [1.0; 0.0; 6.0; 2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [0.5; 4.0; 0.0] atol = 1.0e-9
  end
  @testset "Infeasible instance" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 800,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_lp()
    problem.right_hand_side[3] = 8
    output = FirstOrderLp.optimize(parameters, problem)
    final_stats = output.iteration_stats[end]
    @test output.termination_reason ==
          FirstOrderLp.TERMINATION_REASON_PRIMAL_INFEASIBLE
  end
  @testset "LP without bounds" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 400,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_lp_without_bounds()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [2.0] atol = 1.0e-9
    @test output.dual_solution ≈ [1.0] atol = 1.0e-9
  end
  @testset "Correlation Clustering: triangle plus" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 15,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_cc_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    tol = 1e-14
    final_stats = output.iteration_stats[end]
    @test output.primal_solution ≈ [1.0; 1.0; 0.0; 1.0; 0.0; 0.0] atol = tol
    @test final_stats.convergence_information[1].dual_objective ≈ 1.0 atol = tol
    # There are multiple optimal solutions.
    @test output.dual_solution[1] >= 0.0
    @test output.dual_solution[2] >= 0.0
    @test output.dual_solution[3] >= 0.0
    sum = output.dual_solution[1] + output.dual_solution[2]
    @test sum >= 1.0 - tol
  end
  @testset "Numerical error" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 150,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_cc_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    # Numerical error occurs because no convergence tolerances are set.
    @test output.termination_reason ==
          FirstOrderLp.TERMINATION_REASON_NUMERICAL_ERROR
    tol = 1e-14
    @test output.primal_solution ≈ [1.0; 1.0; 0.0; 1.0; 0.0; 0.0] atol = tol
    final_stats = output.iteration_stats[end]
    @test final_stats.convergence_information[1].dual_objective ≈ 1.0 atol = tol
    # There are multiple optimal solutions.
    @test output.dual_solution[1] >= 0.0
    @test output.dual_solution[2] >= 0.0
    @test output.dual_solution[3] >= 0.0
    sum = output.dual_solution[1] + output.dual_solution[2]
    @test sum >= 1.0 - tol
  end
  @testset "Correlation Clustering: star" begin
    parameters = generate_primal_dual_hybrid_gradient_params(
      iteration_limit = 100,
      primal_importance = 1.0,
      verbosity = 0,
    )
    problem = example_cc_star_lp()
    output = FirstOrderLp.optimize(parameters, problem)
    @test output.primal_solution ≈ [0.5; 0.5; 0.5; 0.0; 0.0; 0.0] atol = 1e-6
    @test output.dual_solution ≈ [0.5; 0.5; 0.5] atol = 1e-6
  end
end


text:
# Based on MATLAB code of Héctor Corte
# B.Sc. in physics 2010

include("new_solution.jl")

mutable struct SA <: AbstractParameters
    N::Int
    x_initial::Vector{Float64}
    tol_fun::Float64
    x::Vector{Float64}
    fx::Float64
end

"""

```julia
    SA(;
        x_initial::Vector = zeros(0),
        N::Int = 500,
        tol_fun::Real= 1e-4,
        information = Information(),
        options = Options()
    )
```

Parameters for the method of Simulated Annealing (Kirkpatrick et al., 1983).

Parameters:

- x_intial: Inital solution. If empty, then SA will generate a random one within the bounds.
- N: The number of test points per iteration.
- tol_fun: tolerance value for the Metropolis condition to accept or reject the test point as current point.

# Example

```jldoctest
julia> f(x) = sum(x.^2)
f (generic function with 1 method)

julia> optimize(f, [-1 -1 -1; 1 1 1.0], SA())
+=========== RESULT ==========+
  iteration: 60
    minimum: 5.0787e-68
  minimizer: [-2.2522059499734615e-34, 3.816133503985569e-36, 6.934348004465088e-36]
    f calls: 29002
 total time: 0.0943 s
+============================+

julia> optimize(f, [-1 -1 -1; 1 1 1.0], SA(N = 100, x_initial = [1, 0.5, -1]))
+=========== RESULT ==========+
  iteration: 300
    minimum: 1.99651e-69
  minimizer: [4.4638292404181215e-35, -1.738939846089388e-36, -9.542441152683457e-37]
    f calls: 29802
 total time: 0.0965 s
+============================+
```

"""
function SA(;
        x_initial::Vector = zeros(0),
        N::Int = 500,
        tol_fun::Real= 1e-4,
        information = Information(),
        options = Options()
    )

    parameters = SA(N, x_initial, tol_fun, zeros(0), Inf)

    Algorithm(
        parameters,
        information = information,
        options = options,
    )

end

function initialize!(
    status,
    parameters::SA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
   )

    l = view( problem.bounds, 1,:)
    u = view( problem.bounds, 2,:)

    if isempty(parameters.x_initial)
        parameters.x_initial = l .+ (u .- l) .* rand(length(u))
	end

    if options.f_calls_limit == 0
        options.f_calls_limit = 10000length(u)
    end

    if options.iterations == 0
        options.iterations = options.f_calls_limit ÷ parameters.N
    end

	# the current point and fx=f(x)
	x = parameters.x_initial
    best_sol = create_solution(x, problem)
	fx = best_sol.f
    status = State(best_sol, [best_sol])
    parameters.x = x
    parameters.fx = fx
    status.f_calls = 1

    return status

end

function update_state!(
    status::State{xf_indiv},
    parameters::SA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )

    nevals = status.f_calls
    max_evals = options.f_calls_limit
    N = parameters.N
    TolFun = parameters.tol_fun

    # T is the inverse of temperature.
    T = nevals / max_evals
    μ = 10.0 ^( 100T )

    l = view( problem.bounds, 1,:)
    u = view( problem.bounds, 2,:)

    # For each temperature we take 500 test points to simulate reach termal
    # equilibrium.
    for i = 1:N
        # We generate new test point using newSol function
        dx = newSol(2rand(length(parameters.x)) .- 1.0 , μ) .* (u-l)

        # the test point and fx1=f(x1)
        x1 = parameters.x + dx

        # Next step is to keep solution within bounds
        #x1 = (x1 .< l).*l+(l .<= x1).*(x1 .<= u).*x1+(u .< x1).*u
        reset_to_violated_bounds!(x1, problem.bounds)
        fx1 = evaluate(x1, problem)

        status.f_calls += 1

        df  = fx1 - parameters.fx

        # If the function variation,df, is <0 we take test point as current
        # point. And if df>0 we use Metropolis condition to accept or
        # reject the test point as current point.
        if (df < 0 || rand() < exp(-T*df/(abs(parameters.fx)) / TolFun))
            parameters.x = x1
            parameters.fx= fx1
        end

        # If the current point is better than current solution, we take
        # current point as cuyrrent solution.
        if fx1 < status.best_sol.f
            status.best_sol.x = x1
            status.best_sol.f = fx1
        end

        stop_criteria!(status, parameters, problem, information, options)

        if status.stop
            break
        end
    end

end

function final_stage!(
    status::State,
    parameters::SA,
    problem::AbstractProblem,
    information::Information,
    options::Options,
    args...;
    kargs...
    )
    status.final_time = time()
    status.population = [status.best_sol]
end


text:
function RosenbrockSchur!(V,dt,Fcn,Jac,CG,Param)
ROS=Param.ROS;
nV1=size(V,1);
nV2=size(V,2);
nV3=size(V,3);
nJ=nV1*nV2*nV3;
nStage=ROS.nStage;
k=Param.k
fV=Param.fV
Vn=Param.Vn

J = Param.J
JacSchur!(J,V,CG,Param)
Vn .= V
for iStage=1:nStage
  V .= Vn;
  for jStage=1:iStage-1
    @views V .= V .+ ROS.a[iStage,jStage] .* k[:,:,:,jStage];
  end
  Fcn(fV,V,CG,Param);
  for jStage=1:iStage-1
      @views fV .= fV .+ (ROS.c[iStage,jStage]/dt) .* k[:,:,:,jStage];
  end
  SchurSolve!(view(k,:,:,:,iStage),fV,J,dt*ROS.Gamma[iStage,iStage],Param);
  end
  V .= Vn;
  for iStage=1:nStage
    @views V .= V .+ ROS.m[iStage] .* k[:,:,:,iStage];
  end
end