text:
@testset "Poiseuille Problem" begin #for q in quadratures
    q = D2Q9()

    scale = 1
    ν = 1.0 / 6.0
    problem = PoiseuilleFlow(ν, scale, static = true)
    f_in = LatticeBoltzmann.initialize(AnalyticalEquilibrium(), q, problem)
    f_out = copy(f_in)
    collision_operator = CollisionModel(SRT, q, problem)

    feq = Array{Float64}(undef, size(f_in, 3))
    f = Array{Float64}(undef, size(f_in, 3))
    u = zeros(dimension(q))

    for x_idx in 2, y_idx in 1:size(f_in, 2)
        f = f_in[x_idx, y_idx, :]

        ρ = LatticeBoltzmann.density(q, f)
        LatticeBoltzmann.velocity!(q, f, ρ, u)
        T = LatticeBoltzmann.temperature(q, f, ρ, u)
    end

    t = 0
    Δt = LatticeBoltzmann.delta_t(problem)
    o = size(f_in, 2)
    p = size(f_in, 2) - 1
    LatticeBoltzmann.collide!(
        collision_operator,
        q,
        time = t * Δt,
        f_old = f_in,
        f_new = f_out,
    )
    LatticeBoltzmann.stream!(q, f_new = f_in, f_old = f_out)

    for x_idx in 2, y_idx in 1:size(f_in, 2)
        f = f_in[x_idx, y_idx, :]

        ρ = LatticeBoltzmann.density(q, f)
        LatticeBoltzmann.velocity!(q, f, ρ, u)
        T = temperature(q, f, ρ, u)
    end

    for x_idx in 2, y_idx in 1:size(f_in, 2)
        f = f_in[x_idx, y_idx, :]

        ρ = LatticeBoltzmann.density(q, f)
        LatticeBoltzmann.velocity!(q, f, ρ, u)
        T = temperature(q, f, ρ, u)
    end

    # TODO: when an a thermal model is considered, the density fluctuates
    # TODO: when a thermal model is considerd we have constant density
end


text:
export update_y_l_parallel

"""
  update l and compute y . This is a subfunction for PARSDMM.jl
"""
function update_y_l_parallel(
  x                     ::Vector{TF},
  i                     ::Integer,
  Blas_active           ::Bool,
  y                     ::Vector{Vector{TF}},
  y_old                 ::Vector{Vector{TF}},
  l                     ::Vector{Vector{TF}},
  l_old                 ::Vector{Vector{TF}},
  rho                   ::Vector{TF},
  gamma                 ::Vector{TF},
  prox                  ::Vector{Any},
  TD_OP                 ::Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joLinearFunction{TF,TF}}},
  P_sub,
  x_hat                 ::Vector{Vector{TF}},
  r_pri                 ::Vector{Vector{TF}},
  s                     ::Vector{Vector{TF}},
  set_feas              ::Vector{TF},
  feasibility_only=false::Bool
  ) where {TF<:Real,TI<:Integer}

rho1 = TF
rho1 = TF(1.0) ./ rho[1];

  copy!(y_old[1],y[1]);
  copy!(l_old[1],l[1]);

  s[1] = TD_OP[1]*x;
  if Blas_active
    if gamma[1]==1 #without relaxation
      copyto!(y[1],s[1])
      #y[1]       = prox[1]( BLAS.axpy!(-rho1[1],l[1],y[1]) )
      BLAS.axpy!(-rho1[1],l[1],y[1])
      y[1] = prox[1](y[1])
      @. r_pri[1] = -s[1]+y[1];
      BLAS.axpy!(rho[1],r_pri[1],l[1]);
    else #relaxed iterations
      @. x_hat[1] = gamma[1]*s[1]
      BLAS.axpy!(TF(1.0)-gamma[1],y[1],x_hat[1]);
      #y[1]       = copy(x_hat[1]);
      copyto!(y[1],x_hat[1])
      #y[1]       = prox[1]( BLAS.axpy!(-rho1[1],l[1],y[1]))
      BLAS.axpy!(-rho1[1],l[1],y[1])
      y[1] = prox[1](y[1])
       @. r_pri[1] = -s[1]+y[1];
      BLAS.axpy!(rho[1],y[1]-x_hat[1],l[1]);
    end
  else
    if gamma[1]==1 #without relaxation
      #y[1]       = prox[1]( s[1]-l[1]*rho1[1] );
       @. y[1]       = s[1]-l[1]*rho1[1]
       y[1] = prox[1](y[1]);
       @. r_pri[1]   = -s[1]+y[1];
       @. l[1]       = l[1]+rho[1]*r_pri[1];
    else #relaxed iterations
       @. x_hat[1]   = gamma[1]*s[1] + ( TF(1.0)-gamma[1] )*y[1]
      #y[1]       = prox[1]( x_hat[1]-l[1]*rho1[1] );
       @. y[1]       = x_hat[1]-l[1]*rho1[1]
       y[1] = prox[1]( y[1] )
       @. r_pri[1]   = -s[1]+y[1]
       @. l[1]       = l[1]+rho[1]*( -x_hat[1]+y[1] )
    end
  end #end blas

#log_PARSDMM.r_pri[i,1]  = norm(r_pri[1])#./(norm(y[1])+(100*eps(TF))); #add 1-14, because y may be 0 for certain initializations of x,y,l
#log r_pri outside this function for parallel use (temporarily)
# if Blas_active
#   log_PARSDMM.r_dual[i,1] = norm( BLAS.scal!(lx,rho[1],TD_OP[1]'*(y[1]-y_old[1]),1 ))
# else
#  log_PARSDMM.r_dual[i,1] = norm( rho[1]*(TD_OP[1]'*(y[1]-y_old[1])) )
# end
#log feasibility
if feasibility_only==false
  if mod(i,10)==0 && myid()<nprocs() #ii<length(P_sub)#log every 10 it, or whatever number is suitable
     temp = P_sub[1](deepcopy(s[1]))
     set_feas[1] = norm(temp-s[1]) ./ (norm(s[1])+(100*eps(TF)))
  end
else
  if mod(i,10)==0 #ii<length(P_sub)#log every 10 it, or whatever number is suitable
     temp = P_sub[1](deepcopy(s[1]))
     set_feas[1] = norm(temp-s[1]) ./ (norm(s[1])+(100*eps(TF)))
  end
end

# if mod(i,10)==0
#   counter+=1;
# end

#return y,l,r_pri,s,log_PARSDMM,counter,y_old,l_old,x_hat,set_feas_vec

end #end of function

ext: jl


text:
include("../psmsdst.jl")

"""
    solvegreedy(problem)

Solves the given PSMSDST problem using a greedy heuristic: for every pair (j,m)
of available jobs and machines, the pair that produces the lowest increment
in the objective function is programmed, until a complete solution is reached.
"""
function solvegreedy(problem::PSMSDSTProblem)::Vector{Vector{Int}}

    machines = 1:problem.nummachines

    # Create resulting vector for the solution
    solution = [ Vector{Int}() for machineid in machines ]

    # Create a set for the jobs that have not been programmed yet
    availablejobs = Set{Int}(1:problem.numjobs)

    # Create a vector for the finish times of each machine
    machinefinishtimes = [ 0 for machineid in machines ]

    # Variable that holds the current objective function
    currentmaxfinishtime = maximum(machinefinishtimes)

    # Matrix that contains how the machine finish time would increase when
    # a job is programmed on it.
    delays = [ problem.durations[job] for job in 1:problem.numjobs,
                machine in machines ]

    # Matrix that indicates if programming a job on a machine would invalidate
    # the pre-calculated delay.
    createsdelayupdate = [ true for job in 1:problem.numjobs,
        machine in machines ]

    # The id of the machine that will be updated because the cached delays
    # where invalidated.
    updatedelays = 0

    # Program jobs until a complete solution is reached
    while !isempty(availablejobs)

        # Cartesian product, (j,m) pairs of jobs and machines.
        options = Iterators.product(availablejobs,machines)

        if updatedelays > 0
            # The delays of a machine need to be recalculated.
            machine = updatedelays
            for job in availablejobs

                # The initial delay is just the duration of the job
                delay = problem.durations[job]

                if length(solution[machine]) > 0
                    lastjob = solution[machine][end]
                    lastfamily = problem.families[lastjob]
                    thisfamily = problem.families[job]

                    # If the family of the next job is different to the
                    # family of this job, this means that all the delays are
                    # invalidated because the setup times will change.
                    createsdelayupdate[job,machine] = thisfamily != lastfamily

                    # Add the setup time to the delay
                    delay += problem.setuptimes[lastfamily,thisfamily]
                end

                delays[job,machine] = delay
            end
        end

        # Look for the best option
        bestoption = nothing
        bestcost = typemax(Int64)

        for (job,machine) in options
            # Finish time of the machine if the job was programmed on it.
            virtualmachinefinish = (machinefinishtimes[machine]
                +delays[job,machine])

            # Objective function if the job was programmed on it.
            cost = max(currentmaxfinishtime,virtualmachinefinish)

            if cost < bestcost
                bestoption = (job,machine,virtualmachinefinish,
                    createsdelayupdate[job,machine])
                bestcost = cost
            end
        end

        job,machine,finishtime,createsupdate = bestoption

        # If the job creates an update in the delays, they have to be
        # re-calculated in the next iteration.
        if createsupdate
            updatedelays = machine
        else
            updatedelays = 0
        end

        machinefinishtimes[machine] = finishtime

        currentmaxfinishtime = maximum(machinefinishtimes)
        push!(solution[machine],job)
        delete!(availablejobs,job)
    end

    return solution
end

#
#    getsolver()
#
#Returns a pair containing the name of the solver, and a pre-parameterized
#function that only takes the problem as argument.
#
#Internally it also solves a small problem using the solver, ensuring that
#the code for the solver is compiled ahead of time and there is no reduced
#performance during evaluation.
#
## Examples
#```julia-repl
#julia> include("solvers/solver.jl")(param1=42.0,param2=64.0)
#"Solver param1=42 param2=64"=>#1 (generic function with 1 method)
#```

function getsolver()::Pair{String,Function}
    solver = problem -> solvegreedy(problem)

    presolve(solver)

    return "Greedy"=>solver
end

max_stars_repo_name: mkualquiera/huepl

max_issues_repo_name: mkualquiera/huepl

max_forks_repo_name: mkualquiera/huepl


text:
# The HORIZONS.jl package is licensed under the MIT "Expat" License:
# Copyright (c) 2017: Jorge Pérez.

module HORIZONS

using Expect, FTPClient

using DelimitedFiles, Dates

export HORIZONS_MACHINE, horizons, vec_tbl, vec_tbl_csv,
    smb_spk, smb_spk_ele

const HORIZONS_MACHINE = "ssd.jpl.nasa.gov"
const HORIZONS_FTP_DIR = "pub/ssd/"
const HORIZONS_DATE_FORMAT = "yyyy-u-dd HH:MM:SS.sss"
const ObjectName = Union{Int, String}
const StartStopTime = Union{DateTime, Date, String}
const StepSize = Union{Period, Int, String}
const VecTable = Union{Int, String}

"""
`horizons()`

Connect to JPL HORIZONS telnet interface

`telnet horizons.jpl.nasa.gov 6775`

"""
function horizons()
run(ignorestatus(`telnet horizons.jpl.nasa.gov 6775`))
end

#auxiliary function which translates a ::Bool to a "YES" or "NO" string
function yesornostring(yesorno::Bool)
    yesorno ? "YES" : "NO"
end

include("vec_tbl.jl")
include("smb_spk.jl")
include("smb_spk_ele.jl")

end # module


text:
# This is basically maple/lda_c_vwn.mpl and maple/vwn.mpl from libxc
#
# LDA correlation according to Vosko Wilk,and Nusair, (DOI 10.1139/p80-159)
#
# This version generates the "standard" parametrization, where a Pade approximation of a
# numerical diffusion Monte-Carlo simulation is used.

as_rational(x) = nsimplify(Sym(x))

A_vwn  = as_rational.([ 0.0310907, 0.01554535, -1/(6*PI^2)])
b_vwn  = as_rational.([ 3.72744,   7.06042,    1.13107  ])
c_vwn  = as_rational.([12.9352,   18.0578,    13.0045   ])
x0_vwn = as_rational.([-0.10498,  -0.32500,   -0.0047584])

A_rpa  = as_rational.([ 0.0310907,  0.01554535,  -1/(6*PI^2)])
b_rpa  = as_rational.([13.0720,    20.1231,      1.06835  ])
c_rpa  = as_rational.([42.7198,   101.578,      11.4813   ])
x0_rpa = as_rational.([-0.409286,  -0.743294,   -0.228344 ])

