text:
function flux_quad_GxI_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)
    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]

    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    sqrt_beta = CUDAnative.sqrt(beta)
    S1 = ut*sqrt_beta
    S2 = un*sqrt_beta
    sqrt_pi_beta = CUDAnative.sqrt(pi*beta)
    B1 = 0.5*CUDAnative.exp(-S1*S1)/sqrt_pi_beta
    B2 = 0.5*CUDAnative.exp(-S2*S2)/sqrt_pi_beta
    A1neg = 0.5*(1.0 - CUDAnative.erf(S1))
    A2neg = 0.5*(1.0 - CUDAnative.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un
    shared[thread_idx] = op((rho*A2neg*(ut*A1neg - B1)), shared[thread_idx])

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1neg-ut*B1
    shared[thread_idx + block_dim] = op((rho*A2neg*temp2), shared[thread_idx + block_dim])

    temp1 = ut*A1neg - B1
    temp2 = un*A2neg - B2
    shared[thread_idx + block_dim * 2] = op((rho*temp1*temp2), shared[thread_idx + block_dim * 2])

    temp1 = (7.0 *pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1neg

    temp1 = (6.0 *pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1
    temp1 = ut*A1neg - B1
    temp4 = 0.5*rho*un*B2*temp1

    shared[thread_idx + block_dim * 3] = op((rho*A2neg*(temp2 - temp3) - temp4), shared[thread_idx + block_dim * 3])
    return nothing
end

function flux_quad_GxII_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)
    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]

    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr

    sqrt_beta = CUDAnative.sqrt(beta)
    S1 = ut*sqrt_beta
    S2 = un*sqrt_beta
    sqrt_pi_beta = CUDAnative.sqrt(pi*beta)
    B1 = 0.5*CUDAnative.exp(-S1*S1)/sqrt_pi_beta
    B2 = 0.5*CUDAnative.exp(-S2*S2)/sqrt_pi_beta
    A1pos = 0.5*(1.0 + CUDAnative.erf(S1))
    A2neg = 0.5*(1.0 - CUDAnative.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un
    shared[thread_idx] = op(rho * A2neg* (ut*A1pos + B1), shared[thread_idx])

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1pos + ut*B1
    shared[thread_idx + block_dim] = op((rho*A2neg*temp2), shared[thread_idx + block_dim])

    temp1 = ut*A1pos + B1
    temp2 = un*A2neg - B2
    shared[thread_idx + block_dim * 2] = op((rho*temp1*temp2), shared[thread_idx + block_dim * 2])

    temp1 = (7 *pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1pos

    temp1 = (6 *pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1

    temp1 = ut*A1pos + B1
    temp4 = 0.5*rho*un*B2*temp1

    shared[thread_idx + block_dim * 3] = op((rho*A2neg*(temp2 + temp3) - temp4), shared[thread_idx + block_dim * 3])
    return nothing
end

function flux_quad_GxIII_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)

    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]
    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr
    sqrt_beta = CUDAnative.sqrt(beta)
    S1 = ut*sqrt_beta
    S2 = un*sqrt_beta
    sqrt_pi_beta = CUDAnative.sqrt(pi*beta)
    B1 = 0.5*CUDAnative.exp(-S1*S1)/sqrt_pi_beta
    B2 = 0.5*CUDAnative.exp(-S2*S2)/sqrt_pi_beta
    A1pos = 0.5*(1.0 + CUDAnative.erf(S1))
    A2pos = 0.5*(1.0 + CUDAnative.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un
    shared[thread_idx] = op(rho*A2pos*(ut*A1pos + B1), shared[thread_idx])

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1pos + ut*B1
    shared[thread_idx + block_dim] = op((rho*A2pos*temp2), shared[thread_idx + block_dim])

    temp1 = ut*A1pos + B1
    temp2 = un*A2pos + B2
    shared[thread_idx + block_dim * 2] = op((rho*temp1*temp2), shared[thread_idx + block_dim * 2])

    temp1 = (7*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1pos
    temp1 = (6*pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1
    temp1 = ut*A1pos + B1
    temp4 = 0.5*rho*un*B2*temp1

    shared[thread_idx + block_dim * 3] = op((rho*A2pos*(temp2 + temp3) + temp4), shared[thread_idx + block_dim * 3])
    return nothing
end

function flux_quad_GxIV_kernel(nx, ny, idx, shared, op::Function, thread_idx, block_dim)

    u1 = shared[thread_idx + block_dim * 4]
    u2 = shared[thread_idx + block_dim * 5]
    rho = shared[thread_idx + block_dim * 6]
    pr = shared[thread_idx + block_dim * 7]
    ut = u1*ny - u2*nx
    un = u1*nx + u2*ny

    beta = 0.5*rho/pr

    sqrt_beta = CUDAnative.sqrt(beta)
    S1 = ut*sqrt_beta
    S2 = un*sqrt_beta
    sqrt_pi_beta = CUDAnative.sqrt(pi*beta)
    B1 = 0.5*CUDAnative.exp(-S1*S1)/sqrt_pi_beta
    B2 = 0.5*CUDAnative.exp(-S2*S2)/sqrt_pi_beta
    A1neg = 0.5*(1.0 - CUDAnative.erf(S1))
    A2pos = 0.5*(1.0 + CUDAnative.erf(S2))

    pr_by_rho = pr/rho
    u_sqr = ut*ut + un*un

    shared[thread_idx] = op((rho*A2pos*(ut*A1neg - B1)), shared[thread_idx])

    temp1 = pr_by_rho + ut*ut
    temp2 = temp1*A1neg - ut*B1
    shared[thread_idx + block_dim] = op((rho*A2pos*temp2), shared[thread_idx + block_dim])

    temp1 = ut*A1neg - B1
    temp2 = un*A2pos + B2
    shared[thread_idx + block_dim * 2] = op((rho*temp1*temp2), shared[thread_idx + block_dim * 2])

    temp1 = (7.0*pr_by_rho) + u_sqr
    temp2 = 0.5*ut*temp1*A1neg

    temp1 = (6.0 *pr_by_rho) + u_sqr
    temp3 = 0.5*B1*temp1

    temp1 = ut*A1neg - B1
    temp4 = 0.5*rho*un*B2*temp1

    shared[thread_idx + block_dim * 3] = op((rho*A2pos*(temp2 - temp3) + temp4), shared[thread_idx + block_dim * 3])
    return nothing
end

# print(flux_quad_GxI(1,2,3,4,5,5))

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


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

struct with inner constructor:
```julia
function PV(;
    tilt::Real,
    array_type::Int=1,
    module_type::Int=0,
    losses::Real=0.14,
    azimuth::Real=180,
    gcr::Real=0.4,
    radius::Int=0,
    name::String="PV",
    location::String="both",
    existing_kw::Real=0,
    min_kw::Real=0,
    max_kw::Real=1.0e9,
    installed_cost_per_kw::Real=1592.0,
    om_cost_per_kw::Real=17.0,
    degradation_pct::Real=0.005,
    macrs_option_years::Int = 5,
    macrs_bonus_pct::Real = 1.0,
    macrs_itc_reduction::Real = 0.5,
    kw_per_square_foot::Float64=0.01,
    acres_per_kw::Float64=6e-3,
    inv_eff::Float64=0.96,
    dc_ac_ratio::Real=1.2,
    prod_factor_series::Union{Missing, Array{Real,1}} = missing,
    federal_itc_pct::Real = 0.26,
    federal_rebate_per_kw::Real = 0.0,
    state_ibi_pct::Real = 0.0,
    state_ibi_max::Real = 1.0e10,
    state_rebate_per_kw::Real = 0.0,
    state_rebate_max::Real = 1.0e10,
    utility_ibi_pct::Real = 0.0,
    utility_ibi_max::Real = 1.0e10,
    utility_rebate_per_kw::Real = 0.0,
    utility_rebate_max::Real = 1.0e10,
    production_incentive_per_kwh::Real = 0.0,
    production_incentive_max_benefit::Real = 1.0e9,
    production_incentive_years::Int = 1,
    production_incentive_max_kw::Real = 1.0e9
    can_net_meter::Bool = off_grid_flag ? false : true,
    can_wholesale::Bool = off_grid_flag ? false : true,
    can_export_beyond_nem_limit::Bool = off_grid_flag ? false : true,
    can_curtail::Bool = true,
    operating_reserve_required_pct::Real = off_grid_flag ? 0.25 : 0.0
)
```
!!! note
    If `tilt` is not provided then it is set to the `Site.latitude`. (Which is handled in the `Scenario` struct.)
"""
struct PV <: AbstractTech
    tilt
    array_type
    module_type
    losses
    azimuth
    gcr
    radius
    name
    location
    existing_kw
    min_kw
    max_kw
    installed_cost_per_kw
    om_cost_per_kw
    degradation_pct
    macrs_option_years
    macrs_bonus_pct
    macrs_itc_reduction
    kw_per_square_foot
    acres_per_kw
    inv_eff
    dc_ac_ratio
    prod_factor_series
    federal_itc_pct
    federal_rebate_per_kw
    state_ibi_pct
    state_ibi_max
    state_rebate_per_kw
    state_rebate_max
    utility_ibi_pct
    utility_ibi_max
    utility_rebate_per_kw
    utility_rebate_max
    production_incentive_per_kwh
    production_incentive_max_benefit
    production_incentive_years
    production_incentive_max_kw
    can_net_meter
    can_wholesale
    can_export_beyond_nem_limit
    can_curtail
    operating_reserve_required_pct

    function PV(;
        off_grid_flag::Bool = false,
        tilt::Real,
        array_type::Int=1,
        module_type::Int=0,
        losses::Real=0.14,
        azimuth::Real=180,
        gcr::Real=0.4,
        radius::Int=0,
        name::String="PV",
        location::String="both",
        existing_kw::Real=0,
        min_kw::Real=0,
        max_kw::Real=1.0e9,
        installed_cost_per_kw::Real=1592.0,
        om_cost_per_kw::Real=17.0,
        degradation_pct::Real=0.005,
        macrs_option_years::Int = 5,
        macrs_bonus_pct::Real = 1.0,
        macrs_itc_reduction::Real = 0.5,
        kw_per_square_foot::Float64=0.01,
        acres_per_kw::Float64=6e-3,
        inv_eff::Float64=0.96,
        dc_ac_ratio::Real=1.2,
        prod_factor_series::Union{Missing, Array{Real,1}} = missing,
        federal_itc_pct::Real = 0.26,
        federal_rebate_per_kw::Real = 0.0,
        state_ibi_pct::Real = 0.0,
        state_ibi_max::Real = 1.0e10,
        state_rebate_per_kw::Real = 0.0,
        state_rebate_max::Real = 1.0e10,
        utility_ibi_pct::Real = 0.0,
        utility_ibi_max::Real = 1.0e10,
        utility_rebate_per_kw::Real = 0.0,
        utility_rebate_max::Real = 1.0e10,
        production_incentive_per_kwh::Real = 0.0,
        production_incentive_max_benefit::Real = 1.0e9,
        production_incentive_years::Int = 1,
        production_incentive_max_kw::Real = 1.0e9,
        can_net_meter::Bool = off_grid_flag ? false : true,
        can_wholesale::Bool = off_grid_flag ? false : true,
        can_export_beyond_nem_limit::Bool = off_grid_flag ? false : true,
        can_curtail::Bool = true,
        operating_reserve_required_pct::Real = off_grid_flag ? 0.25 : 0.0, # if off grid, 25%, else 0%. Applied to each time_step as a % of PV generation.
        )

        if !(off_grid_flag) && !(operating_reserve_required_pct == 0.0)
            @warn "PV operating_reserve_required_pct apply only when off_grid_flag is True. Setting operating_reserve_required_pct to 0.0 for this on-grid analysis."
            operating_reserve_required_pct = 0.0
        end

        if off_grid_flag && (can_net_meter || can_wholesale || can_export_beyond_nem_limit)
            @warn "Net metering, wholesale, and grid exports are not possible for off-grid scenarios. Setting can_net_meter, can_wholesale, and can_export_beyond_nem_limit to False."
            can_net_meter = false
            can_wholesale = false
            can_export_beyond_nem_limit = false
        end

        # validate inputs
        invalid_args = String[]
        if !(0 <= azimuth < 360)
            push!(invalid_args, "azimuth must satisfy 0 <= azimuth < 360, got $(azimuth)")
        end
        if !(array_type in [0, 1, 2, 3, 4])
            push!(invalid_args, "array_type must be in [0, 1, 2, 3, 4], got $(array_type)")
        end
        if !(module_type in [0, 1, 2])
            push!(invalid_args, "module_type must be in [0, 1, 2], got $(module_type)")
        end
        if !(0.0 <= losses <= 0.99)
            push!(invalid_args, "losses must satisfy 0.0 <= losses <= 0.99, got $(losses)")
        end
        if !(0 <= tilt <= 90)
            push!(invalid_args, "tilt must satisfy 0 <= tilt <= 90, got $(tilt)")
        end
        if !(location in ["roof", "ground", "both"])
            push!(invalid_args, "location must be in [\"roof\", \"ground\", \"both\"], got $(location)")
        end
        if !(0.0 <= degradation_pct <= 1.0)
            push!(invalid_args, "degradation_pct must satisfy 0 <= degradation_pct <= 1, got $(degradation_pct)")
        end
        if !(0.0 <= inv_eff <= 1.0)
            push!(invalid_args, "inv_eff must satisfy 0 <= inv_eff <= 1, got $(inv_eff)")
        end
        if !(0.0 <= dc_ac_ratio <= 2.0)
            push!(invalid_args, "dc_ac_ratio must satisfy 0 <= dc_ac_ratio <= 1, got $(dc_ac_ratio)")
        end
        # TODO validate additional args
        if length(invalid_args) > 0
            error("Invalid argument values: $(invalid_args)")
        end

        new(
            tilt,
            array_type,
            module_type,
            losses,
            azimuth,
            gcr,
            radius,
            name,
            location,
            existing_kw,
            min_kw,
            max_kw,
            installed_cost_per_kw,
            om_cost_per_kw,
            degradation_pct,
            macrs_option_years,
            macrs_bonus_pct,
            macrs_itc_reduction,
            kw_per_square_foot,
            acres_per_kw,
            inv_eff,
            dc_ac_ratio,
            prod_factor_series,
            federal_itc_pct,
            federal_rebate_per_kw,
            state_ibi_pct,
            state_ibi_max,
            state_rebate_per_kw,
            state_rebate_max,
            utility_ibi_pct,
            utility_ibi_max,
            utility_rebate_per_kw,
            utility_rebate_max,
            production_incentive_per_kwh,
            production_incentive_max_benefit,
            production_incentive_years,
            production_incentive_max_kw,
            can_net_meter,
            can_wholesale,
            can_export_beyond_nem_limit,
            can_curtail,
            operating_reserve_required_pct
        )
    end
end

function get_pv_by_name(name::String, pvs::AbstractArray{PV, 1})
    pvs[findfirst(pv -> pv.name == name, pvs)]
end

ext: jl

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
using FileIO
using JLD2

struct JLD2Writer <: AbstractWriter end

function write_data(jld::JLD2Writer, filename, dims, varvals, simtime)
    jldopen(full_name(jld, filename), "a+") do ds
        dimnames = collect(keys(dims))
        for di in 1:length(dimnames)
            ds["dim_$di"] = dimnames[di]
        end
        for (dn, dv) in dims
            ds[dn] = dv
        end
        ds["t"] = [1]
        for (vn, vv) in varvals
            ds[vn] = vv[2]
        end
        ds["simtime"] = [simtime]
    end
    return nothing
end


text:
"""
    PolyData

Defines cell types for polygonal datasets.

The following singleton types are defined:

- `PolyData.Verts` for vertices,
- `PolyData.Lines` for lines,
- `PolyData.Strips` for triangular strips,
- `PolyData.Polys` for polygons.

"""
module PolyData

import ..VTKCellTypes: nodes

abstract type CellType end

struct Verts <: CellType end
struct Lines <: CellType end
struct Strips <: CellType end
struct Polys <: CellType end

# All of these cell types can take any number of grid points.
# (This is for compatibility with VTKCellTypes for unstructured datasets.)
nodes(::CellType) = -1

number_attr(::Type{Verts}) = "NumberOfVerts"
number_attr(::Type{Lines}) = "NumberOfLines"
number_attr(::Type{Strips}) = "NumberOfStrips"
number_attr(::Type{Polys}) = "NumberOfPolys"

xml_node(::Type{Verts}) = "Verts"
xml_node(::Type{Lines}) = "Lines"
xml_node(::Type{Strips}) = "Strips"
xml_node(::Type{Polys}) = "Polys"

end

import .PolyData

const PolyCell{T} = MeshCell{T} where {T <: PolyData.CellType}

Base.eltype(::Type{T}) where {T <: PolyCell} = cell_type(T)
cell_type(::Type{<:PolyCell{T}}) where {T} = T
grid_type(::Type{<:PolyCell}) = VTKPolyData()

# Add different types of PolyCell (lines, vertices, ...) recursively.
function add_poly_cells!(vtk, xml_piece, cells::AbstractArray{<:PolyCell},
                         etc...)
    ctype = eltype(eltype(cells))
    add_cells!(vtk, xml_piece, PolyData.number_attr(ctype),
               PolyData.xml_node(ctype), cells, with_types=Val(false))
    add_poly_cells!(vtk, xml_piece, etc...)
end

add_poly_cells!(vtk, xml) = vtk

function vtk_grid(dtype::VTKPolyData, filename::AbstractString,
                  points::UnstructuredCoords,
                  cells::Vararg{AbstractArray{<:PolyCell}}; kwargs...)
    if isempty(cells)
        throw(ArgumentError(
            "constructing PolyData dataset with no cells is not allowed"))
    end

    Npts = num_points(dtype, points)
    Ncls = sum(length, cells)

    xvtk = XMLDocument()
    vtk = DatasetFile(dtype, xvtk, filename, Npts, Ncls; kwargs...)

    xroot = vtk_xml_write_header(vtk)
    xGrid = new_child(xroot, vtk.grid_type)

    xPiece = new_child(xGrid, "Piece")
    set_attribute(xPiece, "NumberOfPoints", Npts)

    xPoints = new_child(xPiece, "Points")
    data_to_xml(vtk, xPoints, points, "Points", 3)

    add_poly_cells!(vtk, xPiece, cells...)

    vtk
end


text:
"""
`gettag(tag, str)` obtain the `tag` from `str` that is embraced by `[]`. For example, `gettag("dt", "stn[NC]dt[20210101].mat")`
"""
function gettag(tag, str)
    expr = Regex("(?<=$tag\\[).*?(?=\\])")
    return match(expr, str).match
end


text:
struct Marks
    n_vert::Int64
    m::Int64
end
Base.eltype(::Type{Marks}) = Vector{Int}
Base.length(mark::Marks) = mark.n_vert^mark.m

#each element of the cycle for is an array of length m such that the position at i is the vector containing the mark I

function Base.iterate(mark::Marks, prev::Union{Nothing, Vector{Int}} = nothing)::Union{Nothing, Tuple{Vector{Int}, Vector{Int}}}

    prev === nothing && return (ones(Int,mark.m), ones(Int,mark.m))

    index = findfirst(x -> prev[x] < mark.n_vert, 1:mark.m)

    index === nothing && return nothing

    prev[index] += 1
    for i in 1:(index-1)
        @inbounds prev[i] = 1
    end

    return (prev, prev)
end

function invert_marks(mark::Vector{Int}, n_vert::Int)::Dict{Int64,Vector{Int64}}

    D = Dict{Int64,Vector{Int64}}([v for v in 1:n_vert] .=> [Int[] for _ in 1:n_vert])

    for t in 1:length(mark)
        push!(D[mark[t]], t)
    end
    return D
end

function num_marks(mark::Vector{Int}, v::Int64)::Int64

    length(mark) == 0 && return 0

    return count(t -> mark[t] == v, 1:length(mark))
end


text:
function compMises3D( σh::Array{Float64,3} )

    σmis = Array{Float64}( size(σh,1), 1, size(σh,3) )
    σmis[:,1,:] = sqrt.( 0.5 * ( (σh[:,1,:] - σh[:,5,:]).^2
                             + (σh[:,5,:] - σh[:,9,:]).^2
                             + (σh[:,9,:] - σh[:,1,:]).^2
                             + 6*(σh[:,2,:].^2 + σh[:,6,:].^2 + σh[:,7,:].^2) ) )

    return σmis

end

function compMises2D( σh::Array{Float64,3} )

    σmis = Array{Float64}( size(σh,1), 1, size(σh,3) )
    σmis[:,1,:] = sqrt.( (σh[:,1,:] - σh[:,4,:]).^2
                         + 3*σh[:,2,:].^2 )

    return σmis

end


text:
"""
    players(env::AbstractEnv)

Return an ordered iterable collection of integer indices for all players, starting with one.

This function is a *static property* of the environment; the value it returns should not change based on the state.

# Example

```julia
@provide players(::MyEnv) = 1:2
```
"""
function players end

"""
    player(env::AbstractEnv)

Return the index of the player who should play next in the environment.
"""
function player end

"""
    all_act!(env::AbstractEnv, actions::AbstractVector)

Take `actions` for all players and advance AbstractEnv `env` forward, and return rewards for all players.

Environments that support simultaneous actions by all players should implement this in addition to or instead of `act!`.
"""
function all_act! end

"""
    all_observe(env::AbstractEnv)

Return observations from the environment for all players.

Environments that support simultaneous actions by all players should implement this in addition to or instead of `observe`.
"""
function all_observe end

"""
    UtilityStyle(env)

Trait that allows an environment to declare certain properties about the relative utility for the players.

Possible returns are:
- `ZeroSum()`
- `ConstantSum()`
- `GeneralSum()`
- `IdenticalUtility()`

See the docstrings for each for more details.
"""
abstract type UtilityStyle end

"""
If `UtilityStyle(env) == ZeroSum()` then the sum of the rewards returned by `act!` is always zero.
"""
struct ZeroSum <: UtilityStyle end
"""
If `UtilityStyle(env) == ConstantSum()` then the sum of the rewards returned by `act!` will always be the same constant.
"""
struct ConstantSum <: UtilityStyle end
"""
If `UtilityStyle(env) == GeneralSum()`, the sum of rewards over a trajectory can take any form.
"""
struct GeneralSum <: UtilityStyle end
"""
If `UtilityStyle(env) == IdenticalUtility()`, all entries of the reward returned by `act!` will be identical for all players.
"""
struct IdenticalUtility <: UtilityStyle end


text:
#
# oldrules.jl --
#
# Implement simplification rules for LazyAlgebra.  This is the old version,
# kept for comparisons and benchmarking.  The new version is simpler and
# faster.
#
module OldRules

using LazyAlgebra

using LazyAlgebra:
    Adjoint,
    Composition,
    Composition,
    Direct,
    Gram,
    Inverse,
    InverseAdjoint,
    LinearMapping,
    Mapping,
    Scaled,
    Sum,
    identifier,
    multiplier,
    unveil,
    terms,
    to_tuple

#------------------------------------------------------------------------------
# RULES FOR SUMS OF MAPPINGS

"""
    OldRules.add(A, B)

yields a simplied sum of mappings `A + B`.

"""
add(A::Mapping, B::Mapping) =
    simplify_sum((split_sum(A)..., split_sum(B)...))

# `split_sum(A)` yields a tuple of the terms of `A` if it is a sum or just
# `(A,)` otherwise.
split_sum(A::Sum) = terms(A)
split_sum(A::Mapping) = (A,)

# `merge_sum(args...)` constructs a fully qualified sum.  It is assumed that
# the argument(s) have already been simplified.  An empty sum is forbidden
# because there is no universal neutral element ("zero") for the addition.
merge_sum(arg::Mapping) = arg
merge_sum(args::Mapping...) = merge_sum(args)
merge_sum(args::Tuple{}) = throw(ArgumentError("empty sum"))
merge_sum(args::Tuple{Mapping}) = args[1]
merge_sum(args::T) where {N,T<:NTuple{N,Mapping}} = Sum{N,T}(args)

# `simplify_sum(args...)` simplifies the sum of all terms in `args...` and
# returns a single term (possibly an instance of `Sum`).  It is assumed that
# the operator `+` is associative and commutative.

# Make a sum out of 0-1 terms.
simplify_sum(args::Tuple{}) = merge_sum(args)
simplify_sum(args::Tuple{Mapping}) = args[1]

# Make a sum out of N terms (with N ≥ 2).
function simplify_sum(args::NTuple{N,Mapping}) where {N}
    # First group terms corresponding to identical mappings (up to an optional
    # multiplier) to produce a single (possibly scaled) term per group.  FIXME:
    # The following algorithm scales as O(N²) which is probably not optimal.
    # Nevertheless, it never compares twice the same pair of arguments.
    terms = Array{Mapping}(undef, 0)
    flags = fill!(Array{Bool}(undef, N), true)
    i = 1
    while i != 0
        # Push next ungrouped argument.
        push!(terms, args[i])
        flags[i] = false

        # Find any other argument which is identical, possibly scaled, term.
        k = i + 1
        i = 0
        for j = k:N
            if flags[j]
                if simplify_sum!(terms, terms[end], args[j])
                    flags[j] = false
                elseif i == 0
                    # Next ungrouped argument to consider.
                    i = j
                end
            end
        end
    end

    # Make a sum out of the terms after having eliminated the zeros.
    return sort_sum(terms)
end

# `simplify_sum!(terms, A, B)` is a helper function for the `simplify_sum`
# method which attempts to make a trivial simplification for `A + B` when `A =
# λ⋅M` and `B = μ⋅M` for any numbers `λ` and `μ` and any mapping `M`.  If such
# a simplification can be done, the result `(λ + μ)⋅M` is stored as the last
# component of `terms` and `true` is returned; otherwise, `false` is returned.
simplify_sum!(terms::Vector{Mapping}, A::Mapping, B::Mapping) = false

function simplify_sum!(terms::Vector{Mapping},
                       A::Scaled{T}, B::Scaled{T}) where {T<:Mapping}
    unscaled(A) === unscaled(B) || return false
    @inbounds terms[end] = (multiplier(A) + multiplier(B))*unscaled(A)
    return true
end

function simplify_sum!(terms::Vector{Mapping},
                       A::Scaled{T}, B::T) where {T<:Mapping}
    unscaled(A) === B || return false
    @inbounds terms[end] = (multiplier(A) + one(multiplier(A)))*B
    return true
end

function simplify_sum!(terms::Vector{Mapping},
                       A::T, B::Scaled{T}) where {T<:Mapping}
    A === unscaled(B) || return false
    @inbounds terms[end] = (one(multiplier(B)) + multiplier(B))*A
    return true
end

function simplify_sum!(terms::Vector{Mapping},
                       A::T, B::T) where {T<:Mapping}
    A === B || return false
    @inbounds terms[end] = 2*A
    return true
end

# Sort the terms of a sum (so that all equivalent expressions eventually yield
# the same result after simplifications) and eliminate the "zeros"
# (if all terms are "zero", the sum simplifies to the first one).
sort_sum(args::Mapping...) = sort_sum(args)
sort_sum(args::Tuple{Mapping}) = args[1]
function sort_sum(terms::Union{Tuple{Vararg{Mapping}},Vector{<:Mapping}})
    perms = sortperm([identifier(terms[i]) for i in 1:length(terms)])
    n = 0
    @inbounds for i in 1:length(perms)
        j = perms[i]
        if ! iszero(terms[j])
            n += 1
            perms[n] = j
        end
    end
    if n ≤ 1
        # All terms are zero or only one term is non-zero, return the first
        # sorted term.
        return terms[perms[1]]
    else
        # Make a sum out of the remaing sorted terms.
        return merge_sum(ntuple(i -> terms[perms[i]], n))
    end
end

#------------------------------------------------------------------------------
# RULES FOR COMPOSITIONS OF MAPPINGS

"""
    OldRules.compose(A, B)

yields a simplified composition of mappings `A*B`.

"""
compose(A::Mapping, B::Mapping) = simplify_mul(A, B)

# `split_mul(A)` yields a tuple of the terms of `A` if it is a composition or
# just `(A,)` otherwise.
split_mul(A::Composition) = terms(A)
split_mul(A::Mapping) = (A,)

# `merge_mul(args...)` constructs a fully qualified composition.  It is assumed
# that the argument(s) have already been simplified.  An empty composition
# yields the identity which is the universal neutral element ("one") for the
# composition.
merge_mul(arg::Mapping) = arg
merge_mul(args::Mapping...) = merge_mul(args)
merge_mul(args::Tuple{}) = Id
merge_mul(args::Tuple{Mapping}) = args[1]
merge_mul(args::T) where {N,T<:NTuple{N,Mapping}} = Composition{N,T}(args)

# `simplify_mul(A,B)` simplifies the product of `A` and `B`.  The result is a
# tuple `C` of the resulting terms if `A` and `B` are both tuples of mappings
# or an instance of `Mapping` if `A` and `B` are both mappings.  If no
# simplification can be made, the result is just the concatenation of the terms
# in `A` and `B`.  To perform simplifications, it is assumed that, if they are
# compositions, `A` and `B` have already been simplified.  It is also assumed
# that the composition is associative but non-commutative.
simplify_mul(A::Mapping, B::Mapping) =
    merge_mul(simplify_mul(split_mul(A), split_mul(B)))

# The following versions of `simplify_mul` are for `A` and `B` in the form of
# tuples and return a tuple.  The algorithm is recursive and should works for
# any non-commutative binary operator.
simplify_mul(A::Tuple{}, B::Tuple{}) = (Id,)
simplify_mul(A::Tuple{Vararg{Mapping}}, B::Tuple{}) = A
simplify_mul(A::Tuple{}, B::Tuple{Vararg{Mapping}}) = B
function simplify_mul(A::NTuple{M,Mapping},
                       B::NTuple{N,Mapping}) where {M,N}
    # Here M ≥ 1 and N ≥ 1.
    @assert M ≥ 1 && N ≥ 1

    # Attempt to simplify the product of the terms at the junction.
    C = split_mul(A[M]*B[1])
    len = length(C)
    if len == 2
        if C === (A[M], B[1])
            # No simplification, just concatenate the 2 compositions.
            return (A..., B...)
        else
            # There have been some changes, but the result of A[M]*B[1] still
            # have two terms which cannot be further simplified.  So we
            # simplify its head with the remaining leftmost operands and its
            # tail with the remaining rightmost operands.
            L = simplify_mul(A[1:M-1], C[1]) # simplify leftmost operands
            R = simplify_mul(C[2], B[2:N])   # simplify rightmost operands
            if L[end] !== C[1] || R[1] !== C[2]
                # At least one of the last of resulting rightmost operands or
                # the first of the resulting leftmost operands has been
                # modified so there may be other possible simplifications.
                return simplify_mul(L, R)
            else
                # No further simplifications possible.
                return (L..., R...)
            end
        end
    elseif len == 1
        # Simplications have occured resulting in a single operand.  This
        # operand can be simplified whith the remaining leftmost operands
        # and/or with the remaining rightmost operands.  To benefit from the
        # maximum simplifications, we can either do:
        #
        #     simplify_mul(A[1:end-1], simplify_mul(C, B[2:end]))
        #
        # that is, simplify right then left, or:
        #
        #     simplify_mul(simplify_mul(A[1:end-1], C), B[2:end])
        #
        # that is simplify left then right.  Since we want to propagate
        # multipliers to the right of compositions, the former is the most
        # appropriate.
        return simplify_mul(A[1:end-1], simplify_mul(C, B[2:end]))
    else
        # len == 0 The result of A[M]*B[1] is the neutral element for * thus
        # eliminating these terms from the merging.  We just have to repeat the
        # process with the truncated associations in case further simplications
        # are possible.
        return simplify_mul(A[1:M-1], B[2:N-1])
    end
end

end # module


text:
# Reference implementations for summary functions.

function area(window)
    -(window.x...) * -(window.y...)
end

function pairdist(xy)
    d((a1, a2), (b1, b2)) = sqrt((a1-b1)^2 + (a2-b2)^2)
    d.(reshape(xy, 1, :), xy)
end
epanechnikov(t, h) = if -h <= t <= h 3/4*(1-(t/h)^2)/h else 0.0 end
function sumpairs(dist, A, r, window, f)
    # Use TwicePrecision because adding a lot of small numbers causes some
    # round off error.
    # Doesn't really matter because there are roundoff errors in other implementations.
    s = Base.TwicePrecision(0.0)
    for j in axes(dist, 2)
        for i in axes(dist, 1)
            if i != j
                s += f(r, dist[i,j])/A[i,j]
            end
        end
    end
    Float64(s)
end
function Atrans(window, d)
    d1, d2 = abs.(d)
    x1, x2 = window.x
    y1, y2 = window.y
    ((x2-x1-d1)*(y2-y1-d2))
end
function Apairs(xy, window)
    ((a, b) -> Atrans(window, a .- b)).(reshape(xy, 1, :), xy)
end
function pcfref(xy, window, r)
    λ = length(xy)/area(window)
    h = 0.15/sqrt(λ)
    dist = pairdist(xy)
    A = Apairs(xy, window)
    (r -> sumpairs(dist, A, r, window, (r, d) -> epanechnikov(r - d, h)) / (2 * pi * r) / λ^2 ).(r)
end

function Kref(xy, window, r)
    dist = pairdist(xy)
    A = Apairs(xy, window)
    n = length(xy)
    λ² = n*(n-1)/area(window)^2
    (r -> sumpairs(dist, A, r, window, (r, d) -> d < r) / λ² ).(r)
end


text:
#
# Proquint.jl by P. Bayer, MIT license
#
"""
# Proquints: Identifiers that are Readable, Spellable, and Pronounceable

This is a Julia implementation of Wilkerson's proquints
described in his article [http://arXiv.org/html/0901.4016](http://arXiv.org/html/0901.4016).

In short 16 bits get represented in a "proquint" as alternating
consonants and vowels.

Four-bits as a consonant:

    0 1 2 3 4 5 6 7 8 9 A B C D E F
    b d f g h j k l m n p r s t v z

Two-bits as a vowel:

    0 1 2 3
    a i o u

Whole 16-bit word, where "con" = consonant, "vo" = vowel:

     0 1 2 3 4 5 6 7 8 9 A B C D E F
    +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
    |con    |vo |con    |vo |con    |
    +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

Proquints are separated by dashes `-`.

## Short Form (extension)

The short form does not encode leading zeros of 16-bit words and
if an entire word is 0, it gets encoded as "x".

## Example

The `UInt32` `0x0a484904` is in bit representation:

    0000 10 1001 00 1000  0100 10 0100 00 0100
      b   o   n   a   m  -  h   o   h   a   h
          o   n   a   m  -  h   o   h   a   h   # short version

The IP address `127.0.0.1` is a `0x7f000001` `UInt32`:

    0111 11 1100 00 0000  0000 00 0000 00 0001
      l   u   s   a   b  -  b   a   b   a   d
      l   u   s   a   b  -                  d   # short version

Proquints are unambiguous and can be generated from integers
and converted back to them:

```julia
julia> using Proquint

julia> uint2quint(0x7f000001)
"lusab-babad"

julia> uint2quint(0x7f000001, short=true) # this gives the short version
"lusab-d"

julia> quint2uint("lusab-babad", UInt32)
0x7f000001

julia> quint2uint("lusab-d", UInt32)
0x7f000001
```
"""
module Proquint

"Gives the package version"
const version = v"0.1.0"

include("quint.jl")

export uint2quint, quint2uint

end


text:
export getSensMat

"""
S = function getSensMat(...)

constructs sensitivity matrix.

WARNING: For large-scale problems this will be prohibively
         expensive. Use with caution

Inputs:

	sigma - model
	pFor  - forward problems

Examples:

	S = getSensMat(sigma, pFor)            # single pFor

    Some methods of getData require further arguments.

"""
function getSensMat(sigma::Vector,
					pFor::ForwardProbType)

	(n, m) = getSensMatSize(pFor)
	sensMat = zeros(n, m)

	if min(m, n) > 1e4
		error("sensitivity matrix is too big to build column- or rowwise.")
	end

	if  n < m # decide which way is less work
		In = Matrix{Float64}(I, n, n);
		for k = 1:n
			tv = vec(getSensTMatVec(vec(In[:, k]), sigma, pFor))
			sensMat[k, :] = tv
		end
	else
		Im = Matrix{Float64}(I, m, m);
		for k = 1:m
			sensMat[:, k] = vec(getSensMatVec(vec(Im[:, k]), sigma, pFor))
		end
	end
	return sensMat
end

function getSensMat(sigma::Union{RemoteChannel, Vector},
                    pFor::ForwardProbType,
                    Mesh2Mesh::Mesh2MeshTypes)

    sig = interpGlobalToLocal(fetch(sigma), fetch(Mesh2Mesh))
    sensMat = getSensMat(sig, pFor)
    sensMat = remotecall(identity, myid(), sensMat)
    return sensMat
end

function getSensMat(sigma::Union{RemoteChannel, Vector},
                    pFor::RemoteChannel,
                    Mesh2Mesh::Mesh2MeshTypes=1.)

    pF = take!(pFor)
    sensMat = getSensMat(fetch(sigma), pF, fetch(Mesh2Mesh))
    put!(pFor, pF)
    return sensMat
end

function getSensMat(sigma::Vector,
                 					   param::Array{RemoteChannel},
                 					   Mesh2Mesh::Array{T}=ones(length(param))) where {T<:Mesh2MeshTypes}

    sensMat = Array{Future}(undef,length(param))
    workerList = getWorkerIds(param)
    sigmaRef = Array{RemoteChannel}(undef,maximum(workers()))
    @sync begin
        for p=workerList
            @async begin
                sigmaRef[p] = initRemoteChannel(identity, p, sigma)  # send model to worker
                for idx = 1:length(param)
                    if p == param[idx].where
                        sensMat[idx] = remotecall_fetch(getSensMat, p, sigmaRef[p], param[idx], Mesh2Mesh[idx])
                    end
                end
            end
        end
    end
    return sensMat
end

function getSensMat(sigma::Vector,
                                                             param::Array{FPT},
                                                             Mesh2Mesh::Array{T}=ones(length(param)),
                                                             workerList::Vector=workers()) where {FPT<:ForwardProbType, T<:Mesh2MeshTypes}
    i=1; nextidx() = (idx = i; i+=1; idx)

    sensMat = Array{Any}(undef,length(param))
    workerList = intersect(workers(),workerList)
    if isempty(workerList)
        error("getSensMat: workers do not exist!")
    end
    @sync begin
        for p=workerList
            @async begin
                while true
                    idx = nextidx()
                    if idx > length(param)
                        break
                    end
                    sensMat[idx] = remotecall_fetch(getSensMat,p,sigma,param[idx],Mesh2Mesh[idx])
                end
            end
        end
    end
    return sensMat
end


text:
# Rewriting the node graph to prepare it for translation.

# Auxiliary nodes.

# Get(over = Get(:a), name = :b) => NameBound(over = Get(:b), name = :a)
mutable struct NameBoundNode <: AbstractSQLNode
    over::SQLNode
    name::Symbol

    NameBoundNode(; over, name) =
        new(over, name)
end

NameBound(args...; kws...) =
    NameBoundNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(NameBound), pats::Vector{Any}) =
    dissect(scr, NameBoundNode, pats)

PrettyPrinting.quoteof(n::NameBoundNode, ctx::QuoteContext) =
    Expr(:call, nameof(NameBound), Expr(:kw, :over, quoteof(n.over, ctx)), Expr(:kw, :name, QuoteNode(n.name)))

# Get(over = q, name = :b) => HandleBound(over = Get(:b), handle = get_handle(q))
mutable struct HandleBoundNode <: AbstractSQLNode
    over::SQLNode
    handle::Int

    HandleBoundNode(; over, handle) =
        new(over, handle)
end

HandleBound(args...; kws...) =
    HandleBoundNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(HandleBound), pats::Vector{Any}) =
    dissect(scr, HandleBoundNode, pats)

PrettyPrinting.quoteof(n::HandleBoundNode, ctx::QuoteContext) =
    Expr(:call, nameof(NameBound), Expr(:kw, :over, quoteof(n.over, ctx)), Expr(:kw, :handle, n.handle))

mutable struct ExtendedBindNode <: AbstractSQLNode
    over::Union{SQLNode, Nothing}
    list::Vector{SQLNode}
    label_map::OrderedDict{Symbol, Int}
    owned::Bool     # Did we find the outer query for this node?

    function ExtendedBindNode(; over = nothing, list, label_map = nothing, owned = false)
        if label_map !== nothing
            return new(over, list, label_map, owned)
        end
        n = new(over, list, OrderedDict{Symbol, Int}())
        for (i, l) in enumerate(n.list)
            n.label_map[label(l)] = i
        end
        n
    end
end

ExtendedBind(args...; kws...) =
    ExtendedBindNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::ExtendedBindNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(ExtendedBind))
    if isempty(n.list)
        push!(ex.args, Expr(:kw, :list, Expr(:vect)))
    else
        append!(ex.args, quoteof(n.list, ctx))
    end
    push!(ex.args, Expr(:kw, :owned, n.owned))
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

rebase(n::ExtendedBindNode, n′) =
    ExtendedBindNode(over = rebase(n.over, n′),
                     list = n.list, label_map = n.label_map, owned = n.owned)

mutable struct ExtendedJoinNode <: TabularNode
    over::Union{SQLNode, Nothing}
    joinee::SQLNode
    on::SQLNode
    left::Bool
    right::Bool
    type::BoxType               # Type of the product of `over` and `joinee`.
    lateral::Vector{SQLNode}    # References from `joinee` to `over` for JOIN LATERAL.

    ExtendedJoinNode(; over, joinee, on, left, right, type = EMPTY_BOX, lateral = SQLNode[]) =
        new(over, joinee, on, left, right, type, lateral)
end

ExtendedJoinNode(joinee, on; over = nothing, left = false, right = false, type = EMPTY_BOX, lateral = SQLNode[]) =
    ExtendedJoinNode(over = over, joinee = joinee, on = on, left = left, right = right, type = type, lateral = lateral)

ExtendedJoin(args...; kws...) =
    ExtendedJoinNode(args...; kws...) |> SQLNode

function PrettyPrinting.quoteof(n::ExtendedJoinNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(ExtendedJoin))
    if !ctx.limit
        push!(ex.args, quoteof(n.joinee, ctx))
        push!(ex.args, quoteof(n.on, ctx))
        if n.left
            push!(ex.args, Expr(:kw, :left, n.left))
        end
        if n.right
            push!(ex.args, Expr(:kw, :right, n.right))
        end
        if n.type !== EMPTY_BOX
            push!(ex.args, Expr(:kw, :type, n.type))
        end
        if !isempty(n.lateral)
            push!(ex.args, Expr(:kw, :lateral, Expr(:vect, quoteof(n.lateral, ctx)...)))
        end
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

rebase(n::ExtendedJoinNode, n′) =
    ExtendedJoinNode(over = rebase(n.over, n′),
                     joinee = n.joinee, on = n.on, left = n.left, right = n.right, type = n.type, lateral = n.lateral)

label(n::Union{NameBoundNode, HandleBoundNode, ExtendedBindNode, ExtendedJoinNode}) =
    label(n.over)

# A SQL subquery with an undetermined SELECT list.
mutable struct BoxNode <: TabularNode
    over::Union{SQLNode, Nothing}
    type::BoxType
    handle::Int
    refs::Vector{SQLNode}

    BoxNode(; over = nothing, type = EMPTY_BOX, handle = 0, refs = SQLNode[]) =
        new(over, type, handle, refs)
end

Box(args...; kws...) =
    BoxNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Box), pats::Vector{Any}) =
    dissect(scr, BoxNode, pats)

function PrettyPrinting.quoteof(n::BoxNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Box))
    if !ctx.limit
        if n.type !== EMPTY_BOX
            push!(ex.args, Expr(:kw, :type, quoteof(n.type)))
        end
        if n.handle != 0
            push!(ex.args, Expr(:kw, :handle, n.handle))
        end
        if !isempty(n.refs)
            push!(ex.args, Expr(:kw, :refs, Expr(:vect, quoteof(n.refs, ctx)...)))
        end
    else
        push!(ex.args, :…)
    end
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::BoxNode) =
    n.type.name

rebase(n::BoxNode, n′) =
    BoxNode(over = rebase(n.over, n′),
            type = n.type, handle = n.handle, refs = n.refs)

box_type(n::BoxNode) =
    n.type

box_type(n::SQLNode) =
    box_type(n[]::BoxNode)

# Annotation context.

# Maps a node in the annotated graph to a path in the original graph (for error reporting).
struct PathMap
    paths::Vector{Tuple{SQLNode, Int}}
    origins::IdDict{Any, Int}

    PathMap() =
        new(Tuple{SQLNode, Int}[], IdDict{Any, Int}())
end

function get_path(map::PathMap, idx::Int)
    path = SQLNode[]
    while idx != 0
        n, idx = map.paths[idx]
        push!(path, n)
    end
    path
end

get_path(map::PathMap, n) =
    get_path(map, get(map.origins, n, 0))

struct AnnotateContext
    path_map::PathMap
    current_path::Vector{Int}
    handles::Dict{SQLNode, Int}
    boxes::Vector{BoxNode}

    AnnotateContext() =
        new(PathMap(), Int[0], Dict{SQLNode, Int}(), BoxNode[])
end

function grow_path!(ctx::AnnotateContext, n::SQLNode)
    push!(ctx.path_map.paths, (n, ctx.current_path[end]))
    push!(ctx.current_path, length(ctx.path_map.paths))
end

function shrink_path!(ctx::AnnotateContext)
    pop!(ctx.current_path)
end

function mark_origin!(ctx::AnnotateContext, n::SQLNode)
    ctx.path_map.origins[n] = ctx.current_path[end]
end

mark_origin!(ctx::AnnotateContext, n::AbstractSQLNode) =
    mark_origin!(ctx, convert(SQLNode, n))

get_path(ctx::AnnotateContext) =
    get_path(ctx.path_map, ctx.current_path[end])

get_path(ctx::AnnotateContext, n::SQLNode) =
    get_path(ctx.path_map, n)

function make_handle!(ctx::AnnotateContext, n::SQLNode)
    get!(ctx.handles, n) do
        length(ctx.handles) + 1
    end
end

function get_handle(ctx::AnnotateContext, n::SQLNode)
    handle = 0
    idx = get(ctx.path_map.origins, n, 0)
    if idx > 0
        n = ctx.path_map.paths[idx][1]
        handle = get(ctx.handles, n, 0)
    end
    handle
end

get_handle(ctx::AnnotateContext, ::Nothing) =
    0

# Rewriting of the node graph.

function annotate(n::SQLNode, ctx)
    grow_path!(ctx, n)
    n′ = convert(SQLNode, annotate(n[], ctx))
    mark_origin!(ctx, n′)
    box = BoxNode(over = n′)
    push!(ctx.boxes, box)
    n′ = convert(SQLNode, box)
    mark_origin!(ctx, n′)
    shrink_path!(ctx)
    n′
end

function annotate_scalar(n::SQLNode, ctx)
    grow_path!(ctx, n)
    n′ = convert(SQLNode, annotate_scalar(n[], ctx))
    mark_origin!(ctx, n′)
    shrink_path!(ctx)
    n′
end

annotate(ns::Vector{SQLNode}, ctx) =
    SQLNode[annotate(n, ctx) for n in ns]

annotate_scalar(ns::Vector{SQLNode}, ctx) =
    SQLNode[annotate_scalar(n, ctx) for n in ns]

function annotate(::Nothing, ctx)
    box = BoxNode()
    push!(ctx.boxes, box)
    n′ = convert(SQLNode, box)
    mark_origin!(ctx, n′)
    n′
end

annotate_scalar(::Nothing, ctx) =
    nothing

annotate(n::AbstractSQLNode, ctx) =
    throw(IllFormedError(path = get_path(ctx)))

function annotate_scalar(n::TabularNode, ctx)
    n′ = convert(SQLNode, annotate(n, ctx))
    mark_origin!(ctx, n′)
    box = BoxNode(over = n′)
    push!(ctx.boxes, box)
    n′ = convert(SQLNode, box)
    n′
end

function rebind(node, base, ctx)
    while @dissect node over |> Get(name = name)
        mark_origin!(ctx, base)
        base = NameBound(over = base, name = name)
        node = over
    end
    if node !== nothing
        handle = make_handle!(ctx, node)
        mark_origin!(ctx, base)
        base = HandleBound(over = base, handle = handle)
    end
    base
end

function annotate_scalar(n::AggregateNode, ctx)
    args′ = annotate_scalar(n.args, ctx)
    filter′ = annotate_scalar(n.filter, ctx)
    n′ = Agg(name = n.name, distinct = n.distinct, args = args′, filter = filter′)
    rebind(n.over, n′, ctx)
end

function annotate(n::AppendNode, ctx)
    over′ = annotate(n.over, ctx)
    list′ = annotate(n.list, ctx)
    Append(over = over′, list = list′)
end

function annotate(n::AsNode, ctx)
    over′ = annotate(n.over, ctx)
    As(over = over′, name = n.name)
end

function annotate_scalar(n::AsNode, ctx)
    over′ = annotate_scalar(n.over, ctx)
    As(over = over′, name = n.name)
end

function annotate(n::BindNode, ctx)
    over′ = annotate(n.over, ctx)
    list′ = annotate_scalar(n.list, ctx)
    ExtendedBind(over = over′, list = list′, label_map = n.label_map)
end

annotate_scalar(n::BindNode, ctx) =
    annotate(n, ctx)

function annotate(n::DefineNode, ctx)
    over′ = annotate(n.over, ctx)
    list′ = annotate_scalar(n.list, ctx)
    Define(over = over′, list = list′, label_map = n.label_map)
end

annotate(n::FromNode, ctx) =
    n

function annotate_scalar(n::FunctionNode, ctx)
    args′ = annotate_scalar(n.args, ctx)
    Fun(name = n.name, args = args′)
end

function annotate_scalar(n::GetNode, ctx)
    rebind(n.over, Get(name = n.name), ctx)
end

function annotate(n::GroupNode, ctx)
    over′ = annotate(n.over, ctx)
    by′ = annotate_scalar(n.by, ctx)
    Group(over = over′, by = by′, label_map = n.label_map)
end

function annotate(n::HighlightNode, ctx)
    over′ = annotate(n.over, ctx)
    Highlight(over = over′, color = n.color)
end

function annotate_scalar(n::HighlightNode, ctx)
    over′ = annotate_scalar(n.over, ctx)
    Highlight(over = over′, color = n.color)
end

function annotate(n::JoinNode, ctx)
    over′ = annotate(n.over, ctx)
    joinee′ = annotate(n.joinee, ctx)
    on′ = annotate_scalar(n.on, ctx)
    ExtendedJoin(over = over′, joinee = joinee′, on = on′, left = n.left, right = n.right)
end

function annotate(n::LimitNode, ctx)
    over′ = annotate(n.over, ctx)
    Limit(over = over′, offset = n.offset, limit = n.limit)
end

annotate_scalar(n::LiteralNode, ctx) =
    n

function annotate(n::OrderNode, ctx)
    over′ = annotate(n.over, ctx)
    by′ = annotate_scalar(n.by, ctx)
    Order(over = over′, by = by′)
end

function annotate(n::PartitionNode, ctx)
    over′ = annotate(n.over, ctx)
    by′ = annotate_scalar(n.by, ctx)
    order_by′ = annotate_scalar(n.order_by, ctx)
    Partition(over = over′, by = by′, order_by = order_by′, frame = n.frame)
end

function annotate(n::SelectNode, ctx)
    over′ = annotate(n.over, ctx)
    list′ = annotate_scalar(n.list, ctx)
    Select(over = over′, list = list′, label_map = n.label_map)
end

function annotate_scalar(n::SortNode, ctx)
    over′ = annotate_scalar(n.over, ctx)
    Sort(over = over′, value = n.value, nulls = n.nulls)
end

annotate_scalar(n::VariableNode, ctx) =
    n

function annotate(n::WhereNode, ctx)
    over′ = annotate(n.over, ctx)
    condition′ = annotate_scalar(n.condition, ctx)
    Where(over = over′, condition = condition′)
end

# Type resolution.

function resolve!(ctx::AnnotateContext)
    for box in ctx.boxes
        over = box.over
        if over !== nothing
            h = get_handle(ctx, over)
            t = resolve(over[])
            t = add_handle(t, h)
            box.handle = h
            box.type = t
        end
    end
end

function resolve(n::AppendNode)
    t = box_type(n.over)
    for m in n.list
        t = intersect(t, box_type(m))
    end
    t
end

function resolve(n::AsNode)
    t = box_type(n.over)
    fields = FieldTypeMap(n.name => t.row)
    row = RowType(fields)
    BoxType(n.name, row, t.handle_map)
end

function resolve(n::DefineNode)
    t = box_type(n.over)
    fields = FieldTypeMap()
    for (f, ft) in t.row.fields
        if f in keys(n.label_map)
            ft = ScalarType()
        end
        fields[f] = ft
    end
    for f in keys(n.label_map)
        if !haskey(fields, f)
            fields[f] = ScalarType()
        end
    end
    row = RowType(fields, t.row.group)
    BoxType(t.name, row, t.handle_map)
end

resolve(n::Union{ExtendedBindNode, HighlightNode, LimitNode, OrderNode, WhereNode}) =
    box_type(n.over)

function resolve(n::ExtendedJoinNode)
    lt = box_type(n.over)
    rt = box_type(n.joinee)
    t = union(lt, rt)
    n.type = t
    t
end

function resolve(n::FromNode)
    fields = FieldTypeMap()
    for f in n.table.columns
        fields[f] = ScalarType()
    end
    row = RowType(fields)
    BoxType(n.table.name, row)
end

function resolve(n::GroupNode)
    t = box_type(n.over)
    fields = FieldTypeMap()
    for name in keys(n.label_map)
        fields[name] = ScalarType()
    end
    row = RowType(fields, t.row)
    BoxType(t.name, row)
end

function resolve(n::PartitionNode)
    t = box_type(n.over)
    row = RowType(t.row.fields, t.row)
    BoxType(t.name, row, t.handle_map)
end

function resolve(n::SelectNode)
    t = box_type(n.over)
    fields = FieldTypeMap()
    for name in keys(n.label_map)
        fields[name] = ScalarType()
    end
    row = RowType(fields)
    BoxType(t.name, row)
end

# Collecting references.

gather!(refs::Vector{SQLNode}, n::SQLNode) =
    gather!(refs, n[])

function gather!(refs::Vector{SQLNode}, ns::Vector{SQLNode})
    for n in ns
        gather!(refs, n)
    end
end

gather!(refs::Vector{SQLNode}, ::Union{AbstractSQLNode, Nothing}) =
    nothing

gather!(refs::Vector{SQLNode}, n::Union{AsNode, BoxNode, HighlightNode, SortNode}) =
    gather!(refs, n.over)

function gather!(refs::Vector{SQLNode}, n::ExtendedBindNode)
    gather!(refs, n.over)
    gather!(refs, n.list)
    n.owned = true
end

gather!(refs::Vector{SQLNode}, n::FunctionNode) =
    gather!(refs, n.args)

function gather!(refs::Vector{SQLNode}, n::Union{AggregateNode, GetNode, HandleBoundNode, NameBoundNode})
    push!(refs, n)
end

# Validating references.

function validate(t::BoxType, ref::SQLNode, ctx)
    if @dissect ref over |> HandleBound(handle = handle)
        if handle in keys(t.handle_map)
            ht = t.handle_map[handle]
            if ht isa AmbiguousType
                throw(ReferenceError(REFERENCE_ERROR_TYPE.AMBIGUOUS_HANDLE,
                                     path = get_path(ctx, ref)))
            end
            validate(ht, over, ctx)
        else
            throw(ReferenceError(REFERENCE_ERROR_TYPE.UNDEFINED_HANDLE,
                                 path = get_path(ctx, ref)))
        end
    else
        validate(t.row, ref, ctx)
    end
end

function validate(t::RowType, ref::SQLNode, ctx)
    while @dissect ref over |> NameBound(name = name)
        ft = get(t.fields, name, EmptyType())
        if !(ft isa RowType)
            type =
                ft isa EmptyType ? REFERENCE_ERROR_TYPE.UNDEFINED_NAME :
                ft isa ScalarType ? REFERENCE_ERROR_TYPE.UNEXPECTED_SCALAR_TYPE :
                ft isa AmbiguousType ? REFERENCE_ERROR_TYPE.AMBIGUOUS_NAME : error()
            throw(ReferenceError(type, name = name, path = get_path(ctx, ref)))
        end
        t = ft
        ref = over
    end
    if @dissect ref nothing |> Get(name = name)
        ft = get(t.fields, name, EmptyType())
        if !(ft isa ScalarType)
            type =
                ft isa EmptyType ? REFERENCE_ERROR_TYPE.UNDEFINED_NAME :
                ft isa RowType ? REFERENCE_ERROR_TYPE.UNEXPECTED_ROW_TYPE :
                ft isa AmbiguousType ? REFERENCE_ERROR_TYPE.AMBIGUOUS_NAME : error()
            throw(ReferenceError(type, name = name, path = get_path(ctx, ref)))
        end
    elseif @dissect ref nothing |> Agg(name = name)
        if !(t.group isa RowType)
            type =
                t.group isa EmptyType ? REFERENCE_ERROR_TYPE.UNEXPECTED_AGGREGATE :
                t.group isa AmbiguousType ? REFERENCE_ERROR_TYPE.AMBIGUOUS_AGGREGATE : error()
            throw(ReferenceError(type, path = get_path(ctx, ref)))
        end
    else
        error()
    end
end

function gather_and_validate!(refs::Vector{SQLNode}, n, t::BoxType, ctx)
    start = length(refs) + 1
    gather!(refs, n)
    for k in start:length(refs)
        validate(t, refs[k], ctx)
    end
end

function route(lt::BoxType, rt::BoxType, ref::SQLNode)
    if @dissect ref over |> HandleBound(handle = handle)
        if get(lt.handle_map, handle, EmptyType()) isa EmptyType
            return 1
        else
            return -1
        end
    end
    return route(lt.row, rt.row, ref)
end

function route(lt::RowType, rt::RowType, ref::SQLNode)
    while @dissect ref over |> NameBound(name = name)
        lt′ = get(lt.fields, name, EmptyType())
        if lt′ isa EmptyType
            return 1
        end
        rt′ = get(rt.fields, name, EmptyType())
        if rt′ isa EmptyType
            return -1
        end
        @assert lt′ isa RowType && rt′ isa RowType
        lt = lt′
        rt = rt′
        ref = over
    end
    if @dissect ref Get(name = name)
        if name in keys(lt.fields)
            return -1
        else
            return 1
        end
    elseif @dissect ref over |> Agg(name = name)
        if lt.group isa RowType
            return -1
        else
            return 1
        end
    else
        error()
    end
end

# Linking references through box nodes.

function link!(ctx::AnnotateContext)
    root_box = ctx.boxes[end]
    for (f, ft) in root_box.type.row.fields
        if ft isa ScalarType
            push!(root_box.refs, Get(f))
        end
    end
    for box in reverse(ctx.boxes)
        box.over !== nothing || continue
        refs′ = SQLNode[]
        for ref in box.refs
            if (@dissect ref over |> HandleBound(handle = handle)) && handle == box.handle
                push!(refs′, over)
            else
                push!(refs′, ref)
            end
        end
        link!(box.over[], refs′, ctx)
    end
end

function link!(n::AppendNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
    for l in n.list
        box = l[]::BoxNode
        append!(box.refs, refs)
    end
end

function link!(n::AsNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    for ref in refs
        if @dissect ref over |> NameBound(name = name)
            @assert name == n.name
            push!(box.refs, over)
        elseif @dissect ref HandleBound()
            push!(box.refs, ref)
        else
            error()
        end
    end
end

function link!(n::DefineNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    seen = Set{Symbol}()
    for ref in refs
        if (@dissect ref (nothing |> Get(name = name))) && name in keys(n.label_map)
            !(name in seen) || continue
            push!(seen, name)
            col = n.list[n.label_map[name]]
            gather_and_validate!(box.refs, col, box.type, ctx)
        else
            push!(box.refs, ref)
        end
    end
end

function link!(n::ExtendedBindNode, refs::Vector{SQLNode}, ctx)
    if !n.owned
        gather_and_validate!(SQLNode[], n.list, EMPTY_BOX, ctx)
    end
    box = n.over[]::BoxNode
    append!(box.refs, refs)
end

function link!(n::ExtendedJoinNode, refs::Vector{SQLNode}, ctx)
    lbox = n.over[]::BoxNode
    rbox = n.joinee[]::BoxNode
    gather_and_validate!(n.lateral, n.joinee, lbox.type, ctx)
    append!(lbox.refs, n.lateral)
    refs′ = SQLNode[]
    gather_and_validate!(refs′, n.on, n.type, ctx)
    append!(refs′, refs)
    for ref in refs′
        turn = route(lbox.type, rbox.type, ref)
        if turn < 0
            push!(lbox.refs, ref)
        else
            push!(rbox.refs, ref)
        end
    end
end

link!(::FromNode, ::Vector{SQLNode}, ctx) =
    nothing

function link!(n::GroupNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    gather_and_validate!(box.refs, n.by, box.type, ctx)
    for ref in refs
        if @dissect ref nothing |> Agg(args = args, filter = filter)
            gather_and_validate!(box.refs, args, box.type, ctx)
            if filter !== nothing
                gather_and_validate!(box.refs, filter, box.type, ctx)
            end
        end
    end
end

function link!(n::Union{HighlightNode, LimitNode}, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
end

function link!(n::OrderNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
    gather_and_validate!(box.refs, n.by, box.type, ctx)
end

function link!(n::PartitionNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    for ref in refs
        if @dissect ref nothing |> Agg(args = args, filter = filter)
            gather_and_validate!(box.refs, args, box.type, ctx)
            if filter !== nothing
                gather_and_validate!(box.refs, filter, box.type, ctx)
            end
        else
            push!(box.refs, ref)
        end
    end
    gather_and_validate!(box.refs, n.by, box.type, ctx)
    gather_and_validate!(box.refs, n.order_by, box.type, ctx)
end

function link!(n::SelectNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    gather_and_validate!(box.refs, n.list, box.type, ctx)
end

function link!(n::WhereNode, refs::Vector{SQLNode}, ctx)
    box = n.over[]::BoxNode
    append!(box.refs, refs)
    gather_and_validate!(box.refs, n.condition, box.type, ctx)
end


text:
type USD <: StepData

	A # Current sum of gradients. Start with a zeros sparse array
	f # Current array of gradient functions. Start with zero.
	  # Has all gradient function in it, no chunking here
 	y # Current array of sparse gradient sums. Initialize empty.
	  # One gradient sum per chunk
	I # Current number of total chunks

	s::Function # Number of chunks to resample
	u::Function # Number of new chunks to sample
	beta::Function
	getStep::Function
	stepString::AbstractString
	chunkSize::Int64
	numVars::Int64
	sagFsagaT::Bool

	function USD(s::Function, u::Function, beta::Function, numVars::Int64, stepString::AbstractString, chunkSize::Int64,sagFsagaT::Bool)
		new(spzeros(numVars,1),Function[], SparseMatrixCSC[], 0, s, u, beta, uComputation, stepString, chunkSize, numVars, sagFsagaT)
	end
end

function onlyAddBeta1(u::Function, uString::AbstractString, numVars::Int64, chunkSize::Int64,sagFsagaT::Bool)

	s = (k,I)-> 0
	beta = (k)->1

	stepString  = "uAb1."*uString*".cs="*string(chunkSize)*"."*string(sagFsagaT)

	USD(s, u , beta, numVars, stepString, chunkSize,sagFsagaT)

end

function onlyUpdateBeta1(s::Function, sString::AbstractString, numVars::Int64, ntp::Int64, chunkSize::Int64,sagFsagaT::Bool)

	u = (k,I)-> round(Int,floor(k==0 ? int(floor(ntp/chunkSize)) : 0))
	beta = (k)->1

	stepString  = "uUb1."*sString*".cs="*string(chunkSize)*"."*string(sagFsagaT)

	USD(s, u , beta, numVars, stepString, chunkSize,sagFsagaT)
end

function uLinBeta1(c::Float64,numVars::Int64, chunkSize::Int64,sagFsagaT::Bool)

	if c<1
		error("C must be >=1")
	end
	s = (k,I)-> round(Int,floor(k==0 ? 0 : c))
	u = (k,I)-> round(Int,c)
	beta = (k)->1

	stepString  = "uLb1"*".c="*string(c)*".cs="*string(chunkSize)*"."*string(sagFsagaT)

	USD(s, u , beta, numVars, stepString, chunkSize,sagFsagaT)
end

function uQuadBeta1(c::Float64,numVars::Int64, chunkSize::Int64,sagFsagaT::Bool)

	s = (k,I)-> round(Int,ceil(k==0 ? 0 : c*k))
	u = (k,I)-> round(Int,ceil(c*(k+1)))
	beta = (k)->1

	stepString  = "uQb1"*".c="*string(c)*".cs="*string(chunkSize)*"."*string(sagFsagaT)

	USD(s, u , beta, numVars, stepString, chunkSize,sagFsagaT)
end

function uExpBeta1(r::Float64, numVars::Int64, chunkSize::Int64,sagFsagaT::Bool)

	c = 1.0 / (r-1)
	s = (k,I)-> round(Int,ceil(k==0 ? 0 : c*(r/(r-1))^(k-1)))
	u = (k,I)-> round(Int,ceil(k==0 ? c*(r-1) : c*(r/(r-1))^(k-1)))
	beta = (k)->1

	stepString  = "uEb1"*".r="*string(r)*".cs="*string(chunkSize)*"."*string(sagFsagaT)

	USD(s, u , beta, numVars, stepString, chunkSize,sagFsagaT)
end

function uExp2Beta1(r::Float64, numVars::Int64, chunkSize::Int64,sagFsagaT::Bool)

	c = 1.0
	s = (k,I)-> round(Int,ceil(k==0 ? 0 : c*(r/(r-1))^(k-1)))
	u = (k,I)-> round(Int,ceil(k==0 ? c*(r-1) : c*(r/(r-1))^(k-1)))
	beta = (k)->1

	stepString  = "uE2b1"*".r="*string(r)*".cs="*string(chunkSize)*"."*string(sagFsagaT)

	USD(s, u , beta, numVars, stepString, chunkSize,sagFsagaT)
end

function uComputation(x, k, gnum, sd::USD, problem::Problem; outputLevel = 0)

	outputLevel > 0 && println("Starting uComputation")

	U = (sd.I+1):(sd.I+sd.u(k,sd.I))

	outputLevel > 0  && println("sd.I=$(sd.I) sd.u(k,sd.I)=$(sd.u(k,sd.I)) sd.s(k,sd.I)=$(sd.s(k,sd.I))")

	for i in U
		try
			for j in 1:sd.chunkSize
				push!(sd.f,consume(problem.getNextSampleFunction))
			end
		catch
			error("Not producing anymore! Probably out of datapoints.")
		end
	end

	# S is composed of chunk indices
	S = StatsBase.sample(1:sd.I,sd.s(k,sd.I);replace=false)

	# B is the current sum of old gradients (for S)
	B = spzeros(sd.numVars,1)
	for i in collect(S)
		B +=sd.y[i]
	end

	# sumy is the new sum of gradients (for both U and S)
	sumy = spzeros(sd.numVars,1)

	for i in collect(U)
		thisSum = spzeros(sd.numVars,1)
		for j in (i-1)*sd.chunkSize+1:i*sd.chunkSize
			(f,sampleG) = (sd.f[j])(x)
			# MIGHT BE SLOW
			thisSum+=sparse(sampleG)
		end
		push!(sd.y,thisSum)
		sumy += thisSum
	end

	for i in collect(S)
		thisSum = spzeros(sd.numVars,1)
		for j in (i-1)*sd.chunkSize+1:i*sd.chunkSize
			(f,sampleG) = (sd.f[j])(x)
			# MIGHT BE SLOW
			thisSum+=sparse(sampleG)
		end
		sd.y[i]=thisSum
		sumy += thisSum
	end

	gnum += sd.chunkSize * (sd.s(k,sd.I)+sd.u(k,sd.I))

	if sd.sagFsagaT == false
		g = (sd.A - B + sumy )/(sd.chunkSize*(sd.I+sd.u(k,sd.I)))
	else
		g = ((sd.s(k,sd.I) > 0  ? sd.s(k,sd.I)*((sd.beta(k)/sd.I)*sd.A) - sd.beta(k)*B : 0) + sumy )/(sd.chunkSize*(sd.s(k,sd.I)+sd.u(k,sd.I)))
	end

	sd.I = sd.I + sd.u(k,sd.I)

	sd.A=sd.A-B+sumy

	(g, gnum)
end


text:
__precompile__(true)

module SharedDataStructures

using Semaphores
using Base: shm_mmap_array, shm_unlink, JL_O_CREAT, JL_O_RDWR
import Base: hash, length, eltype, empty!, delete!, isempty, push!, pop!, unshift!, shift!, getindex, start, next, done, show, close, endof, getindex, setindex!, haskey

include("sysv_shm.jl")
include("shared_circdq.jl")
include("shm_dict.jl")

export SharedCircularDeque, ShmDict
export hash, length, eltype, empty!, delete!, isempty, push!, pop!, unshift!, shift!, getindex, start, next, done, show, capacity, close

end # module