text:
export getFacesStaggeredIndicesOfCellNoPressure,getFacesStaggeredIndicesOfCell,getCellCenteredIndicesOfCell,getNodalIndicesOfCell,getSubMeshOfCell

export getDirichletMassNodal

function getBoxIdxs(UpperLeftCorner::Array{Int64},BottomRightCorner::Array{Int64})
	# this is a type of an ndgrid function
	if length(UpperLeftCorner)==2
		i1 = vec(UpperLeftCorner[1]:BottomRightCorner[1]);
		i2 = vec(UpperLeftCorner[2]:BottomRightCorner[2]);
		I1 = i1*ones(Int64,length(i2))';
		I2 = ones(Int64,length(i1))*i2';
		return I1,I2
	else
		vx = UpperLeftCorner[1]:BottomRightCorner[1];
		vy = UpperLeftCorner[2]:BottomRightCorner[2];
		vz = UpperLeftCorner[3]:BottomRightCorner[3];
		m, n, o = length(vy), length(vx), length(vz)
		vx = reshape(vx, 1, n, 1)
		vy = reshape(vy, m, 1, 1)
		vz = reshape(vz, 1, 1, o)
		om = ones(Int64, m)
		on = ones(Int64, n)
		oo = ones(Int64, o)
		return (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
	end
end

function getOriginalBoundingBoxCells(NumCells::Array{Int64,1},overlap::Array{Int64,1},i,nc)

	cellSize = div.(nc,NumCells);
	originalUpperLeftCorner = (i.-1).*cellSize .+ 1;
	originalBottomRightCorner = originalUpperLeftCorner .+ (cellSize.-1);
	originalBottomRightCorner[i.==NumCells] .== nc[i.==NumCells];
	return (originalUpperLeftCorner,originalBottomRightCorner);
	# if length(nc)==2
		# originalUpperLeftCorner = [(i[1]-1)*cellSize[1] + 1,(i[2]-1)*cellSize[2] + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [cellSize[1]-1,cellSize[2]-1];
	# else
		# originalUpperLeftCorner = [(i[1]-1)*cellSize[1] + 1,(i[2]-1)*cellSize[2] + 1,(i[3]-1)*cellSize[3] + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [cellSize[1]-1,cellSize[2]-1,cellSize[3]-1];
	# end
end

##
# This function just extends the box with overlapp only if possible.
##

function getBoxWithOverlap(UpperLeftCorner::Array{Int64},BottomRightCorner::Array{Int64},nc::Array{Int64},overlap::Array{Int64})
	if length(nc)==2
		newUpperLeft = copy(UpperLeftCorner);
		newBottomRight = copy(BottomRightCorner);
		if UpperLeftCorner[1] > 1
			newUpperLeft[1]-=overlap[1];
		end
		if UpperLeftCorner[2] > 1
			newUpperLeft[2]-=overlap[2];
		end
		if BottomRightCorner[1] < nc[1]
			newBottomRight[1]+=overlap[1];
		end
		if BottomRightCorner[2] < nc[2]
			newBottomRight[2]+=overlap[2];
		end
	else
		newUpperLeft = copy(UpperLeftCorner);
		newBottomRight = copy(BottomRightCorner);
		for i=1:3
			if UpperLeftCorner[i] > 1
				newUpperLeft[i]-=overlap[i];
			end
		end
		for i=1:3
			if BottomRightCorner[i] < nc[i]
				newBottomRight[i]+=overlap[i];
			end
		end
	end
	return (newUpperLeft,newBottomRight);
end

function getCellCenteredIndicesOfCell(NumCells::Array{Int64,1},overlap::Array{Int64,1},i,nc)
	indOver = [];
	if length(nc)==2
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner,nc,overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indOver = I1[:] .+ (I2[:].-1)*(nc[1])
	else
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1,(i[3]-1)*div(nc[3],NumCells[3]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1,div(nc[3],NumCells[3])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# CELLS X CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner,nc,overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indOver = I1[:] .+ (I2[:].-1)*(nc[1]) .+ (I3[:].-1)*(nc[1]*nc[2]);
	end
	return indOver;
end

function getSubMeshOfCell(NumCells::Array{Int64,1},overlap::Array{Int64,1},i,Mesh::RegularMesh)
	indOver = [];
	nc = Mesh.n;
	if length(nc)==2
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner,nc,overlap);
	else
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1,(i[3]-1)*div(nc[3],NumCells[3]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1,div(nc[3],NumCells[3])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# CELLS X CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner,nc,overlap);
	end
	subDomain = copy(Mesh.domain);
	subDomain[1:2:end] = Mesh.domain[1:2:end] .+ (UpperLeftCorner .- 1).*Mesh.h;
	subDomain[2:2:end] = Mesh.domain[2:2:end] .- (Mesh.n - BottomRightCorner).*Mesh.h;
	subMesh = getRegularMesh(subDomain,BottomRightCorner - UpperLeftCorner .+ 1);
	return subMesh
end

function getNodalIndicesOfCell(NumCells::Array{Int64,1},overlap::Array{Int64,1},i,nc)

	if length(nc)==2
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1]),div(nc[2],NumCells[2])];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner.+ [1;1],nc .+ [1;1],overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indNodal = I1[:] + (I2[:].-1)*(nc[1]+1)
	else
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1,(i[3]-1)*div(nc[3],NumCells[3]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1]),div(nc[2],NumCells[2]),div(nc[3],NumCells[3])];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# NODES X NODES X NODES
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner.+ [1;1],nc + [1;1;1],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indNodal = I1[:] + (I2[:].-1)*(nc[1]+1) + (I3[:].-1)*((nc[1]+1)*(nc[1]+1));
	end

	return indNodal
end

function getDirichletMassNodal(NumCells::Array{Int64},overlap::Array{Int64},i::Array{Int64},nc::Array{Int64})
	if length(nc)==2
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1]),div(nc[2],NumCells[2])];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner .+ [1;1],nc .+ [1;1],overlap);
		size = BottomRightCorner-UpperLeftCorner.+1;
		mass = zeros(Float64,tuple(size...));
		if UpperLeftCorner[1]>1
			mass[1,:].=1.0;
		end
		if UpperLeftCorner[2]>1
			mass[:,1].=1.0;
		end
		if BottomRightCorner[1] < nc[1]+1
			mass[end,:].=1.0;
		end
		if BottomRightCorner[2] < nc[2]+1
			mass[:,end].=1.0;
		end
	else
		error("Not fully implemented")
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1,(i[3]-1)*div(nc[3],NumCells[3]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1]),div(nc[2],NumCells[2]),div(nc[3],NumCells[3])];
		# # NODES X NODES X NODES
		# (UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [1;1;1],nc + [1;1;1],overlap);
	end
	return mass;
end

function getFacesStaggeredIndicesOfCell(NumCells::Array{Int64},overlap::Array{Int64},i::Array{Int64},nc::Array{Int64})
	indOver = [];
	if length(nc)==2
		nf = [prod(nc + [1; 0]),prod(nc + [0; 1])];

		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# NODES X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [1;0],nc + [1;0],overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);

		indNC = I1[:] .+ (I2[:] .- 1).*(nc[1]+1)

		# CELLS X NODES
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [0;1],nc + [0;1],overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCN = I1[:] .+ (I2[:].-1).*(nc[1])

		# CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner,nc,overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCC = I1[:] .+ (I2[:].-1).*(nc[1])

		indOver = [indNC; indCN .+ nf[1] ; indCC .+ (nf[1]+nf[2])];
	else
		nf = [prod(nc .+ [1; 0 ; 0]),prod(nc .+ [0; 1 ; 0]),prod(nc .+ [0; 0 ; 1])];
		# overlap = [4,4,4];
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1,(i[3]-1)*div(nc[3],NumCells[3]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1,div(nc[3],NumCells[3])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);

		# NODES X CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner + [1;0;0],nc + [1;0;0],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indNCC = I1[:] .+ (I2[:].-1).*(nc[1]+1) + (I3[:].-1).*((nc[1]+1)*nc[2]);

		# CELLS X NODES X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [0;1;0],nc + [0;1;0],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCNC = I1[:] .+ (I2[:].-1).*(nc[1]) + (I3[:].-1).*(nc[1]*(nc[2]+1));

		# CELLS X CELLS X NODES
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [0;0;1],nc + [0;0;1],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCCN = I1[:] .+ (I2[:].-1).*(nc[1]) + (I3[:].-1).*(nc[1]*nc[2]);

		# CELLS X CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner,nc,overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCCC = I1[:] .+ (I2[:].-1).*(nc[1]) + (I3[:].-1).*(nc[1]*nc[2]);

		indOver = [indNCC; indCNC .+ nf[1] ; indCCN .+ (nf[1]+nf[2]) ; indCCC .+ (nf[1]+nf[2]+nf[3])];
	end

	return indOver
end

function getFacesStaggeredIndicesOfCellNoPressure(NumCells::Array{Int64},overlap::Array{Int64},i::Array{Int64},nc::Array{Int64})
	indOver = [];
	if length(nc)==2
		nf = [prod(nc + [1; 0]),prod(nc + [0; 1])];

		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);

		# NODES X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [1;0],nc + [1;0],overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indNC = I1[:] .+ (I2[:].-1)*(nc[1]+1)

		# CELLS X NODES
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [0;1],nc + [0;1],overlap);
		I1,I2 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCN = I1[:] .+ (I2[:].-1)*(nc[1])

		indOver = [indNC; indCN .+ nf[1]];

	else
		nf = [prod(nc + [1; 0 ; 0]),prod(nc + [0; 1 ; 0]),prod(nc + [0; 0 ; 1])];
		# originalUpperLeftCorner = [(i[1]-1)*div(nc[1],NumCells[1]) + 1,(i[2]-1)*div(nc[2],NumCells[2]) + 1,(i[3]-1)*div(nc[3],NumCells[3]) + 1];
		# originalBottomRightCorner = originalUpperLeftCorner + [div(nc[1],NumCells[1])-1,div(nc[2],NumCells[2])-1,div(nc[3],NumCells[3])-1];
		(originalUpperLeftCorner,originalBottomRightCorner) = getOriginalBoundingBoxCells(NumCells,overlap,i,nc);
		# NODES X CELLS X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner + [1;0;0],nc + [1;0;0],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indNCC = I1[:] + (I2[:].-1)*(nc[1]+1) + (I3[:].-1)*((nc[1]+1)*nc[2]);

		# CELLS X NODES X CELLS
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [0;1;0],nc + [0;1;0],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCNC = I1[:] .+ (I2[:].-1)*(nc[1]) + (I3[:].-1)*(nc[1]*(nc[2]+1));

		# CELLS X CELLS X NODES
		(UpperLeftCorner,BottomRightCorner) = getBoxWithOverlap(originalUpperLeftCorner,originalBottomRightCorner+ [0;0;1],nc + [0;0;1],overlap);
		I1,I2,I3 = getBoxIdxs(UpperLeftCorner,BottomRightCorner);
		indCCN = I1[:] .+ (I2[:].-1)*(nc[1]) + (I3[:].-1)*(nc[1]*nc[2]);

		indOver = [indNCC; indCNC .+ nf[1] ; indCCN .+ (nf[1]+nf[2]) ];

	end
	return indOver
end


text:
####
#### equality operator for unit tests
####

≅(a, b) = a == b

≅(a::Pair, b::Pair) = a.first ≅ b.first && a.second ≅ b.second

≅(::Missing, ::Missing) = true

≅(a::AbstractVector, b::AbstractVector) = length(a) == length(b) && all(a .≅ b)

≅(a::Tuple, b::Tuple) = length(a) == length(b) && all(a .≅ b)

function ≅(a::FunctionalTable, b::FunctionalTable)
    keys(columns(a)) ≡ keys(columns(b)) &&
        ordering(a) == ordering(b) &&
        length(a) == length(b) &&
        all(((ac, bc), ) -> collect(ac) ≅ collect(bc),
            zip(values(columns(a)), values(columns(b))))
end

####
#### random vectors
####

function randvector_fs(; range = 1:1000)
    ranl() = rand() < 0.5 ? 1 : rand(range)
    [() -> fill(missing, ranl()),
     () -> fill(rand(Int8.(1:100)), ranl()),
     () -> fill(rand(Float64.(1:100)), ranl())]
end

function randvector(N, fs = randvector_fs())
    m = length(fs)
    v = Any[]
    for _ in 1:N
        v = vcat(v, fs[rand(1:m)]())
    end
    collect(v)
end

ext: jl


text:
function getAgesAtTime( mpSim::MPsim, timePoint::Float64,
    subpopulations::Vector{Subpopulation}, ageType::Symbol )

    # Get the ids in each subpopulation at the time point.
    result = Vector{Vector{Float64}}( undef, length( subpopulations ) )
    idsInSubpop = getSubpopulationAtTime( mpSim, timePoint, subpopulations )

    # ! No longer necessary after overhaul.
    if ageType === :timeInNode
        queryPrefix = string( "SELECT ", timePoint,
            " - max( timeIndex ) `time in node` FROM `", mpSim.transDBname,
            "` WHERE",
            "\n    timeIndex <= ", timePoint, " AND ",
            "\n    `", mpSim.idKey, "` IN ( '" )
        querySuffix = string( "' )",
            "\n    GROUP BY `", mpSim.idKey, "`" )

        for ii in eachindex( subpopulations )
            if isempty( idsInSubpop[ii] )
                result[ii] = Vector{Float64}()
            else
                queryCmd = string( queryPrefix, join( idsInSubpop[ii],
                    "', '" ), querySuffix )
                result[ii] = DataFrame( DBInterface.execute( mpSim.simDB,
                    queryCmd ) )[:, 1]
            end  # if isempty( idsInSubpop )[ii]
        end  # for ii in eachindex( subpopulations )

        return result
    end  # if ageType === :timeInNode

    queryCmdTmp = string( "SELECT ", timePoint, " - timeEntered",
        ageType === :age ? " + ageAtRecruitment age" : " tenure",
        " FROM `", mpSim.persDBname, "` WHERE ",
        "\n    `", mpSim.idKey, "` IN ( '" )

    for ii in eachindex( subpopulations )
        if isempty( idsInSubpop[ii] )
            result[ii] = Vector{Float64}()
        else
            queryCmd = string( queryCmdTmp, join( idsInSubpop[ii], "', '" ),
                "' )" )
            result[ii] = DataFrame( DBInterface.execute( mpSim.simDB,
                queryCmd ) )[:, 1]
        end  # if isempty( idsInSubpop )[ii]
    end  # for ii in eachindex( subpopulations )

    return result