Q(b, c)      = sqrt(4*c - b^2)
f1(b, c)     = 2*b/Q(b, c)
f2(b, c, x0) = b*x0/(x0^2 + b*x0 + c)
f3(b, c, x0) = 2*(2x0 + b)/Q(b, c)

fpp = 4/(9*(2^(1//3) - 1))

fx(b, c, rs) = rs + b*sqrt(rs) + c

f_aux(A, b, c, x0, rs) = A*(
  + log(rs/fx(b, c, rs))
  + (f1(b, c) - f2(b, c, x0)*f3(b, c, x0))*atan(Q(b, c)/(2*sqrt(rs) + b))
  - f2(b, c, x0)*log((sqrt(rs) - x0)^2/fx(b, c, rs))
)

DMC(rs, z) = (
    + f_aux(A_vwn[2], b_vwn[2], c_vwn[2], x0_vwn[2], rs)
    - f_aux(A_vwn[1], b_vwn[1], c_vwn[1], x0_vwn[1], rs)
)

DRPA(rs, z) = (
    + f_aux(A_rpa[2], b_rpa[2], c_rpa[2], x0_rpa[2], rs)
    - f_aux(A_rpa[1], b_rpa[1], c_rpa[1], x0_rpa[1], rs)
)

f_vwn(rs, z) = (
  + f_aux(A_vwn[1], b_vwn[1], c_vwn[1], x0_vwn[1], rs)
  + f_aux(A_vwn[3], b_vwn[3], c_vwn[3], x0_vwn[3], rs)*f_zeta(z)*(1 - z^4)/fpp
  +  DMC(rs, z)*f_zeta(z)*z^4
)

f(rs, z, junk...) = f_vwn(rs, z)


text:
using PATHSolver
using SparseArrays
using Test

@testset "CheckLicense" begin
    @test PATHSolver.c_api_License_SetString("bad_license") != 0
    @test PATHSolver.c_api_Path_CheckLicense(1, 1) > 0
end

@testset "PathVersion" begin
    s = PATHSolver.c_api_Path_Version()
    @test match(r"Path [0-9]\.[0-9]\.[0-9][0-9]", s) !== nothing
end

@testset "Example" begin
    M = convert(
        SparseArrays.SparseMatrixCSC{Cdouble, Cint},
        SparseArrays.sparse([
            0  0 -1 -1;
            0  0  1 -2;
            1 -1  2 -2;
            1  2 -2  4
        ])
    )
    status, z, info = PATHSolver.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        [0.0, 0.0, 0.0, 0.0];
        output = "yes",
    )
    @test status == PATHSolver.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end

@testset "Example LUSOL" begin
    M = convert(
        SparseArrays.SparseMatrixCSC{Cdouble, Cint},
        SparseArrays.sparse([
            0  0 -1 -1;
            0  0  1 -2;
            1 -1  2 -2;
            1  2 -2  4
        ])
    )
    status, z, info = PATHSolver.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        [0.0, 0.0, 0.0, 0.0];
        output = "yes",
        # TODO(odow): when enabled, I get segfaults :(
        factorization_method = "blu_lusol",
        factorization_library_name = PATHSolver.LUSOL_LIBRARY_PATH,
    )
    @test status == PATHSolver.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end

@testset "Example II" begin
    M = [
        0  0 -1 -1;
        0  0  1 -2;
        1 -1  2 -2;
        1  2 -2  4
    ]

    q = [2; 2; -2; -6]

    function F(n::Cint, x::Vector{Cdouble}, f::Vector{Cdouble})
        @assert n == length(x) == length(f) == 4
        f[1] = -x[3]^2 - x[4] + q[1]
        f[2] =  x[3]^3 - 2x[4]^2 + q[2]
        f[3] = x[1]^5 - x[2] + 2x[3] - 2x[4] + q[3]
        f[4] = x[1] + 2x[2]^3 - 2x[3] + 4x[4] + q[4]
        return Cint(0)
    end

    function J(
        n::Cint,
        nnz::Cint,
        x::Vector{Cdouble},
        col::Vector{Cint},
        len::Vector{Cint},
        row::Vector{Cint},
        data::Vector{Cdouble}
    )
        JAC = [
            0       0     -2x[3]      -1;
            0       0    3x[3]^2  -4x[4];
            5x[1]^4      -1          2      -2;
            1   6x[2]         -2       4
        ]
        @assert n == length(x) == length(col) == length(len) == 4
        @assert nnz == length(row) == length(data)
        i = 1
        for c in 1:n
            col[c] = i
            len[c] = 0
            for r in 1:n
                if !iszero(JAC[r, c])
                    data[i] = JAC[r, c]
                    row[i] = r
                    len[c] += 1
                    i += 1
                end
            end
        end
        return Cint(0)
    end

    status, z, info = PATHSolver.solve_mcp(
        F,
        J,
        fill(0.0, 4),
        fill(10.0, 4),
        [1.0, 1.0, 1.0, 1.0];
        output = "yes"
    )
    @test status == PATHSolver.MCP_Solved
    @test isapprox(z, [1.28475, 0.972916, 0.909376, 1.17304], atol=1e-4)
end

@testset "Name Test" begin
    M = convert(
        SparseArrays.SparseMatrixCSC{Cdouble, Cint},
        SparseArrays.sparse([
            0  0 -1 -1;
            0  0  1 -2;
            1 -1  2 -2;
            1  2 -2  4
        ])
    )

    z0 = [-10.0, 10.0, -5.0, 5.0]

    status, z, info = PATHSolver.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        z0;
        variable_names = ["x1", "x2", "x3", "x4"],
        constraint_names = ["F1", "F2", "F3", "F4"],
        output = "yes",
    )
    @test status == PATHSolver.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])

    status, z, info = PATHSolver.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        z0;
        variable_names = ["x1", "x2", "x3", "x4"],
        constraint_names = ["A2345678901234567890", "B23456789012345678901234567", "C234567890", "D234567890"],
        output = "yes",
    )
    @test status == PATHSolver.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])

    status, z, info = PATHSolver.solve_mcp(
        M,
        Float64[2, 2, -2, -6],
        fill(0.0, 4),
        fill(10.0, 4),
        z0;
        variable_names = ["x1", "x2", "x3", "x4"],
        constraint_names = ["A2345678901234567890", "B23456789012345678901234567", "C🤖λ⚽⚽⚽⚽🏸⚽⚽⚽⚽🏸", "D도레미파솔"],
        output = "yes",
    )
    @test status == PATHSolver.MCP_Solved
    @test isapprox(z, [2.8, 0.0, 0.8, 1.2])
end


text:
@testset "$(_gen("Optimizer: Gamma Loss [Small][Forgiving][Native]"))" begin
    @timeit to "Optimizer: Gamma Loss [Small][Forgiving][Native]" begin
        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 1000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 3000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 5000,
                                      ρ                  = 0,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 1000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 3000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,max_iter = 100) > 0.9
        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 5000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9

        @test unit_test_train_subloss(AbstractGamma(),
                                      gradient_eval      = Loss{AbstractGamma}(),
                                      input_distribution = Distributions.Gamma(5,0.5),
                                      input_size         = 5000,
                                      ρ                  = 0.2,
                                      step_size          = 0.1,
                                      max_iter           = 100) > 0.9
    end
end


text:
using RINEX
using Test

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


text:
module FovWrapper

using JavaCall

const AbstractFieldOfView = @jimport org.orekit.geometry.fov.AbstractFieldOfView
const CircularFieldOfView = @jimport org.orekit.geometry.fov.CircularFieldOfView
const Class = @jimport java.lang.Class
const DoubleDihedraFieldOfView = @jimport org.orekit.geometry.fov.DoubleDihedraFieldOfView
const EllipticalFieldOfView = @jimport org.orekit.geometry.fov.EllipticalFieldOfView
const FieldOfView = @jimport org.orekit.geometry.fov.FieldOfView
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const Object = @jimport java.lang.Object
const OneAxisEllipsoid = @jimport org.orekit.bodies.OneAxisEllipsoid
const PolygonalFieldOfView = @jimport org.orekit.geometry.fov.PolygonalFieldOfView
const PolygonalFieldOfView_DefiningConeType = @jimport org.orekit.geometry.fov.PolygonalFieldOfView$DefiningConeType
const SmoothFieldOfView = @jimport org.orekit.geometry.fov.SmoothFieldOfView
const SphericalPolygonsSet = @jimport org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet
const Transform = @jimport org.orekit.frames.Transform
const Vector3D = @jimport org.hipparchus.geometry.euclidean.threed.Vector3D
const VisibilityTrigger = @jimport org.orekit.propagation.events.VisibilityTrigger

export AbstractFieldOfView
export CircularFieldOfView
export DoubleDihedraFieldOfView
export EllipticalFieldOfView
export FieldOfView
export PolygonalFieldOfView
export SmoothFieldOfView
export equals
export get_center
export get_class
export get_focus1
export get_focus2
export get_footprint
export get_half_aperture
export get_half_aperture_along_x
export get_half_aperture_along_y
export get_margin
export get_x
export get_y
export get_z
export get_zone
export hash_code
export notify
export notify_all
export offset_from_boundary
export project_to_boundary
export to_string
export wait

include("abstract_field_of_view.jl")
include("circular_field_of_view.jl")
include("double_dihedra_field_of_view.jl")
include("elliptical_field_of_view.jl")
include("field_of_view.jl")
include("polygonal_field_of_view.jl")
include("smooth_field_of_view.jl")

end


text:
using Sockets
# @async begin
# Base.

fileName="D:\\fiber\\c\\Data^Ver00^xxxx-001^20190212131313^ID00.RAW3"
Port=19998
IPAddr="127.0.0.1"
IPAddr="192.168.0.109"
glbMachineName="xxxx-001"
include("socketinter.jl")
include("SvrTsk.jl")
const glbtypecode= 0x0003
println("File reader  begin :")
a= @async  SVrReadFileTsk(fileName)

println("Socket Server begin :")
b = SVrSendSocketTskEntry(Port, IPAddr)
println(a)

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
module ClimateModels

using Zarr, AWSCore, DataFrames, CSV, CFTime, Dates, Statistics

export cmip

"""
    cmip(institution_id,source_id,variable_id)

Access CMIP6 climate model archive (https://bit.ly/2WiWmoh) via
`AWSCore.jl` and `Zarr.jl` and compute (1) time mean global map and
(2) time evolving global mean.

This was partly inspired by @rabernat 's https://bit.ly/2VRMgvl notebook

```
using ClimateModels
(mm,gm,meta)=cmip()
nm=meta["long_name"]*" in "*meta["units"]

using Plots
heatmap(mm["lon"], mm["lat"], transpose(mm["m"]),
        title=nm*" (time mean)")
plot(gm["t"][1:12:end],gm["y"][1:12:end],xlabel="time",ylabel=nm,
     title=meta["institution_id"]*" (global mean, month by month)")
display.([plot!(gm["t"][i:12:end],gm["y"][i:12:end], leg = false) for i in 2:12])

```
"""
function cmip(institution_id="IPSL",source_id="IPSL-CM6A-LR",
    variable_id="tas")

    #choose model and variable
    S=[institution_id, source_id, variable_id]

    #initiate cloud connection
    ⅁ = AWSCore.aws_config(creds=nothing, region="",
    service_host="googleapis.com", service_name="storage")

    #get list of contents for cloud storage unit
    β = S3Store("cmip6","", aws=⅁, listversion=1)
    ξ = CSV.read(IOBuffer(β["cmip6-zarr-consolidated-stores.csv"]))

    # get model grid cell areas
    ii=findall( (ξ[!,:source_id].==S[2]).&(ξ[!,:variable_id].=="areacella") )
    μ=ξ[ii,:]
    i1=findfirst("cmip6",μ.zstore[end])[end]+2
    P = μ.zstore[end][i1:end]
    ζ = zopen(S3Store("cmip6", P, aws=⅁, listversion=1))
    Å = ζ["areacella"][:, :];

    # get model solution ensemble list
    i=findall( (ξ[!,:activity_id].=="CMIP").&(ξ[!,:table_id].=="Amon").&
    (ξ[!,:variable_id].==S[3]).&(ξ[!,:experiment_id].=="historical").&
    (ξ[!,:institution_id].==S[1]) )
    μ=ξ[i,:]

    # access one model ensemble member
    i1=findfirst("cmip6",μ.zstore[end])[end]+2
    P = μ.zstore[end][i1:end]
    ζ = zopen(S3Store("cmip6", P, aws=⅁, listversion=1))

    meta=Dict("institution_id" => institution_id,"source_id" => source_id,
        "variable_id" => variable_id, "units" => ζ[S[3]].attrs["units"],
        "long_name" => ζ[S[3]].attrs["long_name"])

    # time mean global map
    m = convert(Array{Union{Missing, Float32},3},ζ[S[3]][:,:,:])
    m = dropdims(mean(m,dims=3),dims=3)

    mm=Dict("lon" => ζ["lon"], "lat" => ζ["lat"], "m" => m)

    # time evolving global mean
    t = ζ["time"]
    t = timedecode(t[:], t.attrs["units"], t.attrs["calendar"])

    y = ζ[S[3]][:,:,:]
    y=[sum(y[:, :, i].*Å) for i in 1:length(t)]./sum(Å)

    gm=Dict("t" => t, "y" => y)

    return mm,gm,meta
end

end # module


text:
export ElectronVolt, Torr

const ElectronVolt = NonSIUnit{typeof(Joule),:eV}()
convert(::Type{SIQuantity},::typeof(ElectronVolt)) = 1.60217656535e-19Joule

const Torr = NonSIUnit{typeof(Pascal),:torr}()
convert(::Type{SIQuantity},::typeof(Torr)) = 133.322368Pascal


text: horner2(coef,x) = foldr((u,v) -> u + x*v, 0, coef)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
path = dirname(@__FILE__)

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:cond_id, 0)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:date_conditional_end, quartertodate("2015-Q4"))
m <= Setting(:use_population_forecast, true)

estroot = normpath(joinpath(dirname(@__FILE__), "..", "reference"))
overrides = forecast_input_file_overrides(m)
overrides[:mode] = joinpath(estroot, "optimize.h5")
overrides[:full] = joinpath(estroot, "metropolis_hastings.h5")

skip_forecast_one_draw = false
if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(ENV["HOME"],".freddatarc"))
    df = load_data(m)
    skip_forecast_one_draw = true
else
    @warn "Skipping forecast_one_draw tests because FRED_API_KEY not present"
end
# Make sure output_vars ignores the untransformed and 4Q things because they are
# computed in compute_meansbands
output_vars = add_requisite_output_vars([:histpseudo, :histobs, :histstdshocks,
                                         :histutpseudo, :histutobs,
                                         :hist4qpseudo, :hist4qobs,
                                         :forecaststates, :forecastpseudo, :forecastobs, :forecaststdshocks,
                                         :forecastutpseudo, :forecastutobs,
                                         :forecast4qpseudo, :forecast4qobs,
                                         :bddforecaststates, :bddforecastshocks, :bddforecastpseudo, :bddforecastobs,
                                         :shockdecpseudo, :shockdecobs,
                                         :trendstates, :trendobs, :trendpseudo,
                                         :dettrendstates, :dettrendobs, :dettrendpseudo,
                                         :irfstates, :irfpseudo, :irfobs])

# Check error handling for input_type = :subset
@testset "Ensure properly error handling for input_type = :subset" begin
    @test_throws ErrorException forecast_one(m, :subset, :none, output_vars,
                                    subset_inds = 1:10, forecast_string = "",
                                    verbose = :none)
    @test_throws ErrorException forecast_one(m, :subset, :none, output_vars,
                                    forecast_string = "test",
                                    verbose = :none)
end

# Run modal forecasts
out = Dict{Symbol, Dict{Symbol, Array{Float64}}}()
for cond_type in [:none, :semi, :full]
    forecast_one(m, :mode, cond_type, output_vars, verbose = :none)
    forecast_one(m, :init, cond_type, output_vars, verbose = :none)

    # Read output
    out[cond_type] = Dict{Symbol, Array{Float64}}()
    output_files = get_forecast_output_files(m, :mode, cond_type, output_vars)
    for var in keys(output_files)
        out[cond_type][var] = load(output_files[var], "arr")
    end
end

# Read expected output
exp_out = JLD2.jldopen("$path/../reference/forecast_one_out.jld2", "r") do file
    read(file, "exp_out")
end

# Test modal forecast outputs
specify_mode!(m, DSGE.get_forecast_input_file(m, :mode); verbose = :none)

@testset "Test modal forecast for all major output_vars" begin
    for cond_type in [:none, :semi, :full]
        # Histories
        @test @test_matrix_approx_eq exp_out[cond_type][:histpseudo]     out[cond_type][:histpseudo]

        # Forecasts
        @test @test_matrix_approx_eq exp_out[cond_type][:forecastobs]    out[cond_type][:forecastobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:forecastpseudo] out[cond_type][:forecastpseudo]

        # Shock decompositions, deterministic trends, trends
        @test @test_matrix_approx_eq exp_out[cond_type][:shockdecobs]    out[cond_type][:shockdecobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:shockdecpseudo] out[cond_type][:shockdecpseudo]
        @test @test_matrix_approx_eq exp_out[cond_type][:dettrendobs]    out[cond_type][:dettrendobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:dettrendpseudo] out[cond_type][:dettrendpseudo]
        @test @test_matrix_approx_eq exp_out[cond_type][:trendobs]       out[cond_type][:trendobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:trendpseudo]    out[cond_type][:trendpseudo]

        # IRFs
        @test @test_matrix_approx_eq exp_out[cond_type][:irfobs]         out[cond_type][:irfobs]
        @test @test_matrix_approx_eq exp_out[cond_type][:irfpseudo]      out[cond_type][:irfpseudo]
    end
end

@testset "Test full-distribution forecasts run" begin
    m <= Setting(:forecast_block_size, 5)
    forecast_one(m, :prior, :none, output_vars, verbose = :none)
    for sampling_method in [:MH]
        m <= Setting(:sampling_method, sampling_method)
        for cond_type in [:none, :semi, :full]
            forecast_one(m, :full, cond_type, output_vars, verbose = :none)
            @test_throws ErrorException forecast_one(m, :subset, cond_type, output_vars, subset_inds = 1:10, verbose = :none)
            forecast_one(m, :subset, cond_type, output_vars, subset_inds = 1:10, forecast_string = "test", verbose = :none)
            forecast_one(m, :init_draw_shocks, cond_type, output_vars, verbose = :none)
            forecast_one(m, :mode_draw_shocks, cond_type, output_vars, verbose = :none)
        end
    end
end

# Test full-distribution blocking
@everywhere using DSGE
m <= Setting(:forecast_block_size, 5)
forecast_one(m, :full, :none, output_vars, verbose = :none)
# Test read_forecast_output
@testset "Test full-distribution blocking" begin
    for input_type in [:mode, :full]
        output_files = get_forecast_output_files(m, input_type, :none, output_vars)
        @test ndims(DSGE.read_forecast_series(output_files[:trendobs], :trend, m.observables[:obs_gdp])) == 2
        @test ndims(DSGE.read_forecast_series(output_files[:forecastobs], :forecast, m.observables[:obs_gdp])) == 2
        @test ndims(DSGE.read_forecast_series(output_files[:irfobs], m.observables[:obs_gdp], m.exogenous_shocks[:rm_sh])) == 2
    end
end

@testset "Test forecast_one_draw" begin
    if !skip_forecast_one_draw
        for input_type in [:mode, :full]
            params = if input_type == :mode
                load_draws(m, input_type)
            else
                load_draws(m, input_type)[1, :]
            end
            m <= Setting(:alternative_policy, AltPolicy(:historical, eqcond, solve))
            @test typeof(DSGE.forecast_one_draw(m, input_type, :none, output_vars, params, df)) == Dict{Symbol, Array{Float64}}

            # Test with alternative policy
            m <= Setting(:alternative_policy, AltPolicy(:taylor93, eqcond, solve))
            @test typeof(DSGE.forecast_one_draw(m, input_type, :none, output_vars, params, df)) == Dict{Symbol, Array{Float64}}
        end
    end
end

nothing


text:
"""
    parse_config_file(filename::AbstractString)::Dict where {T <: SSDFloat}

Reads in a configuration file and returns a parsed dictionary which holds all the
information specified in the configuration file.

Find detailed information on configuration files in [Configuration Files](@ref).

## Arguments
* `filename::AbstractString`: File name of the configuration file. If the file is not
    in the same directory, a path to the file is required.

!!! note
    Currently supported formats for the configuration files:
        - [YAML](https://github.com/JuliaData/YAML.jl): `filename` ends with `.yaml`.
        - [JSON](https://github.com/JuliaIO/JSON.jl): `filename` ends with `.json`.
        - [SigGen](https://github.com/radforddc/icpc_siggen): `filename` ends with `.config`.
"""
function parse_config_file(filename::AbstractString)::Dict{Any,Any}
    if endswith(filename, ".toml")
        error("Currently only .json and .yaml files are supported.")
    elseif endswith(filename, ".json")
        dicttext = read(filename, String)
        dict = JSON.parse(dicttext)
        scan_and_merge_included_json_files!(dict, filename)
    elseif endswith(filename, ".yaml")
        dict = YAML.load_file(filename)
        scan_and_merge_included_json_files!(dict, filename)
    elseif endswith(filename, ".config")
        siggen_dict = readsiggen(filename)
        dict = siggentodict(siggen_dict)
    else
        error("Currently only .json and .yaml files are supported.")
    end
    dict
end

function yaml2json_convert(filename::String)
    data = YAML.load(open(filename), dicttype = DataStructures.OrderedDict)
    if isempty(data)
        @warn "The file $(filename) is empty and will not be converted."
    else
        open(replace(filename, ".yaml" => ".json"),"w") do f
            JSON.print(f, data, 4)
        end
    end
end
function yaml2json(directory::String)# or filename
    if endswith(directory, ".yaml")
        yaml2json_convert(directory)
    else
        yamlfiles = filter(x->endswith(x,".yaml"), joinpath.(directory, readdir(directory)))
        for filename in yamlfiles
            yaml2json_convert(filename)
        end
    end
end

function json2yaml_convert(filename::String)::Nothing
    data = JSON.parsefile(filename, dicttype = DataStructures.OrderedDict)
    if isempty(data)
        @warn "The file $(filename) is empty and will not be converted."
    else
        YAML.write_file(replace(filename, ".json" => ".yaml"), data)
    end

end
function json2yaml(directory::String)# or filename
    if endswith(directory, ".json")
        json2yaml_convert(directory)
    else
        jsonfiles = filter(x->endswith(x,".json"), joinpath.(directory, readdir(directory)))
        for filename in jsonfiles
            json2yaml_convert(filename)
        end
    end
end

function scan_and_merge_included_json_files!(dict, config_filename::AbstractString)
    key_word = "include"
    config_dir = dirname(config_filename)
    for k in keys(dict)
        is_subdict = typeof(dict[k]) <: Dict
        if !is_subdict && string(k) != key_word
            typeof(dict[k]) <: Array ? is_subdict = true : is_subdict = false
        end
        if is_subdict
            scan_and_merge_included_json_files!(dict[k], config_filename)
        elseif string(k) == key_word
            files = []
            if typeof(dict[k]) <: Array
                append!(files, dict[k])
            else
                push!(files, dict[k])
            end
            for file in files
                filepath = joinpath(config_dir, file)
                if isfile(filepath)
                    tmp = parse_config_file(filepath)
                    for sub_k in keys(tmp)
                        dict[sub_k] = tmp[sub_k]
                    end
                end
            end
            delete!(dict, k)
        end
    end
end


text:
using IMDb
using Test

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


text:
# Workaround for JuliaLang/julia/pull/28625
if Base.HOME_PROJECT[] !== nothing
  Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

using
  Documenter,
  FourierFlows

makedocs(
   modules = [FourierFlows],
   doctest = false,
     clean = true,
 checkdocs = :all,
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
   authors = "Gregory L. Wagner and Navid C. Constantinou",
  sitename = "FourierFlows.jl",
     pages = Any[
              "Home" => "index.md",
              "Code Basics" => "basics.md",
              "Forcing" => "forcing.md",
              "DocStrings" => Any[
              "man/types.md",
              "man/functions.md"]
             ]
)

deploydocs(repo = "github.com/FourierFlows/FourierFlows.jl.git")


text:
using LinearAlgebra

a = [1 2;
     3 4]

b = [5 6 7;
     8 9 10;
     10 9 8]

c = [11 12 13 14;
     15 16 17 18;
     19 20 21 22;
     23 24 25 26]

println(a)
println(b)
println(c)


text:
abstract type EDTask end

# output arrays
function initialize! end