end  # getAgesAtTime( mpSim, timePoint, subpopulationulations, ageType )

function processSubpopulationAges( personnelAges::Array{Vector{Float64},2},
    subpopulations::Vector{Subpopulation}, timeGrid::Vector{Float64},
    ageRes::Float64 )

    result = Dict{String, DataFrame}()
    baseNames = [:timePoint, :mean, :median, :stdev, :min, :max]

    for ii in eachindex( subpopulations )
        subpopName = subpopulations[ii].name
        agesOfSubpop = personnelAges[ii, :]
        ageSummary = zeros( Union{Float64, Missing}, length( timeGrid ), 6 )
        ageSummary[:, 1] = timeGrid
        emptySubpop = isempty.( agesOfSubpop )

        if all( emptySubpop )
            ageSummary[:, 2:end] .= missing
            result[subpopName] = DataFrame( ageSummary, baseNames )
        else
            minAge = minimum( minimum.( agesOfSubpop[.!emptySubpop] ) )
            minAge = floor( minAge / ageRes ) * ageRes
            maxAge = maximum( maximum.( agesOfSubpop[.!emptySubpop] ) )
            maxAge = floor( maxAge / ageRes ) * ageRes
            ageGrid = collect( minAge:ageRes:maxAge )

            counts = zeros( Int, length( timeGrid ), length( ageGrid ) )

            for jj in eachindex( timeGrid )
                subpopAgesAtTime = agesOfSubpop[jj]

                if isempty( subpopAgesAtTime )
                    ageSummary[jj, 2:6] .= missing
                else
                    ageSummary[jj, 2:6] = [mean( subpopAgesAtTime ),
                        median( subpopAgesAtTime ),
                        length( subpopAgesAtTime ) == 1 ? missing :
                            std( subpopAgesAtTime ),
                        minimum( subpopAgesAtTime ),
                        maximum( subpopAgesAtTime )]
                end  # if isempty( subpopAgesAtTime )

                tmpCounts = map( age -> count( subpopAgesAtTime .>= age ),
                    ageGrid )
                tmpCounts -= vcat( tmpCounts[2:end], 0 )
                counts[jj, :] = tmpCounts
            end  # for jj in eachindex( timeGrid )

            result[subpopName] = DataFrame( hcat( ageSummary, counts ),
                vcat( baseNames, Symbol.( ageGrid ) ) )
        end  # if all( emptySubpop )
    end  # for ii in eachindex( subpopulations )

    return result

end  # processSubpopulationAges( personnelAges,subpopulations, timeGrid,
     #   ageRes )


text:
using PyCall
using DataFrames
using SMACOF
using Plots
using Statistics
skl = pyimport("sklearn.manifold")

max_iter = 200
function skl_smacof(Δ; anchors = nothing)
    Yhat = skl.smacof(Δ, n_init = 1, eps = eps(), max_iter = max_iter)[1]
    return align(Yhat, anchors)
end

function main(nlst = [10,20, 30, 50, 100, 200, 300, 500], nrepeats = 10)
    first_run = true

    df = DataFrame(n = Int[], method = Symbol[], mse = Float64[], time = Float64[])
    for n in nlst
        for i in 1:nrepeats
            X = 2 * rand(n, 2)

            Δ = dist(X)
            t = time()
            Y_jl = smacof(Δ, anchors = X, ε = eps(), itmax = max_iter, Xinit = rand(n, 2))
            t_jl = time() - t
            # push!(df, [n, :julia, mse(Y, X), time() - t])

            t = time()
            Y_sk = skl_smacof(Δ, anchors = X)
            t_sk = time() - t
            # push!(df, [n, :sklearn, mse(Yhat, X), time() - t])

            if first_run
                first_run = false
                continue
            end
            push!(df, [n, :julia, mse(Y_jl, X), t_jl])
            push!(df, [n, :sklearn, mse(Y_sk, X), t_sk])
        end
    end

    df_mean_times = combine(groupby(df, [:method, :n]), :time => mean => :mean)
    sort!(df_mean_times, :n)

    df_julia_times = df[df[!, :method] .== :julia, [:n, :time]]
    plot(legend = :outertopleft)
    scatter!(df_julia_times[!, :n], df_julia_times[!, :time], label = "julia", color = 1)
    plot!(df_mean_times[df_mean_times[!, :method] .== :julia, :n], df_mean_times[df_mean_times[!, :method] .== :julia, :mean], color = 1, label = nothing)

    df_sklearn_times = df[df[!, :method] .== :sklearn, [:n, :time]]
    scatter!(df_sklearn_times[!, :n], df_sklearn_times[!, :time], label = "sklearn", color = 2)
    plot!(df_mean_times[df_mean_times[!, :method] .== :sklearn, :n], df_mean_times[df_mean_times[!, :method] .== :sklearn, :mean], color = 2, label = nothing)

    title!("First to $max_iter iterations")
    xlabel!("sample size (n)")
    ylabel!("time (sec)")
    savefig(joinpath(dirname(@__FILE__), "results.png"))

end

main()


text:
using Astellarn
using Test

# First test perft on known positions for common bugs
# https://www.chessprogramming.org/Perft_Results

@testset "rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - 0 1" begin
    fen = "rnbqkbnr/pppppppp/8/8/8/8/PPPPPPPP/RNBQKBNR w KQkq - 0 1"
    b = Astellarn.importfen(fen)
    @test cancastlekingside(b) == true
    @test cancastlekingside(b, BLACK) == true
    @test cancastlequeenside(b) == true
    @test cancastlequeenside(b, BLACK) == true
    @test ischeck(b) == false
    @test b.enpass == zero(UInt8)
    @test perft(b, 1) == 20
    @test perft(b, 2) == 400
    @test perft(b, 3) == 8902
    @test perft(b, 4) == 197281
    @test perft(b, 5) == 4865609
    @test perft(b, 6) == 119060324
end

@testset "r3k2r/p1ppqpb1/bn2pnp1/3PN3/1p2P3/2N2Q1p/PPPBBPPP/R3K2R w KQkq -" begin
    fen = "r3k2r/p1ppqpb1/bn2pnp1/3PN3/1p2P3/2N2Q1p/PPPBBPPP/R3K2R w KQkq -"
    b = importfen(fen)
    @test cancastlekingside(b) == true
    @test cancastlekingside(b, BLACK) == true
    @test cancastlequeenside(b) == true
    @test cancastlequeenside(b, BLACK) == true
    @test ischeck(b) == false
    @test b.enpass == zero(UInt8)
    @test perft(b, 1) == 48
    @test perft(b, 2) == 2039
    @test perft(b, 3) == 97862
    @test perft(b, 4) == 4085603
    @test perft(b, 5) == 193690690
    #@test perft(b, 6) == 8031647685
end

@testset "8/2p5/3p4/KP5r/1R3p1k/8/4P1P1/8 w - -" begin
    fen = "8/2p5/3p4/KP5r/1R3p1k/8/4P1P1/8 w - -"
    b = importfen(fen)
    @test cancastlekingside(b) == false
    @test cancastlekingside(b, BLACK) == false
    @test cancastlequeenside(b) == false
    @test cancastlequeenside(b, BLACK) == false
    @test ischeck(b) == false
    @test b.enpass == zero(UInt8)
    @test perft(b, 1) == 14
    @test perft(b, 2) == 191
    @test perft(b, 3) == 2812
    @test perft(b, 4) == 43238
    @test perft(b, 5) == 674624
    @test perft(b, 6) == 11030083
end

@testset "r3k2r/Pppp1ppp/1b3nbN/nP6/BBP1P3/q4N2/Pp1P2PP/R2Q1RK1 w kq - 0 1" begin
    fen1 = "r3k2r/Pppp1ppp/1b3nbN/nP6/BBP1P3/q4N2/Pp1P2PP/R2Q1RK1 w kq - 0 1"
    fen2 = "r2q1rk1/pP1p2pp/Q4n2/bbp1p3/Np6/1B3NBn/pPPP1PPP/R3K2R b KQ - 0 1"
    b1 = importfen(fen1)
    b2 = importfen(fen2)
    @test cancastlekingside(b1, BLACK) == true
    @test cancastlequeenside(b1, BLACK) == true
    @test cancastlekingside(b1, WHITE) == false
    @test cancastlequeenside(b1, WHITE) == false
    @test cancastlekingside(b2, WHITE) == true
    @test cancastlequeenside(b2, WHITE) == true
    @test cancastlekingside(b2, BLACK) == false
    @test cancastlequeenside(b2, BLACK) == false
    @test perft(b1, 1) == perft(b2, 1) == 6
    @test perft(b1, 2) == perft(b2, 2) == 264
    @test perft(b1, 3) == perft(b2, 3) == 9467
    @test perft(b1, 4) == perft(b2, 4) == 422333
    @test perft(b1, 5) == perft(b2, 5) == 15833292
end

@testset "SEE k7/8/6b1/3p4/4r3/3PQ3/5N2/K7 w - - 0 1" begin
    fen = "k7/8/6b1/3p4/4r3/3PQ3/5N2/K7 w - - 0 1"
    b = importfen(fen)
    t = 0
    move = Move(21, 28)
    @test static_exchange_evaluator(b, move, t) == true
    b.turn = !b.turn
    move = Move(28, 20)
    @test static_exchange_evaluator(b, move, t) == true
    move = Move(28, 36)
    @test static_exchange_evaluator(b, move, t) == false
end

@testset "SEE 7k/3r4/3q4/8/3p4/2BK4/8/8 w - - 0 1" begin
    fen = "7k/3r4/3q4/8/3p4/2BK4/8/8 w - - 0 1"
    b = importfen(fen)
    t = 0
    move = Move(22, 29)
    @test static_exchange_evaluator(b, move, t) == false
    fen = "3r3k/8/8/2p1p3/3p4/2BK4/2N5/8 w - - 0 1"
    b = importfen(fen)
    t = 0
    move = Move(22, 29)
    @test static_exchange_evaluator(b, move, t) == false
    fen = "3r3k/3r4/3q4/8/3p4/2BK4/1QN5/8 w - - 0 1"
    b = importfen(fen)
    t = 0
    move = Move(22, 29)
    @test static_exchange_evaluator(b, move, t) == true
    fen = "3r3k/8/8/2p1p3/3n4/2BKP3/8/8 w - - 0 1"
    b = importfen(fen)
    t = 0
    move = Move(22, 29)
    @test static_exchange_evaluator(b, move, t) == true
    fen = "k7/8/8/8/1rrpK3/3R4/8/8 w - - 0 1"
    b = importfen(fen)
    t = 0
    move = Move(21, 29)
    @test static_exchange_evaluator(b, move, t) == false
end

@testset "Wider SEE tests" begin
    fen_data = [
    ("4R3/2r3p1/5bk1/1p1r3p/p2PR1P1/P1BK1P2/1P6/8 b - -", Move(33, 26), 0),
    ("4R3/2r3p1/5bk1/1p1r1p1p/p2PR1P1/P1BK1P2/1P6/8 b - -", Move(33, 26), 0),
    ("4r1k1/5pp1/nbp4p/1p2p2q/1P2P1b1/1BP2N1P/1B2QPPK/3R4 b - -", Move(26, 19), 0),
    ("2r1r1k1/pp1bppbp/3p1np1/q3P3/2P2P2/1P2B3/P1N1B1PP/2RQ1RK1 b - -", Move(45, 36), SEE_VALUES[1]),
    ("7r/5qpk/p1Qp1b1p/3r3n/BB3p2/5p2/P1P2P2/4RK1R w - -", Move(4, 60), 0),
    ("6rr/6pk/p1Qp1b1p/2n5/1B3p2/5p2/P1P2P2/4RK1R w - -", Move(4, 60), -SEE_VALUES[4]),
    ("7r/5qpk/2Qp1b1p/1N1r3n/BB3p2/5p2/P1P2P2/4RK1R w - -", Move(4, 60), -SEE_VALUES[4]),
    ("6RR/4bP2/8/8/5r2/3K4/5p2/4k3 w - -", Move(51, 59, __QUEEN_PROMO), SEE_VALUES[3] - SEE_VALUES[1]),
    ("6RR/4bP2/8/8/5r2/3K4/5p2/4k3 w - -", Move(51, 59, __KNIGHT_PROMO), SEE_VALUES[2] - SEE_VALUES[1]),
    ("7R/5P2/8/8/8/3K2r1/5p2/4k3 w - -", Move(51, 59, __QUEEN_PROMO), SEE_VALUES[5] - SEE_VALUES[1]),
    ("7R/5P2/8/8/8/3K2r1/5p2/4k3 w - -", Move(51, 59, __BISHOP_PROMO), SEE_VALUES[3] - SEE_VALUES[1]),
    ("7R/4bP2/8/8/1q6/3K4/5p2/4k3 w - -", Move(51, 59, __ROOK_PROMO), -SEE_VALUES[1])
    ]
    for fd in fen_data
        @test static_exchange_evaluator(importfen(fd[1]), fd[2], fd[3]) == true
    end