# can preallocate compute space
# Note: This needs to wrap the Task in an additional struct
# The additional memory may not be a part of the same struct
# otherwise the code breaks in a multithreaded (not multiprocessed) setting
# as all threads share the same memory location for their temporary results
initialize_local(task::EDTask) = task

# task, ρindex, shot, fieldindex, eigen
function compute_task! end

# task, ρindex, shot, fieldindex
function failed_task! end

# task, EDDataDescriptor -> Data object
function assemble end


text:
export DynaAgent

"""
    DynaAgent(;kwargs...)

`DynaAgent` is first introduced in: *Sutton, Richard S. "Dyna, an integrated architecture for learning, planning, and reacting." ACM Sigart Bulletin 2.4 (1991): 160-163.*

# Keywords & Fields

- `policy`::[`AbstractPolicy`](@ref): the policy to use
- `model`::[`AbstractEnvironmentModel`](@ref): describe the environment to interact with
- `trajectory`::[`AbstractTrajectory`](@ref): used to store transitions between agent and environment
- `role=:DEFAULT`: used to distinguish different agents
- `plan_step::Int=10`: the count of planning steps

The main difference between [`DynaAgent`](@ref) and [`Agent`](@ref) is that an environment model is involved. It is best described in the book: *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*

![](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/raw/master/docs/src/assets/img/RL_book_fig_8_1.png)
![](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/raw/master/docs/src/assets/img/RL_book_fig_8_2.png)
"""
Base.@kwdef struct DynaAgent{
    P<:AbstractPolicy,
    B<:AbstractTrajectory,
    M<:AbstractEnvironmentModel,
    R,
} <: AbstractAgent
    policy::P
    model::M
    trajectory::B
    role::R = :DEFAULT_PLAYER
    plan_step::Int = 10
end

get_role(agent::DynaAgent) = agent.role

function (agent::DynaAgent{<:AbstractPolicy,<:EpisodicTrajectory})(::PreEpisodeStage, env)
    empty!(agent.trajectory)
    nothing
end

function (agent::DynaAgent{<:AbstractPolicy,<:EpisodicTrajectory})(::PreActStage, env)
    action = agent.policy(env)
    push!(agent.trajectory; state = get_state(env), action = action)
    update!(agent.model, agent.trajectory, agent.policy)  # model learning
    update!(agent.policy, agent.trajectory)  # direct learning
    update!(agent.policy, agent.model, agent.trajectory, agent.plan_step)  # policy learning
    action
end

function (agent::DynaAgent{<:AbstractPolicy,<:EpisodicTrajectory})(::PostActStage, env)
    push!(agent.trajectory; reward = get_reward(env), terminal = get_terminal(env))
    nothing
end

function (agent::DynaAgent{<:AbstractPolicy,<:EpisodicTrajectory})(::PostEpisodeStage, env)
    action = agent.policy(env)
    push!(agent.trajectory; state = get_state(env), action = action)
    update!(agent.model, agent.trajectory, agent.policy)  # model learning
    update!(agent.policy, agent.trajectory)  # direct learning
    update!(agent.policy, agent.model, agent.trajectory, agent.plan_step)  # policy learning
    action
end

"By default, only use trajectory to update model"
RLBase.update!(model::AbstractEnvironmentModel, t::AbstractTrajectory, π::AbstractPolicy) =
    update!(model, t)

# function RLBase.update!(model::AbstractEnvironmentModel, buffer::AbstractTrajectory)
#     transitions = extract_experience(buffer, model)
#     isnothing(transitions) || update!(model, transitions)
# end


text:
module DiffOperatorsTests

using Test
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Fields: MockField
using LinearAlgebra
using FillArrays

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

v = 3.0
f = MockField(v)

@test ∇(f) == gradient(f)

@test divergence(f) == tr(gradient(f))

@test curl(f) == grad2curl(gradient(f))

@test ∇⋅f == divergence(f)

@test cross(∇,f) == curl(f)

@test ∇×f == curl(f)

@test outer(∇,f) == ∇(f)

@test ∇⊗f == ∇(f)

@test outer(f,∇) == transpose(∇(f))

@test f⊗∇ == transpose(∇(f))

@test ε(f) == symmetric_part(gradient(f))

@test Δ(f) == ∇⋅∇(f)

l = 10
f = Fill(f,l)

@test Broadcasting(divergence)(f) == Broadcasting(Operation(tr))(Broadcasting(∇)(f))
@test Broadcasting(curl)(f) == Broadcasting(Operation(grad2curl))(Broadcasting(∇)(f))
@test Broadcasting(ε)(f) == Broadcasting(Operation(symmetric_part))(Broadcasting(∇)(f))

# Test automatic differentiation

u_scal(x) = x[1]^2 + x[2]
∇u_scal(x) = VectorValue( 2*x[1], one(x[2]) )
Δu_scal(x) = 2

u_vec(x) = VectorValue( x[1]^2 + x[2], 4*x[1] - x[2]^2 )
∇u_vec(x) = TensorValue( 2*x[1], one(x[2]), 4*one(x[1]), - 2*x[2] )
Δu_vec(x) = VectorValue( 2, -2 )
εu_vec(x) = SymTensorValue( 2*x[1], 0.5*(one(x[2])+4*one(x[1])), - 2*x[2] )

xs = [ Point(1.,1.), Point(2.,0.), Point(0.,3.), Point(-1.,3.)]
for x in xs
  @test ∇(u_scal)(x) == ∇u_scal(x)
  @test Δ(u_scal)(x) == Δu_scal(x)
  @test ∇(u_vec)(x) == ∇u_vec(x)
  @test Δ(u_vec)(x) == Δu_vec(x)
  @test ε(u_vec)(x) == εu_vec(x)
end

u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
∇u(x) = TensorValue( 2*x[1], 4*x[2], -2*x[1], zero(x[1]) )
Δu(x) = VectorValue( 6, -2 )
εu(x) = SymTensorValue( 2*x[1], 2*x[2]-x[1], zero(x[1]) )

for x in xs
  @test (∇⋅u)(x) == tr(∇u(x))
  @test (∇×u)(x) == grad2curl(∇u(x))
  @test Δ(u)(x) == Δu(x)
  @test ε(u)(x) == εu(x)
end

u(x) = VectorValue( x[1]^2 + 2*x[2]^2, 0 )
∇u(x) = TensorValue( 2*x[1], 4*x[2], 0, 0 )
Δu(x) = VectorValue( 6, 0 )
εu(x) = SymTensorValue( 2*x[1], 2*x[2], 0 )

for x in xs
  @test (∇⋅u)(x) == tr(∇u(x))
  @test (∇×u)(x) == grad2curl(∇u(x))
  @test Δ(u)(x) == Δu(x)
  @test ε(u)(x) == εu(x)
end

end # module


text:
include("../src/experiments/0_setup.jl")
import DataFrames.DataFrame
using DataFrames, CSV
using FastGroupBy
using Base.Threads
K = 100

tries = vcat([Int(2^k-1) for k = 7:31], 3_000_000_000)
for N in tries
    println(N)
    if N < 2_000_000
        by = nothing; val = nothing; gc()
        srand(1)
        by = rand(Int64(1):Int64(round(N/K)), N);
        val = rand(Int32(1):Int32(5), N);
        sp  = @elapsed sumby_sortperm(by, val)
        CSV.write(string("benchmark/out/64/sp$N $(replace(string(now()),":","")).csv"),DataFrame(sp = sp))
    end

    by = nothing; val = nothing; gc()
    srand(1)
    by = rand(Int64(1):Int64(round(N/K)), N);
    val = rand(Int32(1):Int32(5), N);
    srs = @elapsed sumby_multi_rs(by, val)
    CSV.write(string("benchmark/out/64/mrs$N $(replace(string(now()),":","")).csv"),DataFrame(srs = srs))

    if N < 2_000_000_000
        by = nothing; val = nothing; gc()
        srand(1)
        by = rand(Int64(1):Int64(round(N/K)), N);
        val = rand(Int32(1):Int32(5), N);
        srg = @elapsed sumby_radixgroup(by, val)
        CSV.write(string("benchmark/out/64/srg$N $(replace(string(now()),":","")).csv"),DataFrame(srg = srg))
    end

    by = nothing; val = nothing; gc()
    srand(1)
    by = rand(Int64(1):Int64(round(N/K)), N);
    val = rand(Int32(1):Int32(5), N);
    srs = @elapsed sumby_radixsort(by, val)
    CSV.write(string("benchmark/out/64/srs$N $(replace(string(now()),":","")).csv"),DataFrame(srs = srs))
end

tries = vcat([Int(2^k-1) for k = 7:31], 3_000_000_000)
for N in tries
    println(N)
    if N < 400_000
        by = nothing; val = nothing; gc()
        srand(1)
        by = rand(Int32(1):Int32(round(N/K)), N);
        val = rand(Int32(1):Int32(5), N);
        sp  = @elapsed sumby_sortperm(by, val)
        CSV.write(string("benchmark/out/sp$N $(replace(string(now()),":","")).csv"),DataFrame(sp = sp))
    end

    by = nothing; val = nothing; gc()
    srand(1)
    by = rand(Int32(1):Int32(round(N/K)), N);
    val = rand(Int32(1):Int32(5), N);
    srs = @elapsed sumby_multi_rs(by, val)
    CSV.write(string("benchmark/out/mrs$N $(replace(string(now()),":","")).csv"),DataFrame(srs = srs))

    if N < 2_000_000_000
        by = nothing; val = nothing; gc()
        srand(1)
        by = rand(Int32(1):Int32(round(N/K)), N);
        val = rand(Int32(1):Int32(5), N);
        srg = @elapsed sumby_radixgroup(by, val)
        CSV.write(string("benchmark/out/srg$N $(replace(string(now()),":","")).csv"),DataFrame(srg = srg))
    end

    by = nothing; val = nothing; gc()
    srand(1)
    by = rand(Int32(1):Int32(round(N/K)), N);
    val = rand(Int32(1):Int32(5), N);
    srs = @elapsed sumby_radixsort(by, val)
    CSV.write(string("benchmark/out/srs$N $(replace(string(now()),":","")).csv"),DataFrame(srs = srs))
end


text:
using Base: Test, LinAlg
using HDF5

using DSGE
include("../util.jl")
path = dirname(@__FILE__)

m = Model990()
Γ0, Γ1, C, Ψ, Π = eqcond(m)
stake = 1 + 1e-6
G1, C, impact, fmat, fwt, ywt, gev, eu, loose = gensys(Γ0, Γ1, C, Ψ, Π, stake)

h5 = h5open("$path/../reference/gensys.h5")
G1_exp = read(h5, "G1_gensys")
C_exp = reshape(read(h5, "C_gensys"), 60, 1)
impact_exp = read(h5, "impact")
eu_exp = read(h5, "eu")
close(h5)

@test_matrix_approx_eq G1_exp G1

@test_matrix_approx_eq C_exp C

@test_matrix_approx_eq impact_exp impact

@test isequal(eu_exp, eu)