end

@testset "Symmetry eval" begin
    fen = "2q5/b2k3P/4pn2/3pr3/3PR3/4PN2/B2K3p/2Q5 w - - 0 1"
    b = importfen(fen)
    pktable = PawnKingTable()
    eval_white = evaluate(b, pktable)
    b.turn = !b.turn
    eval_black = evaluate(b, pktable)
    @test abs(eval_black) == abs(eval_white)

    fen1 = "k7/P6n/KP5r/7b/7p/8/8/5RBN w - - 0 1"
    b = importfen(fen1)
    pktable = PawnKingTable()
    eval_white = evaluate(b, pktable)
    ei, ea = initEvalInfo(b)
    scorewk = evaluate_knights(b, ei, ea)
    scorewb = evaluate_bishops(b, ei, ea)
    scorewr = evaluate_rooks(b, ei, ea)
    scorewq = evaluate_queens(b, ei, ea)
    scorewki = evaluate_kings(b, ea)
    scorewpin = evaluate_pins(b)
    scorewspace = evaluate_space(b, ea)
    scorewthreat = evaluate_threats(b, ei, ea)
    scorewpass = evaluate_passed(b, ea)
    scorewinit = evaluate_initiative(b, ea.passed, 400)
    fen2 = "5rbn/8/8/7P/7B/kp5R/p6N/K7 b - - 0 1"
    b = importfen(fen2)
    eval_black = evaluate(b, pktable)
    ei, ea = initEvalInfo(b)
    scorebk = evaluate_knights(b, ei, ea)
    scorebb = evaluate_bishops(b, ei, ea)
    scorebr = evaluate_rooks(b, ei, ea)
    scorebq = evaluate_queens(b, ei, ea)
    scorebki = evaluate_kings(b, ea)
    scorebpin = evaluate_pins(b)
    scorebspace = evaluate_space(b, ea)
    scorebthreat = evaluate_threats(b, ei, ea)
    scorebpass = evaluate_passed(b, ea)
    scorebinit = evaluate_initiative(b, ea.passed, 400)

    @test scorebk == -scorewk
    @test scorebb == -scorewb
    @test scorebr == -scorewr
    @test scorebq == -scorewq
    @test scorebki == -scorewki
    @test scorebpin == -scorewpin
    @test scorebspace == -scorewspace
    @test scorebthreat == -scorewthreat
    @test scorebpass == -scorewpass
    @test scorebinit == -scorewinit
    @test eval_black == -eval_white
end

max_stars_repo_name: J-Revell/Astellarn

max_issues_repo_name: J-Revell/Astellarn

max_forks_repo_name: J-Revell/Astellarn


text:
"""
    ConstantAgent

Takes constant actions.
"""
struct ConstantAgent <: AbstractPolicy
    scale::Int
end
create_agent(::Val{ConstantAgent}; scale=50, kwargs...) = ConstantAgent(scale)

function (a::ConstantAgent)(env::AbstractEnv)
    fill(a.scale, length(action_space(env)))
end

"""
    RandomAgent

Takes random actions.
"""
struct RandomAgent <: AbstractPolicy
    rng::AbstractRNG
end
create_agent(::Val{:RandomAgent}; seed = 37, kwargs...) = RandomAgent(StableRNG(seed))

function (a::RandomAgent)(env::AbstractEnv)
    round.(Int, rand.(a.rng, action_space(env)))
end

"""
    SinusAgent

`period` is the period in steps, the whole curve will be between `0` and `2*amplitude`
"""
mutable struct SinusAgent <: AbstractPolicy
    period::Float64
    amplitude::Float64
    step::Int
end
create_agent(::Val{:SinusAgent}; period=1000, amplitude=40.0, kwargs...) = SinusAgent(period, amplitude, 0)

function (a::SinusAgent)(env::AbstractEnv)
    a.step += 1
    round.(Int, a.amplitude .* (1 .+ sin(a.step * 2π / a.period)) .+ ones(length(action_space(env))))
end

"""
    SimpleAgent

Something similar to the k8s scaler
"""
struct SimpleAgent <: AbstractPolicy
    down::Float64
    up::Float64
end
create_agent(::Val{:SimpleAgent}; down = 0.5, up = 0.8, kwargs...) = SimpleAgent(down, up)

function (agent::SimpleAgent)(env::AbstractEnv)
    env = env[!]
    a = Int[ms.running_nodes for ms in env.microservices]
    for (i, ms) in enumerate(env.microservices)
        util = length(ms.queue) / ms.running_nodes
        if util < agent.down
            a[i] -= 1
        elseif util > agent.up
            a[i] += 1
        end
    end
    lowas = [as.left for as in env.action_space]
    highas = [as.right for as in env.action_space]
    return clamp.(a, lowas, highas)
end

# Guided
function create_agent(::Val{:Guided};
    learner, teacher,
    kwargs...
)
    learner = create_agent(Val(learner); kwargs...)
    teacher = create_agent(Val(teacher); kwargs...)
    GuidedPolicy(learner, teacher; kwargs...)
end

"""
    K8sAgent

Something similar to the k8s scaler
"""
mutable struct K8sAgent <: AbstractPolicy
    desired_util::Float64
    step_tolerance::Float64
    scale_memory::Matrix{Int}
    mem_ptr::Int
end
create_agent(::Val{:K8sAgent}; desired_util=0.8, stabilization_steps=300, step_tolerance=0.1, env, kwargs...) = K8sAgent(desired_util, step_tolerance, zeros(Int, length(env[!].microservices), stabilization_steps), 1)

function (agent::K8sAgent)(env::AbstractEnv)
    env = env[!]

    newptr = mod(agent.mem_ptr, size(agent.scale_memory, 2)) + 1
    for (i, ms) in enumerate(env.microservices)
        rn = max(1e-6, ms.running_nodes) # To avoid problems when running_nodes == 0
        util = min(1, length(ms.queue) / rn)
        ratio = util / agent.desired_util
        if abs(ratio - 1) < agent.step_tolerance # No new action if within limit
            agent.scale_memory[i, newptr] = agent.scale_memory[i, agent.mem_ptr]
        else # Record new desired action
            agent.scale_memory[i, newptr] = ceil(Int, rn * ratio)
        end
    end

    agent.mem_ptr = newptr

    a = vec(maximum(agent.scale_memory, dims=2)) # Take max of desired action in window

    lowas = [as.left for as in env.action_space]
    highas = [as.right for as in env.action_space]
    return clamp.(a, lowas, highas)
end

"""
    Oracle

Optimal agent with assumption of knowledge of graph and delays, hardcoded for simpleflip env
"""
struct OracleAgent <: AbstractPolicy
end
create_agent(::Val{:OracleAgent}; kwargs...) = OracleAgent()

function (agent::OracleAgent)(env::AbstractEnv)
    env = env[!]
    a = zeros(Int, length(env.microservices))
    q = map(ms -> length(ms.queue), env.microservices)
    for tidx in 1:length(env.jobtypes)
        q = map(ms -> count(job -> job.typeindex == tidx, ms.queue), env.microservices[env.jobtypes[tidx].path])
        a[env.jobtypes[tidx].path] .+= [q[1] + 1; [max(q[i-1], q[i]) for i in 2:length(q)]]
    end
    lowas = [as.left for as in env.action_space]
    highas = [as.right for as in env.action_space]
    return clamp.(a, lowas, highas)
end

"""
    Loaded

Loads bson saved policy and acts on it using mean values of actions
"""
struct LoadedAgent{P<:GaussianNetwork} <: AbstractPolicy
    policy::P
end
function create_agent(::Val{:LoadedAgent}; loadpath::String, kwargs...)
    # Expects SAC, probably works with others with gaussian net
    # Though make sure the correct env is used then
    dict = BSON.load(loadpath)
    LoadedAgent(dict[:policy].policy.model)
end
function (agent::LoadedAgent)(env::AbstractEnv)
    s = state(env)
    action_mean, action_std = agent.policy(s)
    return action_mean
end

"""
    Wrapped

Wraps an env that returns from the inner action space, to a normalized action space.
"""
struct WrappedAgent{P<:AbstractPolicy,F<:Function} <: AbstractPolicy
    policy::P
    action_transform::F
end
function create_agent(::Val{:WrappedAgent}; env, wrapped_policy, kwargs...)
    A = action_space(env)
    alow = [x isa UnitRange ? x.start : x.left for x in A]
    ahigh = [x isa UnitRange ? x.stop : x.right for x in A]

    action_mapping(a) = 2 .* (a .- alow) ./ (ahigh .- alow) .- 1
    inner_agent = create_agent(Val(wrapped_policy); env=env[!], kwargs...)

    WrappedAgent(inner_agent, action_mapping)
end
function (agent::WrappedAgent)(env::AbstractEnv)
    a = agent.policy(env[!]) # Get action based on base env
    return agent.action_transform(a)
end

max_stars_repo_name: albheim/ServiceMeshControl

max_issues_repo_name: albheim/ServiceMeshControl

max_forks_repo_name: albheim/ServiceMeshControl


text:
mutable struct Genes
  IDs
  PPDE
  logFC
  weights
  colours

  function Genes(IDs, PPDE, logFC, weights, colours)
    new(IDs, PPDE, logFC, weights, colours)
  end

  function Genes(IDs)
    Genes(IDs, [], [], [], [])
  end
end

function nGenes(genes::Genes)
	return length(genes.IDs)
end

"""
Map DE data to model genes
"""
function mapDataToGenes!(genes::Genes, data::Data)
	inData = [findfirst(isequal(name), data.geneIDs) for name in genes.IDs]
	found = [name in data.geneIDs for name in genes.IDs]
	inData = inData[found]
	PPDE = zeros(nGenes(genes))
	logFC = zeros(nGenes(genes))
	PPDE[found] = data.PPDE[inData]
	logFC[found] = data.logFC[inData]

	genes.PPDE = PPDE
	genes.logFC = logFC
end

"""
Determine colours and weights for genes
"""
function colour!(genes::Genes, pThresh, alpha, precision)
	genes.colours = Int.(genes.PPDE .>0.9) .* sign.(genes.logFC)
	genes.weights = weight.(genes.PPDE, pThresh, alpha, precision)
end

max_stars_repo_name: htpusa/moominJL

max_issues_repo_name: htpusa/moominJL

max_forks_repo_name: htpusa/moominJL


text:
module Forms

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element

export form

register_normal_element("q__form", context = @__MODULE__)

function form(args...; noresetfocus::Bool = false, kwargs...)
  q__form(args...; kw(kwargs)...)
end

end


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

mutable struct VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties <: SwaggerModel
    idleTimeoutInMinutes::Any # spec type: Union{ Nothing, Int32 } # spec name: idleTimeoutInMinutes
    dnsSettings::Any # spec type: Union{ Nothing, VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings } # spec name: dnsSettings

    function VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties(;idleTimeoutInMinutes=nothing, dnsSettings=nothing)
        o = new()
        validate_property(VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        setfield!(o, Symbol("idleTimeoutInMinutes"), idleTimeoutInMinutes)
        validate_property(VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties, Symbol("dnsSettings"), dnsSettings)
        setfield!(o, Symbol("dnsSettings"), dnsSettings)
        o
    end
end # type VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties

const _property_map_VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties = Dict{Symbol,Symbol}(Symbol("idleTimeoutInMinutes")=>Symbol("idleTimeoutInMinutes"), Symbol("dnsSettings")=>Symbol("dnsSettings"))
const _property_types_VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties = Dict{Symbol,String}(Symbol("idleTimeoutInMinutes")=>"Int32", Symbol("dnsSettings")=>"VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings")
Base.propertynames(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties }) = collect(keys(_property_map_VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties))
Swagger.property_type(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties[name]))}
Swagger.field_name(::Type{ VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties }, property_name::Symbol) =  _property_map_VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties[property_name]

function check_required(o::VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties)
    true
end

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


text:
using Test
using CoffeeMachine

@testset "CoffeeMachine" begin
    @testset "Make a coffee" begin
        m = Machine(coffee=5.0)

        cup = makecoffee!(m)

        @test cup isa Cup
        @test cup.coffee > 0.0
    end

    @testset "A regular coffee does not have milk" begin
        m = Machine(coffee=5.0)

        cup = makecoffee!(m)

        @test cup.milk == 0.0
    end

    @testset "Out of coffee" begin
        m = Machine(coffee=0.0)

        cup = makecoffee!(m)

        @test cup === nothing
    end

    @testset "Out of coffee, displays message" begin
        m = Machine(coffee=0.0)

        cup = makecoffee!(m)

        @test cup === nothing
        @test readdisplay(m) == "Out of coffee"
    end

    @testset "Coffee with Milk" begin
        m = Machine(coffee=5.0, milk=5.0)

        cup = makecoffee!(m, withmilk=true)

        @test cup.milk > 0.0
    end

    @testset "Out of milk" begin
        m = Machine(coffee=5.0, milk=0.0)

        cup = makecoffee!(m, withmilk=true)

        @test cup === nothing
    end

    @testset "Out of milk displays message" begin
        m = Machine(coffee=5.0, milk=0.0)

        cup = makecoffee!(m, withmilk=true)

        @test readdisplay(m) == "Out of milk"
    end

    @testset "Fill milk" begin
        m = Machine(coffee=5.0, milk=0.0)

        cup = makecoffee!(m, withmilk=true)

        @test readdisplay(m) == "Out of milk"

        fillwith!(m, milk=5.0)

        cup = makecoffee!(m, withmilk=true)

        @test cup.milk > 0.0
        @test readdisplay(m) == "Enjoy"
    end

    @testset "Fill coffee" begin
        m = Machine(coffee=0.0)

        cup = makecoffee!(m)

        @test cup === nothing

        fillwith!(m, coffee=5.0)

        cup = makecoffee!(m)

        @test cup.coffee > 0.0
        @test readdisplay(m) == "Enjoy"
    end

    @testset "Running out of coffee" begin
        m = Machine(coffee=2.0)

        cup1 = makecoffee!(m)
        cup2 = makecoffee!(m)
        cup3 = makecoffee!(m)

        @test cup1 !== nothing
        @test cup2 !== nothing

        @test cup3 === nothing
    end

    @testset "Running out of milk" begin
        m = Machine(coffee=10.0, milk=2.0)

        cup1 = makecoffee!(m, withmilk=true)
        cup2 = makecoffee!(m, withmilk=true)
        cup3 = makecoffee!(m, withmilk=true)

        @test cup1 !== nothing
        @test cup2 !== nothing

        @test cup3 === nothing
        @test readdisplay(m) == "Out of milk"
    end
end


text:
module TestExecutors

include("preamble.jl")
import Distributed

@testset "fcollect" begin
    @test fcollect(x for x in 1:10) == 1:10
    @test fcollect(x for x in Iterators.Stateful(1:10)) == 1:10
    ch = Channel(ch -> foreach(x -> put!(ch, x), 1:10))
    @test fcollect(x for x in ch) == 1:10
    #=
    @sync let ch = Distributed.RemoteChannel()
        try
            @async foreach(x -> put!(ch, x), 1:10)
            @test fcollect(x for x in ch) == 1:10
        finally
            close(ch)
        end
    end
    =#
end

end  # module


text:
# parser

function read_instance(MyFileName::String)
  path = "./" * MyFileName
  # Si le fichier path existe
  if isfile(path)
    # L’ouvrir
    myFile = open(path)
    n = parse(Int64, split(readline(myFile), " ")[3])
    f = Array{Int64,1}(zeros(n-1)) # on met a zero pour l'instant
    c = Array{Int64,1}(zeros(n-1)) # on met a zero pour l'instant
    line1 =  split(readline(myFile), " ")
    line2 =  split(readline(myFile), " ")
    for i in 1:n-1
      f[i] = parse(Int64,line1[i+2])
      if parse(Int64,line2[i+2]) == 66
        c[i] = typemax(Int64)
      end
      if parse(Int64,line2[i+2]) != 66
        c[i] = parse(Int64,line2[i+2])
      end

    end
    w = parse(Int64,line1[n])
    d = parse(Int64, split(readline(myFile), " ")[3])
    return n, f, c, w, d
  end
end

max_stars_repo_name: boumpteryx/MLA

max_issues_repo_name: boumpteryx/MLA

max_forks_repo_name: boumpteryx/MLA


text:
@inline function production_eddy_viscosity(S::SymTen,W::AntiSymTen,c::Real,Δ2::Real)
    a = S:S
    nu = ifelse(a == 0.,zero(eltype(a)),
         @fastmath max(0.0,c*Δ2*((square(W):S) - (square(S):S))/(2a)))
    return nu
end


text:
@doc raw"""
    CanonicalMetric <: AbstractMetric

The Canonical Metric refers to a metric for the [`Stiefel`](@ref)
manifold, see[^EdelmanAriasSmith1998].

[^EdelmanAriasSmith1998]:
    > Edelman, A., Ariar, T. A., Smith, S. T.:
    > _The Geometry of Algorihthms with Orthogonality Constraints_,
    > SIAM Journal on Matrix Analysis and Applications (20(2), pp. 303–353, 1998.
    > doi: [10.1137/S0895479895290954](https://doi.org/10.1137/S0895479895290954)
    > arxiv: [9806030](https://arxiv.org/abs/physics/9806030)
"""
struct CanonicalMetric <: RiemannianMetric end

"""
    ApproximateLogarithmicMap <: ApproximateInverseRetraction

An approximate implementation of the logarithmic map, which is an [`inverse_retract`](@ref)ion.
See [`inverse_retract(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any, ::Any, ::ApproximateLogarithmicMap) where {n,k}`](@ref) for a use case.

# Fields

* `max_iterations` – maximal number of iterations used in the approximation
* `tolerance` – a tolerance used as a stopping criterion

"""
struct ApproximateLogarithmicMap{T} <: ApproximateInverseRetraction
    max_iterations::Int
    tolerance::T
end

function distance(M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p) where {n,k}
    return norm(M, p, log(M, p, q))
end

@doc raw"""
    q = exp(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, X)
    exp!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, q, CanonicalMetric}, p, X)

Compute the exponential map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref).

First, decompose The tangent vector ``X`` into its horizontal and vertical component with
respect to ``p``, i.e.

```math
X = pp^{\mathrm{T}}X + (I_n-pp^{\mathrm{T}})X,
```
where ``I_n`` is the ``n\times n`` identity matrix.
We introduce ``A=p^{\mathrm{T}}X`` and ``QR = (I_n-pp^{\mathrm{T}})X`` the `qr` decomposition
of the vertical component. Then using the matrix exponential ``\operatorname{Exp}`` we introduce ``B`` and ``C`` as

```math
\begin{pmatrix}
B\\C
\end{pmatrix}
\coloneqq
\operatorname{Exp}\left(
\begin{pmatrix}
A & -R^{\mathrm{T}}\\ R & 0
\end{pmatrix}
\right)
\begin{pmatrix}I_k\\0\end{pmatrix}
```

the exponential map reads

```math
q = \exp_p X = pC + QB.
```
For more details, see [^EdelmanAriasSmith1998][^Zimmermann2017].
"""
exp(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, ::Any...) where {n,k}