text:
"""
        electre(decisionMat, weights, fns)

Apply ELECTRE (ELimination Et Choice Translating REality) method
for a given matrix and weights.

# Arguments:
 - `decisionMat::DataFrame`: n × m matrix of objective values for n candidate (or strategy) and m criteria
 - `weights::Array{Float64, 1}`: m-vector of weights that sum up to 1.0. If the sum of weights is not 1.0, it is automatically normalized.
 - `fns::Array{Function, 1}`: m-vector of function that are either minimize or maximize.

# Description
electre() applies the ELECTRE method to rank n strategies subject to m criteria which are supposed to be either maximized or minimized.
C and D values are used to determine the best strategy. If the strategy with the highest C value
is same as the strategy with the lowest D value than the solution is unique. Otherwise, two strategies
are reported as the solution.

# Output
- `::ElectreResult`: ElectreResult object that holds multiple outputs including scores and best index.

# Examples
```julia-repl
julia> w =  [0.110, 0.035, 0.379, 0.384, 0.002, 0.002, 0.010, 0.077];
julia> Amat = [
      100 92 10 2 80 70 95 80 ;
      80  70 8  4 100 80 80 90 ;
      90 85 5 0 75 95 70 70 ;
      70 88 20 18 60 90 95 85
    ];

julia> dmat = makeDecisionMatrix(Amat);
julia> fns = makeminmax([maximum for i in 1:8]);
julia> result = electre(dmat, w, fns)

julia> dmat
4×8 DataFrame
 Row │ Crt1     Crt2     Crt3     Crt4     Crt5     Crt6     Crt7     Crt8
     │ Float64  Float64  Float64  Float64  Float64  Float64  Float64  Float64
─────┼────────────────────────────────────────────────────────────────────────
   1 │   100.0     92.0     10.0      2.0     80.0     70.0     95.0     80.0
   2 │    80.0     70.0      8.0      4.0    100.0     80.0     80.0     90.0
   3 │    90.0     85.0      5.0      0.0     75.0     95.0     70.0     70.0
   4 │    70.0     88.0     20.0     18.0     60.0     90.0     95.0     85.0

julia> result.bestIndex
(4,)

julia> result.C
4-element Array{Float64,1}:
  0.3693693693693696
  0.01501501501501501
 -2.473473473473473
  2.0890890890890894

julia> result.D
4-element Array{Float64,1}:
  0.1914244325928971
 -0.19039293350192432
  2.884307608766315
 -2.885339107857288
```

# References
Celikbilek Yakup, Cok Kriterli Karar Verme Yontemleri, Aciklamali ve Karsilastirmali
Saglik Bilimleri Uygulamalari ile. Editor: Muhlis Ozdemir, Nobel Kitabevi, Ankara, 2018
"""
function electre(decisionMat::DataFrame, weights::Array{Float64,1}, fns::Array{Function,1})::ElectreResult

    w = unitize(weights)

    nalternatives, ncriteria = size(decisionMat)

    normalizedMat = normalize(decisionMat)
    weightednormalizedMat = w * normalizedMat

    fitnessTable = []
    nonfitnessTable = []

    for i in 1:nalternatives
        for j in 1:nalternatives
            if i != j
                betterlist = []
                worselist = []
                for h in 1:ncriteria
                    if fns[h] == maximum
                        if weightednormalizedMat[i, h] >= weightednormalizedMat[j, h]
                            push!(betterlist, h)
                        else
                            push!(worselist, h)
                        end
                    else
                        if weightednormalizedMat[i, h] <= weightednormalizedMat[j, h]
                push!(betterlist, h)
                        else
                            push!(worselist, h)
                        end
                    end
                end
                fitnesstableelement = Dict(:i => i, :j => j, :set => betterlist)
                push!(fitnessTable, fitnesstableelement)

                nonfitnesstableelement = Dict(:i => i, :j => j, :set => worselist)
                push!(nonfitnessTable, nonfitnesstableelement)
            end
        end
    end

    fitnessmatrix = zeros(Float64, nalternatives, nalternatives)
    nonfitnessmatrix = zeros(Float64, nalternatives, nalternatives)

    for elements in fitnessTable
        i = elements[:i]
        j = elements[:j]
        elemlist = elements[:set]

        CC = sum(w[elemlist])
        fitnessmatrix[i, j] = CC
    end

    nonfitnessmatrix = zeros(Float64, nalternatives, nalternatives)
    for elements in nonfitnessTable
        i = elements[:i]
        j = elements[:j]
        elemlist = elements[:set]

        r_ik       = weightednormalizedMat[i, elemlist]
        r_jk       = weightednormalizedMat[j, elemlist]
        r_ik_full  = weightednormalizedMat[i, :]
        r_jk_full  = weightednormalizedMat[j, :]

        if length(r_ik) > 0 && length(r_jk) > 0
            nom = maximum(abs.(r_ik - r_jk))
            dom = maximum(abs.(r_ik_full - r_jk_full))
            nonfitnessmatrix[i, j] = nom / dom
        else
            nonfitnessmatrix[i, j] = Inf64
        end  # end if

    end # end for

    C = zeros(Float64, nalternatives)
    D = zeros(Float64, nalternatives)
    for i in 1:nalternatives
        C[i] = sum(fitnessmatrix[i, :]) - sum(fitnessmatrix[:, i])
        D[i] = sum(nonfitnessmatrix[i, :]) - sum(nonfitnessmatrix[:, i])
    end

    best_C_index = sortperm(C) |> last
    best_D_index = sortperm(D) |> first
    best = nothing

    if best_C_index == best_D_index
        best = (best_C_index,)
    else
        best = (best_C_index, best_D_index)
    end

    result = ElectreResult(
        decisionMat,
        w,
        weightednormalizedMat,
        fitnessTable,
        nonfitnessTable,
        fitnessmatrix,
        nonfitnessmatrix,
        C,
        D,
        best
    )

    return result
end

"""
        electre(setting)

Apply ELECTRE (ELimination Et Choice Translating REality) method
for a given matrix and weights.

# Arguments:
 - `setting::MCDMSetting`: MCDMSetting object.

# Description
electre() applies the ELECTRE method to rank n strategies subject to m criteria which are supposed to be either maximized or minimized.
C and D values are used to determine the best strategy. If the strategy with the highest C value
is same as the strategy with the lowest D value than the solution is unique. Otherwise, two strategies
are reported as the solution.

# Output
- `::ElectreResult`: ElectreResult object that holds multiple outputs including scores and best index.
"""
function electre(setting::MCDMSetting)::ElectreResult
    electre(
        setting.df,
        setting.weights,
        setting.fns
    )
end

max_stars_repo_name: jbytecode/jMCDM

max_issues_repo_name: jbytecode/jMCDM

max_forks_repo_name: jbytecode/jMCDM


text:
# Plot attributes
theme(:mute)
default(; frame=:box)
Plots.scalefontsizes(1.3)
fonts=font("Arial",7)

### Heatmap of disconnected species ###

# Get the MaxEnt degree distribution for a range of species richness and different quantiles of the number of links
Ss = [5:5:100;] # range of species richness
qs = [0:0.05:1;] # quantiles of the number of links
p_disconnected = zeros(length(qs), length(Ss)) # create vector for the probability a species is disconnected given the number of links (rows) and species richness (columns)

for (i, S) in enumerate(Ss)
  Ls = round.(quantile(0:S^2, qs)) # get the number of links L associated to every quantile (proportion of S^2 links below L)
  Ls = convert(Vector{Int64}, Ls)

  degree_dists = degree_dist_maxent.(S, Ls) # get the degree distributions of maximum entropy for a given S and all L
  degree_dists = reduce(hcat, degree_dists)
  p_disconnected[:,i] = degree_dists[1,:]
end

q_minS = Ss ./ (Ss.^2 .+ 1) # get the quantile associated to the minimum number of links S-1 that would allow all species to be connected

# plot the proportion of disconnected species for all combinations of S and L (L in quantile)
heatmap(Ss, qs, log.(p_disconnected),
        c=:viridis,
        colorbar_title="log probability",
        framestyle=:box,
        dpi=1000,
        size=(800,500),
        margin=5Plots.mm,
        guidefont=fonts,
        xtickfont=fonts,
        ytickfont=fonts)
# add the minimum number of links S-1 (quantile) on the plot
scatter!(Ss, q_minS,
      linewidth=1,
      linecolor=:black,
      markerwidth=3,
      markercolor=:black,
      lab="")
plot!(Ss, q_minS,
      linewidth=1,
      linecolor=:black,
      lab="")
xaxis!(xlim=(minimum(Ss)-2.5, maximum(Ss)+2.5), "species richness")
yaxis!(ylim=(0,1), "quantile of the number of links")

savefig(joinpath("figures", "prop_disconnected_species.png"))

################# TK TO DO #########################
## Figure: Mean degree constraints and MaxEnt degree distributions
species015 = Int64(round(species015))
species500 = Int64(round(species500))
species985 = Int64(round(species985))

"""
get_avgk_dist(S::Int64)
    S: number of species
Returns the distribution of mean degrees according to the predictions of the flexible links model
"""
function get_avgk_dist(S::Int64)
  L = predicted_links[:, S-4]
  kavg = 2 .* L ./ S
  return kavg
end

# Get distribution of mean degrees from species richness
kavg015 = get_avgk_dist(species015)
kavg500 = get_avgk_dist(species500)
kavg985 = get_avgk_dist(species985)

# Plot distributions of mean degrees of 3 different levels of species richness
plotA = density(kavg015, linewidth=2, label="$(species015) species",
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white, legendfont=fonts)
density!(kavg500, linewidth=2, label="$(species500) species")
density!(kavg985, linewidth=2, label="$(species985) species")
xaxis!(xlabel="Mean degree (links/species)", xlims=(0,maximum(kavg985)))
yaxis!(ylabel="Density", ylims=(0,0.55))

"""
dd_maxent_prob(S::Int64, q::Float64)
    S: number of species
    q: quantile for the number of links using the predictions of the flexible links model
Returns the degree distribution of MaxEnt (analytical solution) for a network of S species and qth quantile of the number of links from the predictions of the flexible links model
"""
function dd_maxent_prob(S::Int64, q::Float64)
  L = Int64(round(quantile(predicted_links[:, S-4], q)))
  dd_maxent = dd_maxent_prob(S, L)
  return dd_maxent
end

# Get the degree distribution of MaxEnt (analytical solution) for a network with a median number of species and different numbers of links
dd015 = dd_maxent_prob(species500, 0.015)
dd055 = dd_maxent_prob(species500, 0.055)
dd165 = dd_maxent_prob(species500, 0.165)
dd500 = dd_maxent_prob(species500, 0.5)
dd835 = dd_maxent_prob(species500, 0.835)
dd945 = dd_maxent_prob(species500, 0.945)
dd985 = dd_maxent_prob(species500, 0.985)

# Plot different degree distributions of MaxEnt (analytical solution) for a network with a median number of species
a = [1,2,5,10,15,20,25] # specified x-ticks

plotB = plot(dd015, color=:black, alpha=0.7, linewidth=2, linestyle=:dot, label="97% PI", # 97% PI
    framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
    guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
    foreground_color_legend=nothing, background_color_legend=:white, legendfont=fonts)
plot!(dd055, color=:black, alpha=0.7, linestyle=:dash, label="89% PI") # 89% PI
plot!(dd165, color=:black, alpha=0.7, linestyle=:solid, label="67% PI") # 67% PI
plot!(dd835, color=:grey, alpha=0.7, linestyle=:solid, label="") # 67% PI
plot!(dd945, color=:grey, alpha=0.7, linestyle=:dash, label="") # 89% PI
plot!(dd985, color=:grey, alpha=0.7, linestyle=:dot, label="") # 97% PI
plot!(dd500, color=:darkblue, linewidth=4, label="median")
xaxis!(:log, xlabel="Degree k", xticks=(a, a), xlims=(1,species500))
yaxis!(ylabel="p(k)")

plot(plotA, plotB,
    title = ["(a)" "(b)"],
    titleloc=:right, titlefont=fonts)

savefig(joinpath("figures","maxent_degree_distributions"))
###########################################################

################# TK TO DO #########################
# Figure: Measures of empirical and maximum entropy food webs

# Convert food webs archieved on Mangal to UnipartiteNetworks
Ns_mangal = network.(mangal_foodwebs.id)
Ns_mangal = convert.(UnipartiteNetwork, Ns_mangal)

# Measure species richness
S_mangal = richness.(Ns_mangal)
S_maxent_empL = richness.(Ns_maxent_empL)
S_maxent_fl = richness.(Ns_maxent_fl)

# Measure nestedness (spectral radius of the adjacency matrix)
nestedness_mangal = ρ.(Ns_mangal)
nestedness_maxent_empL = ρ.(Ns_maxent_empL)
nestedness_maxent_fl = ρ.(Ns_maxent_fl)

# Measure maximum trophic levels
tl_mangal = values.(trophic_level.(Ns_mangal))
tl_max_mangal = maximum.(tl_mangal)

tl_maxent_empL = values.(trophic_level.(Ns_maxent_empL))
tl_max_maxent_empL = maximum.(tl_maxent_empL)

tl_maxent_fl = values.(trophic_level.(Ns_maxent_fl))
tl_max_maxent_fl = maximum.(tl_maxent_fl)

# Measure network diameter (shortest distance between the two most distant nodes in the network)
shortest_paths_mangal = shortest_path.(Ns_mangal)
diameter_mangal = maximum.(shortest_paths_mangal)

shortest_paths_maxent_empL = shortest_path.(Ns_maxent_empL)
diameter_maxent_empL = maximum.(shortest_paths_maxent_empL)

shortest_paths_maxent_fl = shortest_path.(Ns_maxent_fl)
diameter_maxent_fl = maximum.(shortest_paths_maxent_fl)

entropy_mangal = svd_entropy.(Ns_mangal)
entropy_maxent_empL = svd_entropy.(Ns_maxent_empL)
entropy_maxent_fl = svd_entropy.(Ns_maxent_fl)

# Make vector of measures of MaxEnt food webs (nb of links given by the flexible links model) the same length as the other vectors
n1 = length(Ns_mangal)
n2 = length(Ns_maxent_fl)
i = [findall(in(mangal_foodwebs.S[i]).(S_maxent_fl)) for i in 1:n1]
i[n1-1] = [n2] # The second to last entry has no value. We'll change it manually.
i = reduce(vcat, i)