function exp!(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, q, p, X) where {n,k}
    A = p' * X
    n == k && return mul!(q, p, exp(A))
    QR = qr(X - p * A)
    BC_ext = exp([A -QR.R'; QR.R 0*I])
    @views begin
        mul!(q, p, BC_ext[1:k, 1:k])
        mul!(q, Matrix(QR.Q), BC_ext[(k + 1):(2 * k), 1:k], true, true)
    end
    return q
end

@doc raw"""
    inner(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, X, Y)

compute the inner procuct on the [`Stiefel`](@ref) manifold with respect to the
[`CanonicalMetric`](@ref). The formula reads

```math
g_p(X,Y) = \operatorname{tr}\bigl( X^{\mathrm{T}}(I_n - \frac{1}{2}pp^{\mathrm{T}})Y \bigr).
```
"""
function inner(::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric}, p, X, Y) where {n,k}
    T = Base.promote_eltype(p, X, Y)
    s = T(dot(X, Y))
    if n != k
        s -= dot(p'X, p'Y) / 2
    end
    return s
end

@doc raw"""
    X = inverse_retract(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap)
    inverse_retract!(M::MetricManifold{ℝ, Stiefel{n,k,ℝ}, X, CanonicalMetric}, p, q, a::ApproximateLogarithmicMap)

Compute an approximation to the logarithmic map on the [`Stiefel`](@ref)`(n,k)` manifold with respect to the [`CanonicalMetric`](@ref)
using a matrix-algebraic based approach to an iterative inversion of the formula of the
[`exp`](@ref exp(::MetricManifold{ℝ, Stiefel{n,k,ℝ}, CanonicalMetric}, ::Any...) where {n,k}).

The algorithm is derived in[^Zimmermann2017] and it uses the `max_iterations` and the `tolerance` field
from the [`ApproximateLogarithmicMap`](@ref).

[^Zimmermann2017]:
    > Zimmermann, R.: _A matrix-algebraic algorithm for the Riemannian logarithm on the Stiefel manifold under the canoncial metric.
    > SIAM Journal on Matrix Analysis and Applications 28(2), pp. 322-342, 2017.
    > doi: [10.1137/16M1074485](https://doi.org/10.1137/16M1074485),
    > arXiv: [1604.05054](https://arxiv.org/abs/1604.05054).
"""
inverse_retract(
    ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
    ::Any,
    ::Any,
    ::ApproximateLogarithmicMap,
) where {n,k}

function log(
    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
    p,
    q;
    maxiter::Int=10000,
    tolerance=1e-9,
) where {n,k}
    X = allocate_result(M, log, p, q)
    inverse_retract!(M, X, p, q, ApproximateLogarithmicMap(maxiter, tolerance))
    return X
end

function log!(
    M::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
    X,
    p,
    q;
    maxiter::Int=10000,
    tolerance=1e-9,
) where {n,k}
    inverse_retract!(M, X, p, q, ApproximateLogarithmicMap(maxiter, tolerance))
    return X
end

function inverse_retract!(
    ::MetricManifold{ℝ,Stiefel{n,k,ℝ},CanonicalMetric},
    X,
    p,
    q,
    a::ApproximateLogarithmicMap,
) where {n,k}
    M = p' * q
    QR = qr(q - p * M)
    V = convert(Matrix, qr([M; QR.R]).Q * I)
    Vcorner = @view V[(k + 1):(2 * k), (k + 1):(2k)] #bottom right corner
    Vpcols = @view V[:, (k + 1):(2 * k)] #second half of the columns
    S = svd(Vcorner) # preprocessing: Procrustes
    mul!(Vpcols, Vpcols * S.U, S.V')
    V[1:k, 1:k] .= M
    V[(k + 1):(2 * k), 1:k] .= QR.R

    LV = log_safe!(allocate(V), V)
    C = view(LV, (k + 1):(2 * k), (k + 1):(2 * k))
    expnC = exp(-C)
    i = 0
    new_Vpcols = Vpcols * expnC # allocate once
    while (i < a.max_iterations) && (norm(C) > a.tolerance)
        i = i + 1
        log_safe!(LV, V)
        expnC = exp(-C)
        mul!(new_Vpcols, Vpcols, expnC)
        copyto!(Vpcols, new_Vpcols)
    end
    mul!(X, p, @view(LV[1:k, 1:k]))
    # force the first - Q - to be the reduced form, not the full matrix
    mul!(X, @view(QR.Q[:, 1:k]), @view(LV[(k + 1):(2 * k), 1:k]), true, true)
    return X
end


text:
module TestSamples
    import Benchmarks
    using Base.Test
    using Compat

    s = Benchmarks.Samples()

    push!(s, 1, 1, 1, 1, 1)

    io = IOBuffer()
    show(io, s)

    path = tempname()
    writecsv(path, s)
    bytes = readstring(path)
    rm(path)

    empty!(s)
end


text:
"""
horner(c,x)

Evaluate a polynomial whose coefficients are given in ascending
order in `c`, at the point `x`, using Horner's rule.
"""
function horner(c,x)

n = length(c)
y = c[n]
for k in n-1:-1:1
    y = x*y + c[k]
end

return y
end

max_stars_repo_name: tobydriscoll/fnc

max_issues_repo_name: tobydriscoll/fnc

max_forks_repo_name: tobydriscoll/fnc


text:
#=
struct SAFTgammaMieParam <: EoSParam
    segment::SingleParam{Int}
    shapefactor::SingleParam{Float64}
    lambda_a::PairParam{Float64}
    lambda_r::PairParam{Float64}
    sigma::PairParam{Float64}
    epsilon::PairParam{Float64}
    epsilon_assoc::AssocParam{Float64}
    bondvol::AssocParam{Float64}
end

abstract type SAFTgammaMieModel <: GCSAFTModel end
@newmodelgc SAFTgammaMie SAFTgammaMieModel SAFTgammaMieParam

const SAFTγMie = SAFTgammaMie
export SAFTgammaMie,SAFTγMie
function SAFTgammaMie(components; idealmodel=BasicIdeal, userlocations=String[], ideal_userlocations=String[], verbose=false)
    groups = GroupParam(components, ["SAFT/SAFTgammaMie/SAFTgammaMie_groups.csv"]; verbose=verbose)
    params,sites = getparams(groups, ["SAFT/SAFTgammaMie"]; userlocations=userlocations, verbose=verbose)

    segment = params["vst"]
    shapefactor = params["S"]

    params["sigma"].values .*= 1E-10
    sigma = sigma_LorentzBerthelot(params["sigma"])
    epsilon = epsilon_HudsenMcCoubrey(params["epsilon"], sigma)
    lambda_a = lambda_LorentzBerthelot(params["lambda_a"])
    lambda_r = lambda_LorentzBerthelot(params["lambda_r"])
    epsilon_assoc = params["epsilon_assoc"]
    bondvol = params["bondvol"]

    #sites = SiteParam(Dict("e1" => params["n_e1"], "e2" => params["n_e2"], "H" => params["n_H"]))
    packagedparams = SAFTgammaMieParam(segment, shapefactor, lambda_a, lambda_r, sigma, epsilon, epsilon_assoc, bondvol)
    references = ["10.1063/1.4851455", "10.1021/je500248h"]

    model = SAFTgammaMie(packagedparams, groups, sites, idealmodel; ideal_userlocations=ideal_userlocations, references=references, verbose=verbose)
    return model
end

function a_res(model::SAFTgammaMieModel, V, T, z)
    return @f(a_mono) + @f(a_chain) + @f(a_assoc)
end

function a_mono(model::SAFTgammaMieModel, V, T, z)
    return @f(ÂHS) + @f(Â_1) + @f(Â_2) + @f(Â_3)
end

function a_disp(model::SAFTgammaMieModel, V, T, z)
    return @f(Â_1) + @f(Â_2) + @f(Â_3)
end

function a_chain(model::SAFTgammaMieModel, V, T, z)
    x = z/∑(z)
    v  = model.groups.n_flattenedgroups
    vst = model.params.segment.values
    S = model.params.shapefactor.values
    return -∑(x[i] * (∑(v[i][k]*vst[k]*S[k] for k ∈ @groups(i))-1) * log(@f(g_Mie,i)) for i ∈ @comps)
end

function a_assoc(model::SAFTgammaMieModel, V, T, z)
    x = z/∑(z)
    v = model.groups.n_flattenedgroups
    n = model.sites.n_sites
    X_ = @f(X)
    return ∑(x[i] * ∑(v[i][k] * ∑(n[k][a] * (log(X_[i][k][a])+(1-X_[i][k][a])/2) for a ∈ @sites(k)) for k ∈ @groups(i)) for i ∈ @comps)
end

function ÂHS(model::SAFTgammaMieModel, V, T, z)
    ζ_0   = @f(ζ,0)
    ζ_1   = @f(ζ,1)
    ζ_2   = @f(ζ,2)
    ζ_3   = @f(ζ,3)
    ρ = ∑(z)*N_A/V
    return 6/π/ρ * (3ζ_1*ζ_2/(1-ζ_3) + ζ_2^3/(ζ_3*(1-ζ_3)^2) + (ζ_2^3/ζ_3^2-ζ_0)*log(1-ζ_3))
end

function ζ0123(model::SAFTgammaMieModel, V, T, z)
    ζ_0   = @f(ζ,0)
    ζ_1   = @f(ζ,1)
    ζ_2   = @f(ζ,2)
    ζ_3   = @f(ζ,3)
    return ζ_0,ζ_1,ζ_2,ζ_3
end

function ∑Â_n(model::SAFTgammaMieModel, V, T, z)
    v = model.groups.n_flattenedgroups
    vst = model.params.segment.values
    S = model.params.shapefactor.values
    res = zero(V+T+first(z))
    @inbounds for i ∈ @comps
        res_i = zero(res)
        vi = v[i]
        groups_i = @groups(i)
        for idx in 1:length(groups_i)
            k = groups_i[idx]
            res_i += vi[k]*S[k]*vst[k]
        end
        res += z[i]*res_i
    end
    res1 = res/∑(z)
    #res2 = ∑(z[i]*∑(v[i][k]*vst[k]*S[k] for k ∈ @groups(i)) for i ∈ @comps)/∑(z)
    #       ∑(x[i]*∑(v[i][l]*vst[l]*S[l] for l ∈ @groups(i)) for i ∈ @comps)
    #@show res2
    return res1
end
function Â_1(model::SAFTgammaMieModel, V, T, z)
    return 1/T * @f(∑Â_n) * @f(a_1)
end
function Â_2(model::SAFTgammaMieModel, V, T, z)
    return 1/T^2 * @f(∑Â_n) * @f(a_2)
end
function Â_3(model::SAFTgammaMieModel, V, T, z)
    return 1/T^3 * @f(∑Â_n) * @f(a_3)
end

function a_1(model::SAFTgammaMieModel, V, T, z)
    groups = @groups
    return ∑(@f(x_S,k)*@f(x_S,l)*@f(a_1,k,l) for k ∈ groups for l ∈ groups)
end
function a_2(model::SAFTgammaMieModel, V, T, z)
    groups = @groups
    return ∑(@f(x_S,k)*@f(x_S,l)*@f(a_2,k,l) for k ∈ groups for l ∈ groups)
end
function a_3(model::SAFTgammaMieModel, V, T, z)
    groups = @groups
    return ∑(@f(x_S,k)*@f(x_S,l)*@f(a_3,k,l) for k ∈ groups for l ∈ groups)
end

function a_1(model::SAFTgammaMieModel, V, T, z, k, l)
    σ = model.params.sigma.values[k,l]
    λa = model.params.lambda_a.values[k,l]
    λr = model.params.lambda_r.values[k,l]

    x_0 = σ/@f(d,k,l)
    return @f(C,λa,λr) * (x_0^λa*(@f(aS_1,k,l,λa)+@f(B,k,l,λa)) - x_0^λr*(@f(aS_1,k,l,λr)+@f(B,k,l,λr)))
end
function a_2(model::SAFTgammaMieModel, V, T, z, k, l)
    σ = model.params.sigma.values[k,l]
    ϵ = model.params.epsilon.values[k,l]
    λa = model.params.lambda_a.values[k,l]
    λr = model.params.lambda_r.values[k,l]

    x_0 = σ/@f(d,k,l)
    return 1/2*@f(KHS)*(1+@f(χ,k,l))*ϵ*@f(C,λa,λr)^2 * (
           x_0^(2λa)*(@f(aS_1,k,l,2λa) + @f(B,k,l,2λa))
         - 2x_0^(λa+λr)*(@f(aS_1,k,l,λa+λr) + @f(B,k,l,λa+λr))
         + x_0^(2λr)*(@f(aS_1,k,l,2λr) + @f(B,k,l,2λr)) )
end
function a_3(model::SAFTgammaMieModel, V, T, z, k, l)
    ϵ = model.params.epsilon.values[k,l]
    ζst_X_ = @f(ζst_X)
    return -ϵ^3*@f(f,k,l,4)*ζst_X_ * exp(@f(f,k,l,5)*ζst_X_ + @f(f,k,l,6)*ζst_X_^2)
end

function B(model::SAFTgammaMieModel, V, T, z, k, l, λ)
    ϵ = model.params.epsilon.values[k,l]
    σ = model.params.sigma.values[k,l]
    x_0 = σ/@f(d,k,l)
    ζ_X_ = @f(ζ_X)
    I = (1-x_0^(3-λ))/(λ-3)
    J = (1-x_0^(4-λ)*(λ-3)+x_0^(3-λ)*(λ-4))/((λ-3)*(λ-4))
    return 2π*@f(ρ_S)*@f(d,k,l)^3*ϵ * ((1-ζ_X_/2)/(1-ζ_X_)^3*I-9ζ_X_*(1+ζ_X_)/(2(1-ζ_X_)^3)*J)
end

function aS_1(model::SAFTgammaMieModel, V, T, z, k, l, λ)
    ϵ = model.params.epsilon.values[k,l]
    ζeff_ = @f(ζeff, λ)
    return -2π*@f(ρ_S) * ϵ*@f(d,k,l)^3/(λ-3) * (1-ζeff_/2)/(1-ζeff_)^3
end

function ζ_X(model::SAFTgammaMieModel, V, T, z)
    groups = @groups
    return π*@f(ρ_S)/6 * ∑(@f(x_S,k)*@f(x_S,l)*@f(d,k,l)^3 for k ∈ groups for l ∈ groups)
end

function ζst_X(model::SAFTgammaMieModel, V, T, z)
    groups = @groups
    σ = model.params.sigma.values
    return π*@f(ρ_S)/6 * ∑(@f(x_S,k)*@f(x_S,l)*σ[k,l]^3 for k ∈ groups for l ∈ groups)
end

function ζeff(model::SAFTgammaMieModel, V, T, z, λ)
    A = SAFTγMieconsts.A
    ζ_X_ = @f(ζ_X)
    return A * SA[1; 1/λ; 1/λ^2; 1/λ^3] ⋅ SA[ζ_X_; ζ_X_^2; ζ_X_^3; ζ_X_^4]
end

function KHS(model::SAFTgammaMieModel, V, T, z)
    ζ_X_ = @f(ζ_X)
    return (1-ζ_X_)^4/(1+4ζ_X_+4ζ_X_^2-4ζ_X_^3+ζ_X_^4)
end

function χ(model::SAFTgammaMieModel, V, T, z, k, l)
    ζst_X_ = @f(ζst_X)
    return @f(f,k,l,1)*ζst_X_ + @f(f,k,l,2)*ζst_X_^5 + @f(f,k,l,3)*ζst_X_^8
end

function f(model::SAFTgammaMieModel, V, T, z, k, l, m)
    ϕ = SAFTγMieconsts.ϕ
    λa = model.params.lambda_a.values[k,l]
    λr = model.params.lambda_r.values[k,l]
    α = @f(C,λa,λr)*(1/(λa-3)-1/(λr-3))
    sum(ϕ[i+1][m]*α^i for i ∈ 0:3)/
    (1+sum(ϕ[i+1][m]*α^(i-3) for i ∈ 4:6))
end

function ζ(model::SAFTgammaMieModel, V, T, z, m)
    return π/6*@f(ρ_S)*∑(@f(x_S,k)*@f(d,k)^m for k ∈ @groups)
end

function ρ_S(model::SAFTgammaMieModel, V, T, z)
    #x = z/∑(z)
    #v = model.groups.n_flattenedgroups
    #vst = model.params.segment.values
    #S = model.params.shapefactor.values
    ρ = ∑(z)*N_A/V
    #res1 =  ρ * ∑(x[i] * ∑(v[i][k]*vst[k]*S[k] for k ∈ @groups(i)) for i ∈ @comps)
    res2 = ρ *@f(∑Â_n)
    return res2
end

function x_S(model::SAFTgammaMieModel, V, T, z, k)
    x = z/∑(z)
    v = model.groups.n_flattenedgroups
    vst = model.params.segment.values
    S = model.params.shapefactor.values
    return ∑(x[i]*v[i][k]*vst[k]*S[k] for i ∈ @comps) / ∑(x[i] * ∑(v[i][l]*vst[l]*S[l] for l ∈ @groups(i)) for i ∈ @comps)
end

function d(model::SAFTgammaMieModel, V, T, z, k)
    ϵ = model.params.epsilon.diagvalues[k]
    σ = model.params.sigma.diagvalues[k]
    λa = model.params.lambda_a.diagvalues[k]
    λr = model.params.lambda_r.diagvalues[k]
    u = SAFTγMieconsts.u
    w = SAFTγMieconsts.w
    θ = @f(C,λa,λr)*ϵ/T
    return σ*(1-∑(w[j]*(θ./(θ+u[j]))^(1/λr)*(exp(θ*(1/(θ./(θ+u[j]))^(λa/λr)-1))/(u[j]+θ)/λr) for j ∈ 1:5))
end

function d(model::SAFTgammaMieModel, V, T, z, k, l)
    # Lorentz mixing rule
    if k == l
        return @f(d, k)
    else
        return (@f(d, k) + @f(d, l))/2
    end
end

function C(model::SAFTgammaMieModel, V, T, z, λa, λr)
    return (λr/(λr-λa)) * (λr/λa)^(λa/(λr-λa))
end

function ẑ(model::SAFTgammaMieModel, V, T, z, i, k)
    v = model.groups.n_flattenedgroups
    vst = model.params.segment.values
    S = model.params.shapefactor.values
    return v[i][k]*vst[k]*S[k] / ∑(v[i][l]*vst[l]*S[l] for l ∈ @groups(i))
end

function σ̄(model::SAFTgammaMieModel, V, T, z, i)
    σ = model.params.sigma.values
    return cbrt(∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*σ[k,l]^3 for l ∈ @groups) for k ∈ @groups))
end

function σ̄(model::SAFTgammaMieModel, V, T, z, i, j)
    if i == j
        return @f(σ̄, i)
    else
        return (@f(σ̄, i) + @f(σ̄, j))/2
    end
end

function d̄(model::SAFTgammaMieModel, V, T, z, i)
    return cbrt(∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*@f(d,k,l)^3 for l ∈ @groups) for k ∈ @groups))
end

function ϵ̄(model::SAFTgammaMieModel, V, T, z, i)
    ϵ = model.params.epsilon.values
    return ∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*ϵ[k,l] for l ∈ @groups) for k ∈ @groups)
end

function ϵ̄(model::SAFTgammaMieModel, V, T, z, i, j)
    if i == j
        return @f(ϵ̄, i)
    else
        return sqrt(@f(σ̄,i)*@f(σ̄,j))/@f(σ̄,i,j) * sqrt(@f(ϵ̄,i)*@f(ϵ̄,j))
    end
end

function λ̄a(model::SAFTgammaMieModel, V, T, z, i)
    λa = model.params.lambda_a.values
    return ∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*λa[k,l] for l ∈ @groups) for k ∈ @groups)
end

function λ̄r(model::SAFTgammaMieModel, V, T, z, i)
    λr = model.params.lambda_r.values
    return ∑(∑(@f(ẑ,i,k)*@f(ẑ,i,l)*λr[k,l] for l ∈ @groups) for k ∈ @groups)
end

function g_Mie(model::SAFTgammaMieModel, V, T, z, i)
    ϵ̄_=@f(ϵ̄,i)
    g_HSi = @f(g_HS,i)
    return g_HSi*exp(ϵ̄_/T*@f(g_1,i)/g_HSi+(ϵ̄_/T)^2*@f(g_2,i)/g_HSi)
end

function g_HS(model::SAFTgammaMieModel, V, T, z, i)
    x̄_0 = @f(σ̄,i)/@f(d̄,i)
    ζ_X_ = @f(ζ_X)
    k_0 = -log(1-ζ_X_) + (42ζ_X_-39ζ_X_^2+9ζ_X_^3-2ζ_X_^4)/(6*(1-ζ_X_)^3)
    k_1 = (ζ_X_^4+6ζ_X_^2-12ζ_X_)/(2*(1-ζ_X_)^3)
    k_2 = -3ζ_X_^2/(8(1-ζ_X_)^2)
    k_3 = (-ζ_X_^4+3ζ_X_^2+3ζ_X_)/(6(1-ζ_X_)^3)
    return exp(k_0+k_1*x̄_0+k_2*x̄_0^2+k_3*x̄_0^3)
end

function g_1(model::SAFTgammaMieModel, V, T, z, i)
    λ̄a_ = @f(λ̄a,i)
    λ̄r_ = @f(λ̄r,i)
    σ̄_ = @f(σ̄,i)
    ϵ̄_ = @f(ϵ̄,i)
    d̄_ = @f(d̄,i)
    ρ_S_ = @f(ρ_S)
    x̄_0 = σ̄_/@f(d̄,i)
    C̄_ = @f(C,λ̄a_,λ̄r_)
    return 1/(2π*ϵ̄_*d̄_^3*ρ_S_)*(3*@f(∂ā_1╱∂ρ_S,i) - C̄_*λ̄a_*x̄_0^λ̄a_*(@f(āS_1,i,λ̄a_)+@f(B̄,i,λ̄a_)) + C̄_*λ̄r_*x̄_0^λ̄r_*(@f(āS_1,i,λ̄r_)+@f(B̄,i,λ̄r_)))
end

function g_2(model::SAFTgammaMieModel, V, T, z, i)
    return (1+@f(γ_c,i))*@f(gMCA_2,i)
end

function γ_c(model::SAFTgammaMieModel, V, T, z, i)
    λ̄a_ = @f(λ̄a,i)
    λ̄r_ = @f(λ̄r,i)
    ϵ̄_ = @f(ϵ̄,i)
    ᾱ = @f(C,λ̄a_,λ̄r_)*(1/(λ̄a_-3)-1/(λ̄r_-3))
    θ = exp(ϵ̄_/T)-1
    ζst_X_ = @f(ζst_X)
    return 10 * (-tanh(10*(0.57-ᾱ))+1) * ζst_X_*θ*exp(-6.7*ζst_X_-8ζst_X_^2)
end

function gMCA_2(model::SAFTgammaMieModel, V, T, z, i)
    ϵ̄_ = @f(ϵ̄,i)
    σ̄_ = @f(σ̄,i)
    d̄_ = @f(d̄,i)
    x̄_0 = σ̄_/d̄_
    λ̄a_ = @f(λ̄a,i)
    λ̄r_ = @f(λ̄r,i)
    ρ_S_ = @f(ρ_S)
    C̄_ = @f(C,λ̄a_,λ̄r_)
    KHS_ = @f(KHS)
    return 1/(2π*ϵ̄_^2*d̄_^3*ρ_S_)*(3*@f(∂ā_2╱∂ρ_S,i)-
        ϵ̄_*KHS_*C̄_^2*λ̄r_*x̄_0^(2λ̄r_)*(@f(āS_1,i,2λ̄r_)+@f(B̄,i,2λ̄r_))+
        ϵ̄_*KHS_*C̄_^2*(λ̄a_+λ̄r_)*x̄_0^(λ̄a_+λ̄r_)*(@f(āS_1,i,λ̄a_+λ̄r_)+@f(B̄,i,λ̄a_+λ̄r_))-
        ϵ̄_*KHS_*C̄_^2*λ̄a_*x̄_0^(2λ̄a_)*(@f(āS_1,i,2λ̄a_)+@f(B̄,i,2λ̄a_)))
end

function B̄(model::SAFTgammaMieModel, V, T, z, i, λ̄)
    ϵ̄_ = @f(ϵ̄,i)
    σ̄_ = @f(σ̄,i)
    d̄_ = @f(d̄,i)
    x̄_0 = σ̄_/d̄_
    ζ_X_ = @f(ζ_X)
    I = (1-x̄_0^(3-λ̄))/(λ̄-3)
    J = (1-x̄_0^(4-λ̄)*(λ̄-3)+x̄_0^(3-λ̄)*(λ̄-4))/((λ̄-3)*(λ̄-4))
    return 2π*@f(ρ_S)*d̄_^3*ϵ̄_ * ((1-ζ_X_/2)/(1-ζ_X_)^3*I-9ζ_X_*(1+ζ_X_)/(2(1-ζ_X_)^3)*J)
end

function ∂B╱∂ρ_S(model::SAFTgammaMieModel, V, T, z, i, λ̄)
    ϵ̄_ = @f(ϵ̄,i)
    σ̄_ = @f(σ̄,i)
    d̄_ = @f(d̄,i)
    x̄_0 = σ̄_/d̄_
    ζ_X_ = @f(ζ_X)
    I = (1-x̄_0^(3-λ̄))/(λ̄-3)
    J = (1-x̄_0^(4-λ̄)*(λ̄-3)+x̄_0^(3-λ̄)*(λ̄-4))/((λ̄-3)*(λ̄-4))
    return @f(B̄,i,λ̄) + 2π*@f(ρ_S)*d̄_^3*ϵ̄_ * ζ_X_*((3*(1-ζ_X_/2)*(1-ζ_X_)^2
            -0.5*(1-ζ_X_)^3)*I/(1-ζ_X_)^6-9*J*((1+2*ζ_X_)*(1-ζ_X_)^3+ζ_X_*(1+ζ_X_)*3*(1-ζ_X_)^2)/(2*(1-ζ_X_)^6))
end

function ā_1(model::SAFTgammaMieModel, V, T, z, i)
    λ̄a_ = @f(λ̄a,i)
    λ̄r_ = @f(λ̄r,i)
    x̄_0 = @f(σ̄,i)/@f(d̄,i)
    return @f(C,λ̄a_,λ̄r_)*(x̄_0^λ̄a_*(@f(āS_1,i,λ̄a_)+@f(B̄,i,λ̄a_))-x̄_0^λ̄r_*(@f(āS_1,i,λ̄r_)+@f(B̄,i,λ̄r_)))
end

function ∂ā_1╱∂ρ_S(model::SAFTgammaMieModel, V, T, z, i)
    λ̄a_ = @f(λ̄a,i)
    λ̄r_ = @f(λ̄r,i)
    x̄_0 = @f(σ̄,i)/@f(d̄,i)
    return @f(C,λ̄a_,λ̄r_)*(x̄_0^λ̄a_*(@f(∂āS_1╱∂ρ_S,i,λ̄a_)+@f(∂B╱∂ρ_S,i,λ̄a_))
                          - x̄_0^λ̄r_*(@f(∂āS_1╱∂ρ_S,i,λ̄r_)+@f(∂B╱∂ρ_S,i,λ̄r_)))
end

function āS_1(model::SAFTgammaMieModel, V, T, z, i, λ̄)
    ϵ̄_ = @f(ϵ̄,i)
    σ̄_ = @f(σ̄,i)
    d̄_ = @f(d̄,i)
    ζ̄eff_ = @f(ζeff, λ̄)
    return -2π*@f(ρ_S) * ϵ̄_*d̄_^3/(λ̄-3) * (1-ζ̄eff_/2)/(1-ζ̄eff_)^3
end

function ∂āS_1╱∂ρ_S(model::SAFTgammaMieModel, V, T, z, i, λ̄)
    ϵ̄_ = @f(ϵ̄,i)
    σ̄_ = @f(σ̄,i)
    d̄_ = @f(d̄,i)
    ζ̄eff_ = @f(ζeff, λ̄)
    A = SAFTγMieconsts.A
    ζ_X_ = @f(ζ_X)
    ∂ζ̄eff╱∂ρ_S = A * SA[1; 1/λ̄; 1/λ̄^2; 1/λ̄^3] ⋅ SA[1; 2ζ_X_; 3ζ_X_^2; 4ζ_X_^3]
    return @f(āS_1,i,λ̄) - 2π*(ϵ̄_*d̄_^3)/(λ̄-3) *@f(ρ_S)* ((3*(1-ζ̄eff_/2)*(1-ζ̄eff_)^2-1/2*(1-ζ̄eff_)^3)/(1-ζ̄eff_)^6 * ∂ζ̄eff╱∂ρ_S*ζ_X_)
end

function ∂ā_2╱∂ρ_S(model::SAFTgammaMieModel, V, T, z, i)
    ϵ̄_ = @f(ϵ̄,i)
    σ̄_ = @f(σ̄,i)
    d̄_ = @f(d̄,i)
    x̄_0 = σ̄_/d̄_
    ζ_X_ = @f(ζ_X)
    ρ_S_ = @f(ρ_S)
    ∂KHS╱∂ρ_S = -(4*(1-ζ_X_)^3*(1+4ζ_X_+4ζ_X_^2-4ζ_X_^3+ζ_X_^4)+(1-ζ_X_)^4*(4+8ζ_X_-12ζ_X_^2+4ζ_X_^3))/(1+4ζ_X_+4ζ_X_^2-4ζ_X_^3+ζ_X_^4)^2*ζ_X_
    λ̄a_ = @f(λ̄a,i)
    λ̄r_ = @f(λ̄r,i)
    return 1/2*ϵ̄_*@f(C,λ̄a_,λ̄r_)^2*(∂KHS╱∂ρ_S*(x̄_0^(2λ̄a_)*(@f(āS_1,i,2λ̄a_)+@f(B̄,i,2λ̄a_))
        -2x̄_0^(λ̄a_+λ̄r_)*(@f(āS_1,i,λ̄a_+λ̄r_)+@f(B̄,i,λ̄a_+λ̄r_))
        +x̄_0^(2λ̄r_)*(@f(āS_1,i,2λ̄r_)+@f(B̄,i,2λ̄r_)))
        +@f(KHS)*(x̄_0^(2λ̄a_)*(@f(∂āS_1╱∂ρ_S,i,2λ̄a_)+@f(∂B╱∂ρ_S,i,2λ̄a_))
        -2x̄_0^(λ̄a_+λ̄r_)*(@f(∂āS_1╱∂ρ_S,i,λ̄a_+λ̄r_)+@f(∂B╱∂ρ_S,i,λ̄a_+λ̄r_))
        +x̄_0^(2λ̄r_)*(@f(∂āS_1╱∂ρ_S,i,2λ̄r_)+@f(∂B╱∂ρ_S,i,2λ̄r_))))
end

function X(model::SAFTgammaMieModel, V, T, z)
    _1 = one(V+T+first(z))
    _0 = zero(_1)
    x = z/∑(z)
    ρ = ∑(z)*N_A/V
    v = model.groups.n_flattenedgroups
    n = model.sites.n_sites
    itermax = 1000
    damping_factor = 0.5
    error = 1.
    error_old = 1.
    ndamping = 2
    tol = model.absolutetolerance
    iter = 1
    X_ = [[[_1 for a ∈ @sites(k)] for k ∈ @groups] for i ∈ @comps]
    X_old = deepcopy(X_)
    while error > tol
        if iter > itermax
            throw("X has failed to converge after $itermax iterations")
        end
        @inbounds for i ∈ @comps, k ∈ @groups(i), a ∈ @sites(k)
            #rhs = 1/(1+ρ*∑(x[j] * ∑(v[j][l] * ∑(n[l][b] * X_old[j][l][b] * @f(Δ,i,j,k,l,a,b) for b ∈ @sites(l),init=_0) for l ∈ @groups(j),init=_0) for j ∈ @comps,init=_0))
            res = _0
            for j ∈ @comps
                res_j = _0
                for l ∈ @groups(j)
                    res_l = _0
                    for b ∈ @sites(l)
                        res_l += n[l][b] * X_old[j][l][b] * @f(Δ,i,j,k,l,a,b)
                    end
                    res_j += v[j][l]*res_l
                end
                res += x[j]*res_j
            end

            rhs = (1/(1+ρ*res))
            @inbounds begin
            X_[i][k][a] = (1-damping_factor)*X_old[i][k][a] + damping_factor*rhs
            end
        end
        #error = sqrt(∑(∑(∑((X_[i][k][a]-X_old[i][k][a])^2 for a ∈ @sites(k)) for k ∈ @groups(i)) for i ∈ @comps))
        _error = _0
        for i = 1:length(X_)
            for j in 1:length(X_[i])
                for k in 1:length(X_[i][j])
                    new = X_[i][j][k]
                    old = X_old[i][j][k]
                    _error += abs2(old-new)
                   X_old[i][j][k] = new
                end
            end
        end
        error_old = error
        error = sqrt(_error)
        if error < 1e-4
            damping_factor = 1.0
        end
        iter += 1
    end

    return X_
end

function Δ(model::SAFTgammaMieModel, V, T, z, i, j, k, l, a, b)
    σ = model.params.sigma.values
    σ3_x = ∑(∑(@f(x_S,k)*@f(x_S,l)*σ[k,l]^3 for k ∈ @groups) for l ∈ @groups)
    ϵ̄_ = @f(ϵ̄,i,j)
    ρR = @f(ρ_S)*σ3_x
    TR = T/ϵ̄_
    c  = SAFTγMieconsts.c
    I = ∑(∑(c[p+1,q+1]*ρR^p*TR^q for q ∈ 0:(10-p)) for p ∈ 0:10)

    ϵHB = model.params.epsilon_assoc.values
    K = model.params.bondvol.values
    F = (exp(ϵHB[k,l][a,b]/T)-1)
    return F*K[k,l][a,b]*I
end