nestedness_maxent_fl_long = [nestedness_maxent_fl[j] for j in i]
tl_max_maxent_fl_long = [tl_max_maxent_fl[j] for j in i]
diameter_maxent_fl_long = [diameter_maxent_fl[j] for j in i]
entropy_maxent_fl_long = [entropy_maxent_fl[j] for j in i]

# Plot predicted as a function of empirical measures
plotA = scatter(nestedness_mangal, nestedness_maxent_empL, alpha=0.3, markersize=3, label="Empirical L",
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:bottomright, legendfont=fonts,
      xlabel="Nestedness (empirical webs)",
      ylabel="Nestedness (MaxEnt webs)")
scatter!(nestedness_mangal, nestedness_maxent_fl_long, alpha=0.3, markersize=3, label="Median L")
plot!(LinRange(0.4, 0.9, 100), LinRange(0.4, 0.9, 100), lab="", color="grey", linestyle=:dot)

plotB = scatter(tl_max_mangal, tl_max_maxent_empL, alpha=0.3, markersize=3, label="Empirical L",
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:bottomright, legendfont=fonts,
      ticks=[2,4,6,8,10,12],
      xlabel="Maximum trophic level (empirical webs)",
      ylabel="Maximum trophic level (MaxEnt webs)")
scatter!(tl_max_mangal, tl_max_maxent_fl_long, alpha=0.3, markersize=3, label="Median L")
plot!(LinRange(2, 12, 100), LinRange(2, 12, 100), lab="", color="grey", linestyle=:dot)

plotC = scatter(diameter_mangal, diameter_maxent_empL, alpha=0.3, markersize=3, label="Empirical L",
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:bottomright, legendfont=fonts,
      ticks=[1,3,5,7,9,11,13],
      xlabel="Diameter (empirical webs)",
      ylabel="Diameter (MaxEnt webs)")
scatter!(diameter_mangal .+ 0.1, diameter_maxent_fl_long, alpha=0.3, markersize=3, label="Median L")
plot!(LinRange(1, 13, 100), LinRange(1, 13, 100), lab="", color="grey", linestyle=:dot)

plotD = scatter(entropy_mangal, entropy_maxent_empL, alpha=0.3, markersize=3, label="Empirical L",
framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
foreground_color_legend=nothing, background_color_legend=:white,
legend=:bottomright, legendfont=fonts,
xlabel="SVD-entropy (empirical webs)",
ylabel="SVD-entropy (MaxEnt webs)")
scatter!(entropy_mangal, entropy_maxent_fl_long, alpha=0.3, markersize=3, label="Median L")
plot!(LinRange(0.77, 1, 100), LinRange(0.77, 1, 100), lab="", color="grey", linestyle=:dot)

plot(plotA, plotB, plotC, plotD,
     title = ["(a)" "(b)" "(c)" "(d)"],
     titleloc=:right, titlefont=fonts)

savefig(joinpath("figures", "measures_mangal_maxent.png"))
#########################################################

################# TK TO DO #########################
# Figure: Plot measures as a function of species richness and of other measures
a = [5,10,20,50,100,400] # specified x-ticks

plotA = scatter(S_mangal, nestedness_maxent_empL, alpha=0.3, label="Empirical L", smooth=true, markersize=3, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      xlabel="Species richness",
      ylabel="Nestedness")