const SAFTγMieconsts =(
    A = SA[ 0.81096    1.7888   -37.578   92.284;
          1.02050  -19.341    151.26  -463.50 ;
         -1.90570   22.845   -228.14   973.92 ;
          1.08850   -6.1962   106.98  -677.64  ],

    ϕ = [[7.5365557, -359.440,  1550.9, -1.199320, -1911.2800,  9236.9  ],
        [-37.604630,  1825.60, -5070.1,  9.063632,  21390.175, -129430 ],
        [71.745953, -3168.00,  6534.6, -17.94820, -51320.700,  357230  ],
        [-46.835520,  1884.20, -3288.7,  11.34027,  37064.540, -315530 ],
        [-2.4679820,- 0.82376, -2.7171,  20.52142,  1103.7420,  1390.2 ],
        [-0.5027200, -3.19350,  2.0883, -56.63770, -3264.6100, -4518.2 ],
        [8.0956883,  3.70900,  0.0000,  40.53683,  2556.1810,  4241.6  ]],

    u = [0.26356031971814109102031,1.41340305910651679221800,3.59642577104072208122300,7.08581000585883755692200,12.6408008442757826594300],
    w = [0.5217556105828086524759,0.3986668110831759274500,7.5942449681707595390e-2,3.6117586799220484545e-3,2.3369972385776227891e-5],

         c  = [0.0756425183020431	-0.128667137050961	 0.128350632316055	-0.0725321780970292	   0.0257782547511452  -0.00601170055221687	  0.000933363147191978  -9.55607377143667e-05  6.19576039900837e-06 -2.30466608213628e-07 3.74605718435540e-09
                0.134228218276565	    -0.182682168504886 	 0.0771662412959262	-0.000717458641164565 -0.00872427344283170	0.00297971836051287	 -0.000484863997651451	 4.35262491516424e-05 -2.07789181640066e-06	4.13749349344802e-08 0
               -0.565116428942893	     1.00930692226792   -0.660166945915607	 0.214492212294301	  -0.0388462990166792	0.00406016982985030	 -0.000239515566373142	 7.25488368831468e-06 -8.58904640281928e-08	0	                 0
               -0.387336382687019	    -0.211614570109503	 0.450442894490509	-0.176931752538907	   0.0317171522104923  -0.00291368915845693	  0.000130193710011706  -2.14505500786531e-06  0	                0	                 0
                2.13713180911797	    -2.02798460133021 	 0.336709255682693	 0.00118106507393722  -0.00600058423301506	0.000626343952584415 -2.03636395699819e-05	 0	                   0	                0	                 0
               -0.300527494795524	     2.89920714512243   -0.567134839686498	 0.0518085125423494	  -0.00239326776760414	4.15107362643844e-05  0	                     0	                   0	                0                    0
               -6.21028065719194	    -1.92883360342573	 0.284109761066570	-0.0157606767372364	   0.000368599073256615	0 	                  0	                     0	                   0	                0	                 0
                11.6083532818029	     0.742215544511197  -0.0823976531246117	 0.00186167650098254   0	                0	                  0	                     0	                   0	                0	                 0
               -10.2632535542427	    -0.125035689035085	 0.0114299144831867	 0	                   0	                0	                  0	                     0	                   0	                0	                 0
                4.65297446837297	    -0.00192518067137033 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0
               -0.867296219639940	     0	                 0	                 0	                   0	                0	                  0	                     0	                   0	                0	                 0],
)
=#


text:
@static if Sys.isapple()
    struct Cpasswd
        pw_name::Cstring
        pw_passwd::Cstring
        pw_uid::Cint
        pw_gid::Cint
        pw_change::Cint
        pw_class::Cstring
        pw_gecos::Cstring
        pw_dir::Cstring
        pw_shell::Cstring
        pw_expire::Cint
        pw_fields::Cint
    end
elseif Sys.islinux()
    struct Cpasswd
       pw_name::Cstring
       pw_passwd::Cstring
       pw_uid::Cint
       pw_gid::Cint
       pw_gecos::Cstring
       pw_dir::Cstring
       pw_shell::Cstring
    end
else
    struct Cpasswd
        pw_name::Cstring
        pw_uid::Cint
        pw_gid::Cint
        pw_dir::Cstring
        pw_shell::Cstring
    end

    Cpasswd() = Cpasswd(pointer("NA"), 0, 0, pointer("NA"), pointer("NA"))
end

struct Cgroup
    gr_name::Cstring
    gr_passwd::Cstring
    gr_gid::Cint
end

Cgroup() = Cgroup(pointer("NA"), pointer("NA"), 0)

struct User
    name::String
    uid::UInt64
    gid::UInt64
    dir::String
    shell::String
end

function User(ps::Cpasswd)
    User(
        unsafe_string(ps.pw_name),
        UInt64(ps.pw_uid),
        UInt64(ps.pw_gid),
        unsafe_string(ps.pw_dir),
        unsafe_string(ps.pw_shell)
    )
end

User(passwd::Ptr{Cpasswd}) = User(unsafe_load(passwd))

function Base.show(io::IO, user::User)
    print(io, "$(user.uid) ($(user.name))")
end

function User(name::String)
    ps = @static if Sys.isunix()
        ccall(:getpwnam, Ptr{Cpasswd}, (Ptr{UInt8},), name)
    else
        Cpasswd()
    end

    User(ps)
end

function User(uid::UInt)
    ps = @static if Sys.isunix()
        ccall(:getpwuid, Ptr{Cpasswd}, (UInt64,), uid)
    else
        Cpasswd()
    end

    User(ps)
end

function User()
    uid = @static Sys.isunix() ? ccall(:geteuid, Cint, ()) : 0
    User(UInt64(uid))
end

struct Group
    name::String
    gid::UInt64
end

Group(gr::Cgroup) = Group(unsafe_string(gr.gr_name), UInt64(gr.gr_gid))
Group(group::Ptr{Cgroup}) = Group(unsafe_load(group))

function Base.show(io::IO, group::Group)
    print(io, "$(group.gid) ($(group.name))")
end

function Group(name::String)
    ps = @static if Sys.isunix()
        ccall(:getgrnam, Ptr{Cgroup}, (Ptr{UInt8},), name)
    else
        Cgroup()
    end

    Group(ps)
end

function Group(gid::UInt)
    gr = @static if Sys.isunix()
        ccall(:getgrgid, Ptr{Cgroup}, (UInt64,), gid)
    else
        Cgroup()
    end

    Group(gr)
end

function Group()
    gid = @static Sys.isunix() ? ccall(:getegid, Cint, ()) : 0
    Group(UInt64(gid))
end


text:
struct DenseBatchedCache{Ta,Tb,Tc,Td,Te}
    σ::Tc
    out::Tb
    δℒℒ::Td
    δℒ::Te

    θ::Ta
    out2::Tb
    valid::Bool
end

function cache(l::Dense{Ta,Tb}, arr_T, in_T ,in_sz, batch_sz) where {Ta,Tb}
    in_T2 = promote_type(in_T, eltype(l.W))

    c = DenseBatchedCache(similar(l.W, in_T2, size(l.W,2), batch_sz),
                          similar(l.b, in_T2, size(l.b,1), batch_sz),
                          similar(l.W, in_T2, size(l.W,1), batch_sz),
                          similar(l.W, in_T2, size(l.W,2), batch_sz),

                          similar(l.b, in_T2, size(l.b,1), batch_sz),
                          similar(l.b, in_T2, size(l.b,1), batch_sz),
                          false)
    return c
end

batch_size(c::DenseBatchedCache) = size(c.out, 2)

function (l::Dense)(c::DenseBatchedCache, x)
    # The preallocated caches
    out  = c.out
    θ = c.θ

    # Store the input to this layer for the backpropagation
    σ = copyto!(c.σ, x)

    #θ .= net.b .+ net.W * x
    mul!(θ, l.W, σ)
    θ .+= l.b

    # Apply the nonlinear function
    out  .= l.σ.(θ)
    return out
end

function backprop(∇, l::Dense, c::DenseBatchedCache, δℒ)
    # The preallocated caches
    θ = c.θ
    δℒℒ = c.δℒℒ

    # Compute the actual sensitivity
    copyto!(δℒℒ, δℒ)
    δℒℒ .*= fwd_der.(l.σ, θ)

    _batched_outer_prod_noconj!(∇.W, δℒℒ, c.σ)
    ∇.b .= δℒℒ

    δℒ = mul!(c.δℒ, transpose(l.W), δℒℒ)
    return δℒ
end


text:
## define initial values ##

function defineY(N, y0, β)

    Y = zeros(N + 1)
    Y[1] = y0
    Y

end

function defineY(N, y0::Vector{<:Real}, β)

    if size(y0) == size(β)

        Y = zeros(N + 1, length(β))
        Y[1, :] .= y0
        Y

    elseif size(y0) != size(β)

        Y = zeros(N + 1)
        Y[1] = y0[1]
        Y

    end

end

function defineY(N, y0::Matrix{<:Real}, β)

    Y = zeros(N + 1, size(y0)[2])
    Y[1, :] .= y0[1, :]
    Y

end

## indexY ##
function indexY(n, Y::Vector{<:Real}, Ynext)

    Y[n + 1] = Ynext
    Y

end

function indexY(n, Y::Matrix{<:Real}, Ynext)

    Y[n + 1, :] = Ynext
    Y

end

## taylor_expansion ##
taylor_expansion(t0, t, y0, β, m) = y0

function taylor_expansion(t0, t, y0::Vector{<:Real}, β, m)

    if size(y0) == size(β)

        return y0

    elseif size(y0) != size(β)

        return sum([(t - t0) ^ k / factorial(k) * y0[k + 1] for k in 0:m - 1])

    end

end

taylor_expansion(t0, t, y0::Matrix{<:Real}, β, m) = sum([(t - t0) ^ k / factorial(k) * y0[:, k + 1] for k in 0:m - 1])

## function ϕ ##
function Phi(y, F, β, t, n, par...)

    alpha = zeros(n - 1)
    fun = zeros(n - 1)

    alpha = map((x) -> α(x, β), n - 1:-1:1)
    fun = map((x) -> F(t, x, β, y, par...), 2:n)

    a_n0(n, β) .* F(t, 1, β, y, par...) .+ sum([alpha[i] .* fun[i] for i in 1:n - 1])

end

a_n0(n, β) = ((n - 1) .^ (β .+ 1) .- n .^ β .* (n .- β .- 1)) ./ Γ(β .+ 2)

## function α ##
function α(n, β)

    if n == 0

        return 1 ./ Γ(β .+ 2)

    else

        return ((n .- 1) .^ (β .+ 1) .- 2 .* n .^ (β .+ 1) .+ (n .+ 1) .^ (β .+ 1)) ./ Γ(β .+ 2)

    end

end

## Gamma function for vectors ##
Γ(b) = map(gamma, b)


text:
# in this overload we will work on some meta/debug information management
function CC.typeinf(interp::JETInterpreter, frame::InferenceState)
    linfo = frame.linfo

    #= logging start =#
    local sec, depth
    logger_activated = !isnothing(JETLogger(interp).inference_logger)
    if logger_activated
        sec = time()
        with_inference_logger(interp, ==(DEBUG_LOGGER_LEVEL)) do io
            depth = interp.depth

            print_rails(io, depth)
            printstyled(io, "┌ @ "; color = RAIL_COLORS[(depth+1)%N_RAILS+1])
            print(io, linfo)
            if is_constant_propagated(frame)
                printstyled(io, " (constant prop': ", frame.result.argtypes, ')'; color = NOERROR_COLOR)
            end
            file, line = get_file_line(linfo)
            print(io, ' ', file, ':', line)
            println(io)
            interp.depth += 1 # manipulate this only in debug mode
        end
    end
    #= logging end =#

    prev_frame = interp.current_frame
    interp.current_frame = frame

    ret = @invoke typeinf(interp::AbstractInterpreter, frame::InferenceState)

    interp.current_frame = prev_frame

    #= logging start =#
    if logger_activated
        sec = round(time() - sec; digits = 3)
        with_inference_logger(interp, ==(INFO_LOGGER_LEVEL)) do io
            println(io, "inference on $linfo finished in $sec sec")
        end
        with_inference_logger(interp, ==(DEBUG_LOGGER_LEVEL)) do io
            print_rails(io, depth)
            printstyled(io, "└─→ "; color = RAIL_COLORS[(depth+1)%N_RAILS+1])
            result = get_result(frame)
            isa(result, InferenceState) || printstyled(io, result; color = TYPE_ANNOTATION_COLOR)
            println(io, " (",
                        join(filter(!isnothing, (
                             linfo,
                             ret ? nothing : "in cycle",
                             "$(length(interp.reports)) reports",
                             "$sec sec"
                             )), ", "),
                        ')')
            interp.depth -= 1 # manipulate this only in debug mode
        end
    end
    #= logging end =#

    return ret
end

# TODO: disable optimization for better performance, only do necessary analysis work by ourselves

# in this overload we can work on `frame.src::CodeInfo` (and also `frame::InferenceState`)
# where type inference (and also optimization if applied) already ran on
function CC._typeinf(interp::JETInterpreter, frame::InferenceState)
    linfo = frame.linfo
    parent = frame.parent
    isentry = parent === nothing
    iscp = is_constant_propagated(frame)
    reports = interp.reports

    # some methods like `getproperty` can't propagate accurate types without actual values,
    # and constant prop' plays a somewhat critical role in those cases by overwriteing the
    # previous non-constant inference result (under the current design constant prop' always
    # happens after inference with non-constant abstract elements)
    # JET also needs that in order to reduce false positive reports, and here we will
    # throw-away previously-collected error reports that are "lineage" of this frame,
    # when it is being re-inferred with constants
    # NOTE `frame.linfo` is the exactly same object as that of the previous non-constant inference
    # IDEA we may still want to keep some "serious" error reports like `GlobalUndefVarErrorReport`
    # even when constant prop' reveals it never happens given the current constant arguments
    iscp && !isentry && filter!(!is_from_same_frame(parent.linfo, linfo), reports)

    reports_before             = Set(reports)
    uncaught_exceptions_before = Set(interp.uncaught_exceptions)

    ret = @invoke _typeinf(interp::AbstractInterpreter, frame::InferenceState)

    stmts = frame.src.code

    # report (local) undef var error
    # this only works when optimization is enabled, just because `:throw_undef_if_not` and
    # `:(unreachable)` are introduced by `optimize`
    for (idx, stmt) in enumerate(stmts)
        if isa(stmt, Expr) && stmt.head === :throw_undef_if_not
            sym, _ = stmt.args

            # slots in toplevel frame may be a abstract global slot
            istoplevel(interp, frame) && is_global_slot(interp, sym) && continue

            next_idx = idx + 1
            if checkbounds(Bool, stmts, next_idx) && is_unreachable(@inbounds stmts[next_idx])
                # the optimization so far has found this statement is never "reachable";
                # JET reports it since it will invoke undef var error at runtime, or will just
                # be dead code otherwise
                report!(interp, LocalUndefVarErrorReport(interp, frame, sym))
            # else
                # by excluding this pass, JET accepts some false negatives (i.e. don't report
                # those that may actually happen on actual execution)
            end
        end
    end

    # XXX this is a dirty fix for performance problem, we need more "proper" fix
    # https://github.com/aviatesk/JET.jl/issues/75
    unique!(report_identity_key, reports)

    reports_after = Set(reports)

    # report `throw` calls "appropriately"
    if get_result(frame) === Bottom
        # if the return type here is `Bottom` annotated, this _may_ mean there're uncaught
        # `throw` calls
        # XXX: well, it's possible that the `throw` calls within them are all caught but the
        # other critical errors make the return type `Bottom`
        # NOTE: to reduce the false positive `UncaughtExceptionReport`s described above, we count
        # `throw` calls here after optimization, since it may have eliminated "unreachable"
        # `throw` calls
        codelocs    = frame.src.codelocs
        linetable   = frame.src.linetable::Vector
        throw_locs  = LineInfoNode[]
        throw_calls = Expr[]
        for r in reports
            if isa(r, ExceptionReport) && last(r.st).linfo === linfo
                push!(throw_locs, r.lin)
            end
        end
        for (i, stmt) in enumerate(stmts)
            is_throw_call_expr(interp, frame, stmt) || continue
            # if this `throw` is already reported, don't duplciate
            linetable[codelocs[i]]::LineInfoNode in throw_locs && continue
            push!(throw_calls, stmt)
        end
        if !isempty(throw_calls)
            stash_uncaught_exception!(interp, UncaughtExceptionReport(interp, frame, throw_calls))
        end
    else
        # the non-`Bottom` result here may mean `throw` calls from the children frames
        # (if exists) are caught and not propagated here;
        # we don't want to cache `UncaughtExceptionReport`s for those calls for this frame
        # and its parents, so just filter them away
        empty!(interp.uncaught_exceptions)
    end

    uncaught_exceptions_after = Set(interp.uncaught_exceptions)

    # compute JET analysis results that should be cached for this linfo
    this_caches = union!(setdiff!(reports_after, reports_before),
                         setdiff!(uncaught_exceptions_after, uncaught_exceptions_before))

    if !isempty(this_caches)
        if iscp
            result = frame.result
            argtypes = result.argtypes
            cache = interp.cache
            @static JET_DEV_MODE && @assert jet_cache_lookup(linfo, argtypes, cache) === nothing "invalid local caching $linfo, $argtypes"
            local_cache = InferenceErrorReportCache[]
            for report in this_caches
                # # TODO make this hold
                # @assert first(report.st).linfo === linfo "invalid local caching"
                cache_report!(local_cache, report)
            end
            # branching on https://github.com/JuliaLang/julia/pull/39972
            given_argtypes, overridden_by_const = @static if VERSION ≥ v"1.7.0-DEV.705"
                def = result.linfo.def
                va_overwride = isa(def, Method) && def.is_for_opaque_closure
                matching_cache_argtypes(linfo, argtypes, va_overwride)
            else
                matching_cache_argtypes(linfo, argtypes)
            end
            push!(cache, AnalysisResult(linfo, given_argtypes, overridden_by_const, local_cache))
        elseif frame.cached # only cache when `NativeInterpreter` does
            cache = jet_report_cache(interp)
            @static JET_DEV_MODE && @assert !haskey(cache, linfo) || isentry "invalid global caching $linfo"
            global_cache = InferenceErrorReportCache[]
            for report in this_caches
                # # TODO make this hold
                # @assert first(report.st).linfo === linfo "invalid global caching"
                cache_report!(global_cache, report)
            end
            cache[linfo] = global_cache
        end
    end

    interp.to_be_updated = this_caches

    if !iscp && !isentry
        # refinement for this `linfo` may change analysis result for parent frame
        # XXX: is this okay from performance perspective ?
        add_backedge!(linfo, parent)
    end

    return ret
end

function CC.typeinf_local(interp::JETInterpreter, frame::InferenceState)
    ret = @invoke CC.typeinf_local(interp::AbstractInterpreter, frame::InferenceState)
    stmts = frame.src.code

    iter_types = Dict{Int64, Vector{CC.VarState}}()  # a map from variable ids to their possible types
    for (stmt, typ) in zip(stmts, frame.stmt_types)
        if isa(stmt, Expr) && length(stmt.args) >= 2
            sym = stmt.args[2]
            if (isa(sym, Expr) && !isempty(sym.args) &&
                isa(sym.args[1], GlobalRef) && sym.args[1].name === :iterate)
                iter_var = sym.args[2].id  # e.g. the `_1` in `_1 = Base.iterate(%1, %8)`
                types = filter(t -> isa(t, VarState), typ)
                @info(iter_var)
                if !haskey(iter_types, iter_var)
                    iter_types[iter_var] = types
                else
                    append!(iter_types[iter_var], types)
                end
            elseif isa(sym, Core.SlotNumber) && haskey(iter_types, sym.id)
                append!(iter_types[sym.id], typ)
            end
        end
    end
    for (var, types) in pairs(iter_types)
        if (!isempty(types) &&
            !any(t -> typeintersect(t.typ, Union{Nothing, Tuple{Any, Any}}) !== Bottom, types))
            report!(interp, InfiniteIterationErrorReport(interp, frame))
        end
    end

    return ret
end

"""
    is_from_same_frame(parent_linfo::MethodInstance, current_linfo::MethodInstance) ->
        (report::InferenceErrorReport) -> Bool

Returns a function that checks if a given `InferenceErrorReport` is generated from `current_linfo`.
It also checks `current_linfo` is a "lineage" of `parent_linfo` (i.e. entered from it).

This function is supposed to be used to filter out reports collected from analysis on `current_linfo`
  without using constants when entering into the constant analysis. As such, this function
  assumes that when a report should be filtered out, the first elment of its virtual stack
  frame `st` is for `parent_linfo` and the second element of that is for `current_linfo`.

Example: Assume `linfo2` will produce a report for some reason.
```
entry
└─ linfo1
   ├─ linfo2 (report1: linfo2)
   ├─ linfo3 (report1: linfo1->linfo2, report2: linfo3->linfo2)
   │  └─ linfo2 (report1: linfo1->linfo2, report2: linfo2)
   └─ linfo3′ (report1: linfo1->linfo2, ~~report2: linfo1->linfo3->linfo2~~)
```
In the example analysis above, `report2` will be filtered out on re-entering into `linfo3′`
  (i.e. we're analyzing `linfo3` with constants argument), because
  `is_from_same_frame(linfo1, linfo3)(report2)` returns `true`.
Note that `report1` is still kept there because of the lineage check, i.e.
  `is_from_same_frame(linfo1, linfo3)(report1)` returns `false`.
"""
function is_from_same_frame(parent_linfo::MethodInstance,
                            current_linfo::MethodInstance,
                            )
    function (report::InferenceErrorReport)
        @inbounds begin
            st = report.st
            length(st) > 1 || return false
            st[1].linfo === parent_linfo || return false
            return st[2].linfo === current_linfo
        end
    end
end

is_unreachable(@nospecialize(x)) = isa(x, ReturnNode) && !isdefined(x, :val)

@withmixedhash struct ReportIdentityKey
    T::Type{<:InferenceErrorReport}
    sig::Vector{Any}
    # entry_frame::VirtualFrame
    error_frame::VirtualFrame
end

report_identity_key(report::T) where {T<:InferenceErrorReport} =
    ReportIdentityKey(T, report.sig, #=first(report.st),=# last(report.st))

# basically same as `is_throw_call`, but also toplevel module handling added
function is_throw_call_expr(interp::JETInterpreter, frame::InferenceState, @nospecialize(e))
    if isa(e, Expr)
        if e.head === :call
            f = e.args[1]
            if istoplevel(interp, frame) && isa(f, Symbol)
                f = GlobalRef(interp.toplevelmod, f)
            end
            if isa(f, GlobalRef)
                ff = CC.abstract_eval_global(f.mod, f.name)
                if isa(ff, Const) && ff.val === Core.throw
                    return true
                end
            end
        end
    end
    return false
end


text:
module SrkMethodsForDegenerateLagrangianSystems

    include("common.jl")
    include("tableau_lists.jl")

    export tableaus_vprk_glrk,
           tableaus_vprk_lobatto_IIIA_IIIB,
           tableaus_vprk_lobatto_IIIB_IIIA,
           tableaus_vprk_radau,
           tableaus_srk_glrk,
           tableaus_firk_glrk

end


text:
using Convex

m = 4; n = 5
A = randn(m,n); b = randn(m, 1)


text:
render_depwarn_msg = "This version of the render function is deprecated since v0.8. Please use render!(renderables) instead."

function render(scene::EntityFrame{S,D,I}, roadway::R, overlays::AbstractVector{O};
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{I,C}=Dict{I,Colorant}(),
    surface::CairoSurface = CairoSVGSurface(IOBuffer(), canvas_width, canvas_height)
) where {S,D,I,O<:SceneOverlay,R,C<:Colorant}
    Base.depwarn(render_depwarn_msg, :render)
    renderables = [roadway]
    for (i, veh) in enumerate(scene)
        c = car_colors[i]
        if rendermode == :fancy
            r = (
                class(veh.def) == AgentClass.PEDESTRIAN
                ? FancyPedestrian(ped=veh, color=c)
                : FancyCar(car=veh, color=c)
            )
        else
            x, y, theta = posg(veh)
            r = ArrowCar(x, y, theta, length=length(veh.def), width=width(veh.def), color=c, id=veh.id)
        end
        push!(renderables, r)
    end
    for o in overlays
        push!(renderables, RenderableOverlay(o, scene, roadway))
    end
    update_camera!(cam, scene)
    render(renderables, camera=cam, surface=surface)
    return surface
end

function render!(
    rendermodel::RenderModel,
    scene::EntityFrame{S,D,I};
    car_color::Colorant=colortheme["COLOR_CAR_OTHER"], # default color
    car_colors::Dict{I,C}=Dict{I,Colorant}(), #  id -> color
) where {S,D,I,C<:Colorant}
    Base.depwarn(render_depwarn_msg, :render)
    renderables = []
    for (i, veh) in enumerate(scene)
        c = car_colors[i]
        if rendermode == :fancy
            r = (
                class(veh.def) == AgentClass.PEDESTRIAN
                ? FancyPedestrian(ped=veh, color=c)
                : FancyCar(car=veh, color=c)
            )
        else
            x, y, theta = posg(veh)
            r = ArrowCar(x, y, theta, length=length(veh.def), width=width(veh.def), color=c, id=veh.id)
        end
        push!(renderables, r)
    end
    render(renderables, camera=cam, surface=surface)
    return rendermodel
end

function render(roadway::R;
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel = RenderModel(),
    cam::Camera = SceneFollowCamera(),
    surface::CairoSurface = CairoSVGSurface(IOBuffer(), canvas_width, canvas_height)
) where {R<:Roadway}
    Base.depwarn(render_depwarn_msg, :render)
    renderables = [roadway]
    render(renderables, camera=cam, surface=surface)
    return surface
end

function render(ctx::CairoContext, scene::EntityFrame{S,D,I}, roadway::R;
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{I,C}=Dict{I,Colorant}(),
) where {S,D,I,R,C<:Colorant}
    Base.depwarn(render_depwarn_msg, :render)
    renderables = [roadway]
    for (i, veh) in enumerate(scene)
        c = car_colors[i]
        if rendermode == :fancy
            r = (
                class(veh.def) == AgentClass.PEDESTRIAN
                ? FancyPedestrian(ped=veh, color=c)
                : FancyCar(car=veh, color=c)
            )
        else
            x, y, theta = posg(veh)
            r = ArrowCar(x, y, theta, length=length(veh.def), width=width(veh.def), color=c, id=veh.id)
        end
        push!(renderables, r)
    end
    reset_instructions!(rendermodel)
    update_camera!(cam, scene)
    for renderable in renderables
        add_renderable!(rendermodel, renderable)
    end
    render_to_canvas(rendermodel, cam, ctx)
    ctx
end
function render(scene::EntityFrame{S,D,I}, roadway::R;
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel::RenderModel=RenderModel(),
    cam::Camera=SceneFollowCamera(),
    car_colors::Dict{I,C}=Dict{I,Colorant}(), # id
    surface::CairoSurface = CairoSVGSurface(IOBuffer(), canvas_width, canvas_height)
) where {S,D,I,R, C<:Colorant}
    Base.depwarn(render_depwarn_msg, :render)
    ctx = creategc(surface)
    render(ctx, scene, roadway, rendermodel=rendermodel, cam=cam, car_colors=car_colors)
    return surface
end

function render(roadway::StraightRoadway;
    canvas_width::Int=DEFAULT_CANVAS_WIDTH,
    canvas_height::Int=DEFAULT_CANVAS_HEIGHT,
    rendermodel = RenderModel(),
    cam::Camera = SceneFollowCamera(),
    surface::CairoSurface = CairoSVGSurface(IOBuffer(), canvas_width, canvas_height)
)
    Base.depwarn(render_depwarn_msg, :render)
    renderables = [roadway]
    update_camera!(cam, scene)
    render(renderables, camera=cam, surface=surface)
    return surface
end

Base.show(io::IO, ::MIME"image/png", roadway::StraightRoadway) = show(io, MIME"image/png"(), render(roadway))


text:
simple_network = QNetwork()
A = BasicNode("A")
C = BasicNode("C")
B = BasicNode("B")
AB = BasicChannel(A, B, exp_cost=false)
AC = BasicChannel(A, C, exp_cost=false)
CB = BasicChannel(C, B, exp_cost=false)
AB.costs = unit_costvector()
AC.costs = unit_costvector()
CB.costs = unit_costvector()

for i in [A, B, C, AB, AC, CB]
    add(simple_network, i)
end

refresh_graph!(simple_network)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
#
# Function that tries to read a number as an integer given a number
# or, perhaps, an hexadecimal representation string
#
function parse_int(s::String)
  try
    i = parse(Int, s)
    return i
  catch
    i = parse(Int,s,base=16)
    return i
  end
  error("Could not read integer from string: \"$s\"")
end