scatter!(S_mangal, nestedness_maxent_fl_long, alpha=0.3, label="Median L", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
scatter!(S_mangal, nestedness_mangal, alpha=0.3, label="Mangal", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
xaxis!(:log, xticks=(a,a))

plotB = scatter(S_mangal, tl_max_maxent_empL, alpha=0.3, label="Empirical L", smooth=true, markersize=3, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      xlabel="Species richness",
      ylabel="Maximum trophic level",
      yticks=[2,4,6,8,10,12])
scatter!(S_mangal, tl_max_maxent_fl_long, alpha=0.3, label="Median L", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
scatter!(S_mangal, tl_max_mangal, alpha=0.3, label="Mangal", smooth=true,markersize=3, linestyle=:dot, linealpha=1)
xaxis!(:log, xticks=(a,a))

plotC = scatter(S_mangal, diameter_maxent_empL, alpha=0.3, label="Empirical L", smooth=true, markersize=3, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      xlabel="Species richness",
      ylabel="Diameter",
      yticks=[1,3,5,7,9,11,13])
scatter!(S_mangal, diameter_maxent_fl_long, alpha=0.3, label="Median L", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
scatter!(S_mangal, diameter_mangal, alpha=0.3, label="Mangal", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
xaxis!(:log, xticks=(a,a))

plotD = scatter(S_mangal, entropy_maxent_empL, alpha=0.3, label="Empirical L", smooth=true, markersize=3, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:bottomright, legendfont=fonts,
      xlabel="Species richness",
      ylabel="SVD-entropy")
scatter!(S_mangal, entropy_maxent_fl_long, alpha=0.3, label="Median L", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
scatter!(S_mangal, entropy_mangal, alpha=0.3, label="Mangal", smooth=true, markersize=3, linestyle=:dot, linealpha=1)
xaxis!(:log, xticks=(a,a))

plot(plotA, plotB, plotC, plotD,
     title = ["(a)" "(b)" "(c)" "(d)"],
     titleloc=:right, titlefont=fonts)

savefig(joinpath("figures", "measures_richness.png"))

# Plot nestedness as a function of maximum entropy
scatter(tl_max_maxent_empL, nestedness_maxent_empL, alpha=0.5, label="Empirical L", smooth=true, markersize=4, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      xticks=[2,4,6,8,10,12],
      xlabel="Maximum trophic level",
      ylabel="Nestedness")
scatter!(tl_max_maxent_fl, nestedness_maxent_fl, alpha=0.5, label="Median L", smooth=true, markersize=4, linestyle=:dot, linealpha=1)
scatter!(tl_max_mangal, nestedness_mangal, alpha=0.5, label="Mangal", smooth=true, markersize=4, linestyle=:dot, linealpha=1)

savefig(joinpath("figures", "maxtrophiclevel_nestedness.png"))
#############################################################

################# TK TO DO #########################
## Figure: Entropy and divergence of degree distributions
"""
get_divergence_dd(N1::UnipartiteNetwork{Bool, MangalNode}, N2::UnipartiteNetwork{Bool, String})
    N1: Empirical food web
    N2: MaxEnt food web
Returns the divergence in probabilistic degree distributions between the two networks
"""
function get_divergence_dd(N1::UnipartiteNetwork{Bool, MangalNode}, N2::UnipartiteNetwork{Bool, String})
  # Maximum in species richness between the two networks (since some predicted networks didn't have the exact same number of species as their empirical counterpart)
  S = maximum(vcat(richness(N1), richness(N2)))

  # Get the probabilistic degree distributions (proportion of species with degree k) of the two networks
  dd_N1 = get_ranked_dd(N1)
  dd_N2 = get_ranked_dd(N2)

  dd_prob_N1 = zeros(Float64, S)
  dd_prob_N2 = zeros(Float64, S)

  for k in 1:S
    dd_prob_N1[k] = sum(dd_N1 .== k) / S
    dd_prob_N2[k] = sum(dd_N2 .== k) / S
  end

  # Compute the divergence in degree distributions between the two networks (sum of absolute differences)
  divergence = sum(abs.(dd_prob_N1 .- dd_prob_N2))

  return divergence
end

# Compute the divergence in degree distributions between food webs archived on Mangal and MaxEnt food webs (empirical L)
divergence_dd = get_divergence_dd.(Ns_mangal, Ns_maxent_empL)

# Plot the divergence in degree distributions as a function of entropy
plotA = scatter(entropy_maxent_empL, divergence_dd, alpha=0.3, label="", smooth=true, markersize=3, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      xlabel="SVD-entropy (MaxEnt food web)",
      ylabel="Divergence in degree distributions")

# Plot the divergence in degree distributions as a function of species richness
plotB = scatter(S_mangal, divergence_dd, alpha=0.3, label="", smooth=true, markersize=3, linestyle=:dot, linealpha=1,
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      xlabel="Species richness",
      ylabel="Divergence in degree distributions")
xaxis!(:log, xticks=(a,a))

plot(plotA, plotB,
     title = ["(a)" "(b)"],
     titleloc=:right, titlefont=fonts)

savefig(joinpath("figures", "divergence_degree_distributions.png"))
################################################################################

################# TK TO DO #########################
## Figure: Distribution of entropy and z-scores of empirical food webs

plotA = density(entropy_maxent_empL, label="Empirical L",
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      ylims=(0,30),
      xlabel="SVD-entropy",
      ylabel="Density")
density!(entropy_maxent_fl, label="Median L")
density!(entropy_mangal, label="Mangal")

entropy_maxent_fl_avg = mean(entropy_maxent_fl)
entropy_maxent_fl_std = std(entropy_maxent_fl)
entropy_mangal_zscores = (entropy_mangal .- entropy_maxent_fl_avg) ./ entropy_maxent_fl_std

entropy_mangal_zscores_500 = quantile(entropy_mangal_zscores, 0.500)
entropy_mangal_zscores_015 = quantile(entropy_mangal_zscores, 0.015)
entropy_mangal_zscores_985 = quantile(entropy_mangal_zscores, 0.985)

plotB = density(entropy_mangal_zscores, label="",
      framestyle=:box, dpi=1000, size=(800,500), margin=5Plots.mm,
      guidefont=fonts, xtickfont=fonts, ytickfont=fonts,
      foreground_color_legend=nothing, background_color_legend=:white,
      legend=:topright, legendfont=fonts,
      color=:grey,
      ylims=(0,0.16),
      xlabel="z-score of SVD-entropy",
      ylabel="Density")
plot!([entropy_mangal_zscores_500], seriestype=:vline, color=:grey, ls=:dash, lab="")

plot(plotA, plotB,
     title = ["(a)" "(b)"],
     titleloc=:right, titlefont=fonts)

savefig(joinpath("figures", "entropy_distribution.png"))
#####################################################################


text:
module DistributionsAD

using PDMats,
      ForwardDiff,
      Zygote,
      LinearAlgebra,
      Distributions,
      Random,
      Combinatorics,
      SpecialFunctions,
      StatsFuns,
      Compat

using Tracker: Tracker, TrackedReal, TrackedVector, TrackedMatrix, TrackedArray,
                TrackedVecOrMat, track, @grad, data
using SpecialFunctions: logabsgamma, digamma
using ZygoteRules: ZygoteRules, @adjoint, pullback
using LinearAlgebra: copytri!, AbstractTriangular
using Distributions: AbstractMvLogNormal,
                     ContinuousMultivariateDistribution
using DiffRules, SpecialFunctions, FillArrays
using ForwardDiff: @define_binary_dual_op # Needed for `eval`ing diffrules here
using Base.Iterators: drop
using ReverseDiff

import StatsFuns: logsumexp,
                  binomlogpdf,
                  nbinomlogpdf,
                  poislogpdf,
                  nbetalogpdf
import Distributions: MvNormal,
                      MvLogNormal,
                      poissonbinomial_pdf_fft,
                      logpdf,
                      quantile,
                      PoissonBinomial

export TuringScalMvNormal,
       TuringDiagMvNormal,
       TuringDenseMvNormal,
       TuringMvLogNormal,
       TuringPoissonBinomial,
       TuringWishart,
       TuringInverseWishart,
       arraydist,
       filldist

include("common.jl")
include("univariate.jl")
include("multivariate.jl")
include("matrixvariate.jl")
include("flatten.jl")
include("arraydist.jl")
include("filldist.jl")
include("reversediff.jl")

end


text:
using Test
using LogicCircuits

include("../helper/little_circuits.jl")

@testset "SDD as LogicCircuit parser" begin

  circuit = zoo_sdd("random.sdd")

  @test circuit isa PlainLogicCircuit

  @test num_nodes(circuit) == 1676
  @test isdecomposable(circuit)
  @test !issmooth(circuit)
  @test any(isfalse, linearize(circuit))
  @test any(istrue, linearize(circuit))
  @test num_variables(circuit) == 30

  # not sure what this is testing for
  @test issomething(canonical_literals(circuit))
  @test issomething(canonical_constants(circuit))

end

@testset "SDD write and save" begin

  sdd = readme_sdd()

  mktempdir() do tmp

    # write as a unstructured logic circuit
    sdd_path = "$tmp/example.sdd"
    write(sdd_path, sdd)

    # read as a unstructured logic circuit
    sdd2 = read(sdd_path, LogicCircuit, SddFormat())

    @test sdd2 isa PlainLogicCircuit

    @test num_nodes(sdd) == num_nodes(sdd2)
    @test prob_equiv(sdd, sdd2, 10)

    # write with vtree
    vtree_path = "$tmp/example.vtree"
    paths = (sdd_path, vtree_path)
    write(paths, sdd)

    # read as a structured logic circuit
    formats = (SddFormat(), VtreeFormat())
    sdd3 = read(paths, StructLogicCircuit, formats)

    @test sdd3 isa PlainStructLogicCircuit

    @test num_nodes(sdd) == num_nodes(sdd3)
    @test prob_equiv(sdd, sdd3, 10)

    @test Vtree(mgr(sdd)) == vtree(sdd3)

    # read as an SDD
    sdd4 = read(paths, Sdd, formats)

    @test sdd4 isa Sdd

    @test num_nodes(sdd) == num_nodes(sdd4)
    @test prob_equiv(sdd, sdd4, 10)

    @test Vtree(mgr(sdd)) == Vtree(mgr(sdd4))

    # write/read compressed
    paths = ("$sdd_path.gz", vtree_path)
    write(paths, sdd)
    sdd4 = read(paths, StructLogicCircuit)

    @test sdd4 isa PlainStructLogicCircuit
    @test num_nodes(sdd) == num_nodes(sdd4)
    @test prob_equiv(sdd, sdd4, 10)
    @test Vtree(mgr(sdd)) == vtree(sdd4)

  end

  mktempdir() do tmp
    mgr = SddMgr(7, :balanced)
    v = Dict([(i => compile(mgr, Lit(i))) for i=1:7])
    c = (v[1] | !v[2] | v[3]) &
        (v[2] | !v[7] | v[6]) &
        (v[3] | !v[4] | v[5]) &
        (v[1] | !v[4] | v[6])
    @test num_edges(c) == 147

    @test_nowarn write("$tmp/temp.sdd", c)
    l = read("$tmp/temp.sdd", LogicCircuit)
    @test num_edges(l) == 147
    @test num_variables(l) == 7

  end

end

max_issues_repo_name: Juice-jl/LogicCircuits

max_forks_repo_name: Juice-jl/LogicCircuits


text:
using Colors

export normalcolors, lambert, basic, phong, visible, edges

abstract MaterialPrimitive

immutable Material{P<:MaterialPrimitive} <: Compose3DNode
    primitives::Vector{P}
end

isscalar(material::Material) = length(material.primitives) == 1

immutable MeshColor <: MaterialPrimitive
    color::Color
end

function mesh_color(color::Color)
    Material([MeshColor(color)])
end

function mesh_color(colorString::AbstractString)
    Material([MeshColor(parse(Colorant, colorString))])
end

immutable WireFrameMesh <: MaterialPrimitive
    wireframe::Bool
end

function wireframe(wireframe::Bool)
    Material([WireFrameMesh(wireframe)])
end

immutable MaterialKind <: MaterialPrimitive
    kind::AbstractString
end

function normalcolors()
    Material([MaterialKind("normal")])
end

function lambert()
    Material([MaterialKind("lambert")])
end

function basic()
    Material([MaterialKind("basic")])
end

function phong()
    Material([MaterialKind("phong")])
end

immutable Visiblity <: MaterialPrimitive
    visible::Bool
end

function visible(visible::Bool)
    Material([Visiblity(visible)])
end

immutable Edges <: MaterialPrimitive
    color::Color
end

function edges(color::Color)
    Material([Edges(color)])
end

function edges(colorString::AbstractString="black")
    Material([Edges(parse(Colorant, colorString))])
end

function resolve(box::Absolute3DBox, prop::Material)
    prop
end


text:
@testset "Words" begin

    @test Word(Int[]) isa KnuthBendix.AbstractWord
    @test Word(Int[]) isa Word
    @test Word(Int[]) isa Word{UInt16}
    @test one(Word{UInt16}) isa Word{UInt16}
    @test one(Word{Int}) isa Word{Int}

    w = Word([1, 2])
    W = Word{Int}([1, 2])
    @test w isa Word{UInt16}
    @test W isa Word{Int}

    @test w == Word([1, 2])
    @test w == W
    @test w !== Word([1, 2])
    @test w != one(w)
    @test !isone(w)
    @test isone(one(w))
    @test hash(w) isa UInt
    @test hash(w) == hash(W)
    @test hash(w, UInt(1)) != hash(w, UInt(0))
    @test hash(w) != hash(one(w))
    @test_throws AssertionError Word([1, -2])

    @test push!(w, 3) == Word([1, 2, 3])
    @test w[3] == 3
    @test w == Word([1, 2, 3])
    @test_throws AssertionError push!(w, -3)

    @test pushfirst!(w, 4) == Word([4, 1, 2, 3])
    @test w[1] == 4
    @test w == Word([4, 1, 2, 3])
    @test_throws AssertionError pushfirst!(w, -4)

    @test w[2:3] isa Word{UInt16}
    @test w[2:3] == Word([1,2])
    @test W[2:2] isa Word{Int}
    @test W[2:2][1] == W[2]

    @test append!(w, Word([6])) == Word([4, 1, 2, 3, 6])
    @test prepend!(w, Word([5])) == Word([5, 4, 1, 2, 3, 6])
    @test collect(w) == [5, 4, 1, 2, 3, 6]
    @test collect(w) isa Vector{UInt16}
    @test w[1] == 5
    @test length(w) == 6
    @test_throws BoundsError w[-1]
    @test_throws BoundsError w[7]

    @test similar(w) isa Word{UInt16}
    @test similar(W) isa Word{Int}

    @test Word([1, 2]) * Word([2, 3]) == Word([1, 2, 2, 3])
    @test Word([1, 2]) * W == Word([1, 2, 1, 2])

    u1 = Word([1,2,3,4])
    u2 = Word([1,2])
    u3 = Word([4,1,2,3,4])

    @test pop!(u2) == 2
    @test u2 == Word([1])
    @test popfirst!(u3) == 4
    @test u3 == u1

    @test sprint(show, Word([1,2])) == "1·2"
    @test sprint(show, MIME"text/plain"(), Word([1,2])) == "Word{UInt16}: 1·2"
    @test sprint(show, one(Word{UInt16})) == "(id)"
end

@testset "SubWords" begin
    w = Word([1,2,3])
    @test @view(w[1:2]) isa KnuthBendix.SubWord
    vw = @view w[2:3]
    @test vw == w[2:3]
    @test vw[1] == w[2]
    vw[1] = 10
    @test vw[1] == w[2] == 10

    @test @view(vw[1:1]) isa KnuthBendix.SubWord

    u = Word([5,6,7])
    v = @view deepcopy(u)[2:3]
    @test u*v == Word([5,6,7,6,7])
    @test append!(u, v) isa Word
    @test prepend!(u,v) isa Word

    @test isone(@view(u[1:0]))
end

abstract_word_conformance_test(Word)


text:
using FinancialMonteCarlo
path1=joinpath(dirname(pathof(FinancialMonteCarlo)),"..","test","af")
test_listTmp=readdir(path1);
BlackList=["REQUIRE","runtests.jl","runner.jl","test_black_cuda.jl","cuda","af"]
test_list=[test_element for test_element in test_listTmp if !Bool(sum(test_element.==BlackList))]
println("Running tests:\n")
for (current_test,i) in zip(test_list,1:length(test_list))
	println("------------------------------------------------------------")
    println("  * $(current_test) *")
    include(joinpath(path1,current_test))
	println("------------------------------------------------------------")
	if (i<length(test_list))
		println("")
	end
end


text:
module ExploreFunction

include("MCMCMapping.jl")
include("VietorisRips.jl")

end


text:
include("peps_util.jl")

function reindexA(A::fPEPS, new_A::fPEPS)
    Ny, Nx   = size(new_A)
    midpoint = div(Nx, 2)
    old_column = Vector{Int}(undef, Nx)
    old_column[1:midpoint-1]  = 1:midpoint-1
    old_column[midpoint]      = midpoint
    old_column[midpoint+1]    = midpoint - 1
    old_column[midpoint+2:Nx] = midpoint:Nx-2
    new_right_bonds = Matrix{Index}(undef, Ny, Nx)
    new_left_bonds  = Matrix{Index}(undef, Ny, Nx)
    new_up_bonds    = Matrix{Index}(undef, Ny, Nx)
    new_down_bonds  = Matrix{Index}(undef, Ny, Nx)
    for col in 1:Nx, row in 1:Ny
        oldSpin = findindex(new_A[row, col], "Site")
        newSpin = Index(dim(oldSpin), tags(oldSpin))
        replaceind!(new_A[row, col], oldSpin, newSpin)
        if row != Ny
            oldUp = commonind(A[row, old_column[col]], A[row+1, old_column[col]])
            newUp = Index(dim(oldUp), "Link,u,c$col,r$row")
            new_up_bonds[row, col] = newUp
            new_down_bonds[row+1, col] = newUp
            replaceind!(new_A[row, col], oldUp, newUp)
        end
        if row != 1
            oldDown = commonind(A[row, old_column[col]], A[row-1, old_column[col]])
            newDown = new_down_bonds[row, col]
            replaceind!(new_A[row, col], oldDown, newDown)
        end
        if col != Nx
            oldRight = commonind(A[row, old_column[col]+1], A[row, old_column[col]])
            newRight = Index(dim(oldRight), "Link,r,c$col,r$row")
            new_right_bonds[row, col]  = newRight
            new_left_bonds[row, col+1] = newRight
            replaceind!(new_A[row, col], oldRight, newRight)
        end
        if col != 1
            oldLeft = commonind(A[row, old_column[col]-1], A[row, old_column[col]])
            newLeft = new_left_bonds[row, col]
            replaceind!(new_A[row, col], oldLeft, newLeft)
        end
    end
    return new_A
end

function InsertTensors(A::fPEPS, Rs::Vector{Environments}, Ls::Vector{Environments}, H; kwargs...)
    maxdim      = get(kwargs, :maxdim, 1)
    Ny, Nx      = size(A)
    midpoint    = div(Nx, 2)
    A, Q_f, R_f, finds = gaugeColumnForInsert(A, midpoint, :right; kwargs...)
    for row in 1:Ny
        QRind = findindex(Q_f[row], "QR")
        ci    = commonind(A[row, midpoint], A[row, midpoint + 1])
        replaceind!(Q_f[row], QRind, ci)
        A[row, midpoint] = Q_f[row]
    end
    oldH = deepcopy(H)
    Ar = deepcopy(A[:, midpoint + 1])
    Al = deepcopy(A[:, midpoint])
    new_A = fPEPS(Nx+2, Ny, hcat(deepcopy(A[:, 1:midpoint]), Ar, Al, deepcopy(A[:, midpoint+1:Nx])))
    new_A = reindexA(A, new_A)
    A     = new_A
    Hr = H[:, midpoint + 1]
    Hl = H[:, midpoint]
    H = hcat(H[:, 1:midpoint], Hr, Hl, H[:, midpoint+1:Nx])
    Nx += 2
    A, Q_l, R_l, linds = gaugeColumnForInsert(A, midpoint + 1, :right; kwargs...)
    A, Q_r, R_r, rinds = gaugeColumnForInsert(A, midpoint + 2, :left; kwargs...)
    A[:, midpoint + 1] = tensors(Q_l)
    A[:, midpoint + 2] = tensors(Q_r)
    Λı = deepcopy(R_l) # Lambdan
    Λȷ = MPO(Ny) # Lambdanmo
    X  = MPO(Ny)
    X_up_inds = [Index(maxdim, "Link,Xup,l=$row") for row in 1:Ny-1]
    midpoint += 1
    for row in 1:Ny
        Λȷ[row] = R_f[row]
        Xinds   = Index[findindex(R_l[row], "QR"); rinds[row]]
        if row > 1
            push!(Xinds, X_up_inds[row-1])
        end
        if row < Ny
            push!(Xinds, X_up_inds[row])
        end
        X[row]  = cuITensor(dense(diagITensor(1.0, IndexSet(Xinds...))))
        X[row] /= norm(X[row])
        i_ind = findindex(Λı[row], "DM")
        j_ind = findindex(Λȷ[row], "QR")
        replaceind!(Λȷ[row], j_ind, i_ind)
        replaceind!(Λȷ[row], finds[row], rinds[row])
    end
    γ = 0.05
    inner_product = 0.1
    iter   = 1
    Λk     = deepcopy(Λȷ) #Lambdanmo_
    Λℓ     = deepcopy(Λȷ) #Lambdanmo__
    Λm     = deepcopy(Λı) #Lambdan_
    # F = ||A*x - b||^2 - ||x||^2
    # Grad F = 2A^T(Ax - b) - 2x
    # grad_const is 2A^T*b
    #foreach(println, inds.(Λk))
    #foreach(println, inds.(Λm))
    Δconst = multMPO(Λk, Λm; kwargs...)
    println("Begin grad descent")
    flush(stdout)
    while abs(inner_product) < 0.95 && iter < 100
        tmp_Δ  = multMPO(X, Λȷ; kwargs...)
        tmp_Δ2 = multMPO(tmp_Δ, Λℓ; kwargs...)
        Δ      = sum([2.0*Δconst, -2.0*X, 2.0*tmp_Δ2])
        X      = sum(X, -γ*Δ; kwargs...)
        tmp_Δ  = multMPO(Λℓ, X; kwargs...)
        inner_product = overlap(tmp_Δ, Λm)
        iter  += 1
    end
    println("Done grad descent")
    flush(stdout)
    XR  = multMPO(X, R_r; kwargs...)
    # must combine inds of QR
    Q_r_cmbs = Vector{ITensor}(undef, Ny)
    for row in 1:Ny
        cmb_inds = IndexSet(uniqueindex(findinds(Q_r[row], "Site"), findinds(XR[row], "Site")), findindex(Q_r[row], "Link,r"))
        Q_r_cmb, ci = combiner(cmb_inds, tags="Qcmb,Site,r$row")
        Q_r_cmbs[row] = Q_r_cmb
        Q_r[row] *= Q_r_cmb
    end
    XRQ = multMPO(XR, Q_r; kwargs...)
    for row in 1:Ny
        XRQ[row] *= Q_r_cmbs[row]
    end
    A[:, midpoint+1] = tensors(XRQ)
    for row in 1:Ny
        qr_ci = commonind(A[row, midpoint], A[row, midpoint+1])
        new_link = Index(dim(qr_ci), "Link,r,c$midpoint,r$row")
        replaceind!(A[row, midpoint], qr_ci, new_link)
        replaceind!(A[row, midpoint+1], qr_ci, new_link)
    end
    return A, H
end

function rightwardSweepToInsert(A::fPEPS, Ls::Vector{Environments}, Rs::Vector{Environments}, H; kwargs...)
    simple_update_cutoff = get(kwargs, :simple_update_cutoff, 4)
    Ny, Nx = size(A)
    dummyI = MPS(Ny, fill(ITensor(1.0), Ny), 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny))
    prev_cmb_r = Vector{ITensor}(undef, Ny)
    next_cmb_r = Vector{ITensor}(undef, Ny)
    sweep::Int = get(kwargs, :sweep, 0)
    stop_col::Int = get(kwargs, :stop_col, div(Nx, 2))
    sweep_width::Int = get(kwargs, :sweep_width, Nx)
    offset    = mod(Nx, 2)
    midpoint  = div(Nx, 2)
    rightmost = midpoint
    leftmost  = sweep_width == Nx ? 1 : midpoint - div(sweep_width, 2)
    if leftmost > 1
        for row in 1:Ny
            prev_cmb_r[row] = reconnect(commonind(A[row, leftmost], A[row, leftmost-1]), Ls[leftmost-1].I[row])
        end
    end
    @inbounds for col in leftmost:rightmost
        L = col == 1 ? dummyEnv : Ls[col - 1]
        @debug "Sweeping col $col"
        if sweep >= simple_update_cutoff
            @timeit "sweep" begin
                A = sweepColumn(A, L, Rs[col+1], H, col; kwargs...)
            end
        end
        if sweep < simple_update_cutoff
            # Simple update...
            A = simpleUpdate(A, col, col+1, H; do_side=(col < Nx), kwargs...)
        end
        if sweep >= simple_update_cutoff
            # Gauge
            A = gaugeColumn(A, col, :right; kwargs...)
        end
        if col == 1
            left_H_terms = getDirectional(H[1], Horizontal)
            @timeit "left edge env" begin
                Ls[col] = buildEdgeEnvironment(A, H, left_H_terms, prev_cmb_r, :left, 1; kwargs...)
            end
        else
            @timeit "left next env" begin
                I, H_, IP = buildNextEnvironment(A, Ls[col-1], H, prev_cmb_r, next_cmb_r, :left, col; kwargs...)
                Ls[col] = Environments(I, H_, IP)
            end
            prev_cmb_r = deepcopy(next_cmb_r)
        end
    end
    return A, Ls, Rs
end

function rightwardSweepFromInsert(A::fPEPS, Ls::Vector{Environments}, Rs::Vector{Environments}, H; kwargs...)
    simple_update_cutoff = get(kwargs, :simple_update_cutoff, 4)
    Ny, Nx = size(A)
    dummyI = MPS(Ny, fill(ITensor(1.0), Ny), 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny))
    prev_cmb_r = Vector{ITensor}(undef, Ny)
    next_cmb_r = Vector{ITensor}(undef, Ny)
    sweep::Int = get(kwargs, :sweep, 0)
    start_col::Int = get(kwargs, :start_col, div(Nx, 2))
    sweep_width::Int = get(kwargs, :sweep_width, Nx)
    offset    = mod(Nx, 2)
    midpoint  = div(Nx, 2)
    rightmost = midpoint + div(sweep_width, 2) + offset
    leftmost  = midpoint + 1
    @inbounds for col in leftmost:rightmost
        L = col == 1 ? dummyEnv : Ls[col - 1]
        @debug "Sweeping col $col"
        if sweep >= simple_update_cutoff
            @timeit "sweep" begin
                A = sweepColumn(A, L, Rs[col+1], H, col; kwargs...)
            end
        end
        if sweep < simple_update_cutoff
            # Simple update...
            A = simpleUpdate(A, col, col+1, H; do_side=(col < Nx), kwargs...)
        end
        if sweep >= simple_update_cutoff
            # Gauge
            A = gaugeColumn(A, col, :right; kwargs...)
        end
        if col == 1
            left_H_terms = getDirectional(H[1], Horizontal)
            @timeit "left edge env" begin
                Ls[col] = buildEdgeEnvironment(A, H, left_H_terms, prev_cmb_r, :left, 1; kwargs...)
            end
        else
            @timeit "left next env" begin
                I, H_, IP = buildNextEnvironment(A, Ls[col-1], H, prev_cmb_r, next_cmb_r, :left, col; kwargs...)
                Ls[col] = Environments(I, H_, IP)
            end
            prev_cmb_r = deepcopy(next_cmb_r)
        end
    end
    return A, Ls, Rs
end

function rightwardSweepFromInsert(A::fPEPS, Ls::Vector{Environments}, Rs::Vector{Environments}, H; kwargs...)
    simple_update_cutoff = get(kwargs, :simple_update_cutoff, 4)
    Ny, Nx = size(A)
    dummyI = MPS(Ny, fill(ITensor(1.0), Ny), 0, Ny+1)
    dummyEnv = Environments(dummyI, dummyI, fill(ITensor(), 1, Ny))
    prev_cmb_r = Vector{ITensor}(undef, Ny)
    next_cmb_r = Vector{ITensor}(undef, Ny)
    sweep::Int = get(kwargs, :sweep, 0)
    @inbounds for col in 1:Nx-1
        L = col == 1 ? dummyEnv : Ls[col - 1]
        @debug "Sweeping col $col"
        if sweep >= simple_update_cutoff
            @timeit "sweep" begin
                A = sweepColumn(A, L, Rs[col+1], H, col; kwargs...)
            end
        end
        if sweep < simple_update_cutoff
            # Simple update...
            A = simpleUpdate(A, col, col+1, H; do_side=(col < Nx), kwargs...)
        end
        if sweep >= simple_update_cutoff
            # Gauge
            A = gaugeColumn(A, col, :right; kwargs...)
        end
        if col == 1
            left_H_terms = getDirectional(H[1], Horizontal)
            @timeit "left edge env" begin
                Ls[col] = buildEdgeEnvironment(A, H, left_H_terms, prev_cmb_r, :left, 1; kwargs...)
            end
        else
            @timeit "left next env" begin
                I, H_, IP = buildNextEnvironment(A, Ls[col-1], H, prev_cmb_r, next_cmb_r, :left, col; kwargs...)
                Ls[col] = Environments(I, H_, IP)
            end
            prev_cmb_r = deepcopy(next_cmb_r)
        end
    end
    return A, Ls, Rs
end

function doSweepsInsert(A::fPEPS, Ls::Vector{Environments}, Rs::Vector{Environments}, H; mindim::Int=1, maxdim::Int=1, simple_update_cutoff::Int=4, sweep_start::Int=1, sweep_count::Int=10, cutoff::Float64=0., insert_interval::Int=20)
    Ny, Nx      = size(A)
    sweep_width = Nx
    for sweep in sweep_start:sweep_count
        insert_flag = mod(sweep, insert_interval) == 0
        if iseven(sweep)
            if !insert_flag
                println("SWEEP RIGHT $sweep")
                A, Ls, Rs = rightwardSweep(A, Ls, Rs, H; sweep_width=sweep_width, sweep=sweep, mindim=mindim, maxdim=maxdim, simple_update_cutoff=simple_update_cutoff, overlap_cutoff=0.999, cutoff=cutoff)
            else
                println("SWEEP RIGHT TO INSERT $sweep")
                A, Ls, Rs = rightwardSweepToInsert(A, Ls, Rs, H; sweep_width=sweep_width, stopcol=div(Nx, 2), sweep=sweep, mindim=mindim, maxdim=maxdim, simple_update_cutoff=simple_update_cutoff, overlap_cutoff=0.999, cutoff=cutoff)
            end
        else
            println("SWEEP LEFT $sweep")
            A, Ls, Rs = leftwardSweep(A, Ls, Rs, H; sweep_width=sweep_width, sweep=sweep, mindim=mindim, maxdim=maxdim, simple_update_cutoff=simple_update_cutoff, overlap_cutoff=0.999, cutoff=cutoff)
        end
        flush(stdout)
        if sweep == simple_update_cutoff - 1
            for col in reverse(2:Nx)
                A = gaugeColumn(A, col, :left; mindim=1, maxdim=maxdim, cutoff=cutoff)
            end
            Ls = buildLs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
            Rs = buildRs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
        end
        if insert_flag
            A, H = InsertTensors(A, Rs, Ls, H; mindim=mindim, maxdim=maxdim, overlap_cutoff=0.999, cutoff=cutoff)
            Ny, Nx = size(A)
            println(" Insert complete at Nx = $Nx, Ny = $Ny")
            flush(stdout)
            sweep_width = Nx
            Ls = buildLs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
            Rs = buildRs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
            println(" Rebuild of envs complete at Nx = $Nx, Ny = $Ny")
            flush(stdout)
            rightwardSweepFromInsert(A, Ls, Rs, H; sweep_width=sweep_width, startcol=div(Nx, 2), sweep=sweep, mindim=mindim, maxdim=maxdim, simple_update_cutoff=simple_update_cutoff, overlap_cutoff=0.999, cutoff=cutoff)
            sweep_width = Nx
            if sweep_width < Nx
                midpoint = div(Nx, 2)
                barrier = midpoint + div(sweep_width, 2) + 2
                barrier = barrier > Nx ? Nx : barrier
                for col in reverse((midpoint + div(sweep_width, 2) + 1):barrier)
                    A = gaugeColumn(A, col, :left; mindim=1, maxdim=maxdim, cutoff=cutoff)
                end
                #Ls = buildLs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
                #Rs = buildRs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
            end
        end
        Ls = buildLs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
        Rs = buildRs(A, H; mindim=1, maxdim=maxdim, cutoff=cutoff)
        #=if sweep == sweep_count
            A_ = deepcopy(A)
            L_s = buildLs(A_, H; mindim=mindim, maxdim=maxdim)
            R_s = buildRs(A_, H; mindim=mindim, maxdim=maxdim)
            x_mag = measureXmag(A_, L_s, R_s; mindim=mindim, maxdim=maxdim)
            z_mag = measureZmag(A_, L_s, R_s; mindim=mindim, maxdim=maxdim)
            display(z_mag)
            println()
            v_mag = measureSmagVertical(A_, L_s, R_s; mindim=mindim, maxdim=maxdim)
            display(v_mag)
            println()
            h_mag = measureSmagHorizontal(A_, L_s, R_s; mindim=mindim, maxdim=maxdim)
            display(h_mag)
            println()
        end=#
    end
    return A
end