text:
#=
Code related with input output (IO) of .nc files directly to/from ClimArrays
An initial version of parts of this code was taken from:
https://github.com/rafaqz/GeoData.jl
=#
using NCDatasets
export NCDataset
export nckeys, ncdetails, globalattr
export ncread, ncwrite

dim_to_commonname(::Lat) = "lat"
dim_to_commonname(::Lon) = "lon"
dim_to_commonname(::Time) = "time"
dim_to_commonname(::Pre) = "level"
dim_to_commonname(D::Dim) = string(DimensionalData.name(D))

#########################################################################
# Utilities
#########################################################################
"""
    nckeys(file::String)
Return all keys of the `.nc` file in `file`.
"""
function nckeys(path::String)
    NCDataset(path) do ds
        return keys(ds)
    end
end
nckeys(a::NCDataset) = keys(a)

"""
    ncdetails(file::String, io = stdout)
Print details about the `.nc` file in `file` on `io`.
"""
function ncdetails(file::String, io = stdout)
    NCDataset(file) do ds
        show(io, MIME"text/plain"(), ds)
    end
end
ncdetails(ds::NCDataset, io = stdout) = show(io, MIME"text/plain"(), ds)

"""
    globalattr(file::String) → Dict
Return the global attributes of the .nc file.
"""
function globalattr(file::String)
    NCDataset(file) do ds
        return Dict(ds.attrib)
    end
end

#########################################################################
# Reading
#########################################################################
"""
    ncread(file, var; name, kwargs...) → A
Load the variable `var` from the `file` and convert it into a [`ClimArray`](@ref)
with proper dimension mapping and also containing the variable attributes as a dictionary.
Dimension attributes are also given to the dimensions of `A`, if any exist.
See keywords below for specifications for unstructured grids.

`file` can be a string to a `.nc` file. Or, it can be an
`NCDataset`, which allows you to lazily combine different
`.nc` data (typically split by time), e.g.
```julia
using Glob # for getting all files
alldata = glob("toa_fluxes_2020_*.nc")
file = NCDataset(alldata; aggdim = "time")
A = ClimArray(file, "tow_sw_all")
```

`var` is a `String` denoting which variable to load.
For `.nc` data containing groups `var` can also be a tuple `("group_name", "var_name")`
that loads a specific variable from a specific group.
In this case, the attributes of both the group and the CF-variable are attributed to
the created [`ClimArray`](@ref).

See also [`ncdetails`](@ref), [`nckeys`](@ref) and [`ncwrite`](@ref).

We do two performance improvements while loading the data:
1. If there are no missing values in the data (according to CF standards), the
   returned array is automatically converted to a concrete type (i.e. `Union{Float32, Missing}`
   becomes `Float32`).
2. Dimensions that are ranges (i.e. sampled with constant step size) are automatically
   transformed to a standard Julia `Range` type (which makes sub-selecting faster).

## Keywords
* `name` optionally rename loaded array.
* `grid = nothing` optionally specify whether the underlying grid is `grid = LonLatGrid()`
  or `grid = UnstructuredGrid()`, see [Types of spatial coordinates](@ref).
  If `nothing`, we try to deduce automatically based on
  the names of dimensions and other keys of the `NCDataset`.
* `lon, lat`. These two keywords are useful in unstructured grid data where the grid
  information is provided in a *separate .nc file*. What we need is the user to
  provide vectors of the central longitude and central latitude of each grid point.
  This is done e.g. by
  ```julia
  ds = NCDataset("path/to/grid.nc");
  lon = Array(ds["clon"]);
  lat = Array(ds["clat"]);
  ```
  If `lon, lat` are given, `grid` is automatically assumed `UnstructuredGrid()`.
"""
function ncread(path::Union{String, Vector{String}}, args...; kwargs...)
    NCDataset(path) do ds
        data = ncread(ds, args...; kwargs...)
        return data
    end
end

# TODO: Allow this function to take as input a tuple of indices, e.g. (:, :, 1:5)
# and only load this part, and correctly and instantly make it a ClimArray, which
# can solve "large memory" or "large data" problems. This funcionality
# must be sure to load the correct ranges of dimensions as well though!
#
# TODO: Allow reading multiple variables at once. This has the performance benefit
# of not re-creating dimensions all the time.

function ncread(ds::NCDatasets.AbstractDataset, var;
        name = var2name(var), grid = nothing, lon = nothing, lat = nothing,
    )
    if lon isa Vector && lat isa Vector
        gridtype = UnstructuredGrid()
    elseif isnothing(grid)
        gridtype = autodetect_grid(ds)
    else
        gridtype = grid
    end

    if gridtype == UnstructuredGrid()
        return ncread_unstructured(ds, var, name, lon, lat)
    else
        return ncread_lonlat(ds, var, name)
    end
end

function autodetect_grid(ds)
    if haskey(ds, "reduced_points") || haskey(ds.dim, "ncells") || haskey(ds, "clon")
        return UnstructuredGrid()
    else
        return LonLatGrid()
    end
end

var2name(var) = string(var)
var2name(var::Tuple) = join(var, "_")

#########################################################################
# Reading: LonLatGrid and main reading functionality
#########################################################################
# Notice that this function properly loads even without any spatial coordinate
function ncread_lonlat(ds::NCDatasets.AbstractDataset, var, name)
    cfvar = var isa String ? ds[var] : ds.group[var[1]][var[2]]
    attrib = get_attributes_from_var(ds, cfvar, var)
    A = cfvar |> Array
    dnames = Tuple(NCDatasets.dimnames(cfvar))
    if !any(ismissing, A)
        A = nomissing(A)
    end
    data = ClimArray(A, create_dims(ds, dnames); name = Symbol(name), attrib = attrib)
    return data
end

get_attributes_from_var(ds, cfvar, var::String) = Dict(cfvar.attrib)
function get_attributes_from_var(ds, cfvar, var::Tuple)
    d1 = Dict(cfvar.attrib)
    d2 = Dict(ds.group[var[1]].attrib)
    return merge!(d2, d1)
end

"""
    create_dims(ds::NCDatasets.AbstractDataset, dnames)
Create a tuple of `Dimension`s from the `dnames` (tuple of strings).
"""
function create_dims(ds::NCDatasets.AbstractDataset, dnames)
    # true_dims = getindex.(Ref(COMMONNAMES), dnames)
    true_dims = to_proper_dimensions(dnames)
    dim_values = Array.(getindex.(Ref(ds), dnames))
    # Some stupid datasets return a union{Missing} type for dimensions.
    # this is of course nonsense, a dimension cannot have "missing" values.
    if any(d -> Missing <: eltype(d), dim_values)
        dim_values = nomissing.(dim_values)
    end
    optimal_values = vector2range.(dim_values)
    attribs = [
        ds[d].attrib isa NCDatasets.BaseAttributes ? Dict(ds[d].attrib) : nothing
        for d in dnames
    ]
    out = []
    for i in 1:length(true_dims)
        push!(out, true_dims[i](optimal_values[i]; metadata = attribs[i]))
    end
    return (out...,)
end

function to_proper_dimensions(dnames)
    r = []
    for n in dnames
        if haskey(COMMONNAMES, n)
            push!(r, COMMONNAMES[n])
        else
            @warn """
            Dimension name "$n" not in common names. Strongly recommended to ask for
            adding this name to COMMONNAMES on GitHub. Making generic dimension for now...
            """
            push!(r, Dim{Symbol(n)})
        end
    end
    return (r...,)
end

export Dim # for generic dimensions this must be exported

#########################################################################
# Making vectors → ranges
#########################################################################
function vector2range(x::Vector{<:Real})
    length(x) < 3 && return x
    dx = x[2]-x[1]
    for i in 3:length(x)
        x[i]-x[i-1] ≠ dx && return x # if no constant step, return array as is
    end
    # do not check value equality, only difference equality
    return x[1]:dx:x[end]
end

function vector2range(t::Vector{<:Dates.AbstractTime})
    tsamp = temporal_sampling(t)
    period = tsamp2period(tsamp)
    isnothing(period) && return t
    t1 = period == :hourly ? t[1] : Date(t[1])
    tf = period == :hourly ? t[end] : Date(t[end])
    r = t1:period:tf
    return r == t ? r : t # final safety check to ensure equal values
end

vector2range(r::AbstractRange) = r

#########################################################################
# Reading: UnstructuredGrid
#########################################################################
export SVector

function ncread_unstructured(ds::NCDatasets.AbstractDataset, var::String, name, lon, lat)
    cfvar = var isa String ? ds[var] : ds.group[var[1]][var[2]]
    attrib = get_attributes_from_var(ds, cfvar, var)

    # Here we generate the longitudes and latitudes based on whether we have
    # reduced points or not, and obtain the name of the coord dimension in the `ds`
    if isnothing(lon)
        lonlat, original_grid_dim = load_coordinate_points(ds)
        lonlat = convert_to_degrees(lonlat, ds)
    else
        lonlat = [SVector(lo, la) for (lo, la) in zip(lon, lat)]
        original_grid_dim = intersect(NCDatasets.dimnames(cfvar), POSSIBLE_CELL_NAMES)[1]
        lonlat = convert_to_degrees(lonlat)
    end

    # TODO: I've noticed that this converts integer dimension (like pressure)
    # into Float64, but I'm not sure why...
    alldims = [NCDatasets.dimnames(cfvar)...]

    # Set up the remaining dimensions of the dataset
    @assert original_grid_dim ∈ alldims
    i = findfirst(x -> x == original_grid_dim, alldims)
    remainingdims = deleteat!(copy(alldims), i)
    actualdims = Any[create_dims(ds, remainingdims)...]

    # Make coordinate dimension
    si = sortperm(lonlat, by = reverse)
    coords = Coord(lonlat, (Lon, Lat))
    insert!(actualdims, i, coords)

    # Create array, sort by latitude and remove missings
    A = Array(cfvar)
    X = ClimArray(A, Tuple(actualdims))
    X = X[Coord(si)]
    if !any(ismissing, X)
        X = nomissing(X)
    end
    return ClimArray(X; name = Symbol(name), attrib)
end

const POSSIBLE_CELL_NAMES = ["ncells", "cell", "rgrid", "grid"]

function has_unstructured_key(ds)
    any(x -> haskey(ds.dim, x), POSSIBLE_CELL_NAMES) ||
    haskey(ds, "lon") || haskey(ds, "clon")
end

function load_coordinate_points(ds)
    if haskey(ds, "reduced_points")
        lonlat = reduced_grid_to_points(ds["lat"], ds["reduced_points"])
        original_grid_dim = "rgrid" # TODO: Specific to CDO Gaussian grid
    elseif has_unstructured_key(ds)
        if haskey(ds, "lon")
            lons = ds["lon"] |> Array .|> wrap_lon
            lats = ds["lat"] |> Array
            original_grid_dim = NCDatasets.dimnames(ds["lon"])[1]
        elseif haskey(ds, "clon")
            lons = ds["clon"] |> Array .|> wrap_lon
            lats = ds["clat"] |> Array
            original_grid_dim = NCDatasets.dimnames(ds["clon"])[1]
        else
            error("""
            We didn't find key `"lon"` or `"clon"` that represents the longitude of each
            polygon in an unstructured grid.
            """)
        end
        lonlat = [SVector(lo, la) for (lo, la) in zip(lons, lats)]
    else
        error("""
        We didn't find any of the following keys: `"ncells", "cells", "reduced_points",
        "clon", "lon"`, at least one of which is mandatory for unstructured grid.
        If your data stores the "cell" information with a different name, please open
        an issue and let us know!
        """)
    end
    return lonlat, original_grid_dim
end

function reduced_grid_to_points(lat, reduced_points)
    lonlat = SVector{2, Float32}[]
    for (i, θ) in enumerate(lat)
        n = reduced_points[i]
        dλ = Float32(360/n)
        for j in 0:n-1
            push!(lonlat, SVector(0 + dλ*j, θ))
        end
    end
    return lonlat
end

function convert_to_degrees(lonlat, ds)
    !(haskey(ds, "clat") || haskey(ds, "lat")) && return convert_to_degrees(lonlat)
    x = haskey(ds, "clat") ? ds["clat"] : ds["lat"]
    if get(x.attrib, "units", nothing) == "radian" || !any(ll -> abs(ll[2]) > π/2, lonlat)
        lonlat = [SVector(lo*180/π, la*180/π) for (lo, la) in lonlat]
    end
    return lonlat
end
function convert_to_degrees(lonlat)
    if !any(ll -> abs(ll[2]) > π/2, lonlat)
        lonlat = [SVector(lo*180/π, la*180/π) for (lo, la) in lonlat]
    end
    return lonlat
end

#########################################################################
# Saving to .nc files
#########################################################################
const DEFAULT_ATTRIBS = Dict(
    "time" => Dict(
        "units" => "days since 0000-00-01 00:00:00",
        "standard_name" => "time"
    ),
    "lon" => Dict(
        "units" => "degrees_east",
        "standard_name" => "longitude",
        "valid_range" => Float32[-180.0, 360.0]
    ),
    "lat" => Dict(
        "units" => "degrees_north",
        "standard_name" => "latitude",
        "valid_range" => Float32[-90.0, 90.0]
    ),
    "level" => Dict(
        "units" => "millibars",
        "long_name" => "pressure_level",
    ),
)

"""
    ncwrite(file::String, Xs; globalattr = Dict())
Write the given `ClimArray` instances (any iterable of `ClimArray`s or a single `ClimArray`)
to a `.nc` file following CF standard conventions using NCDatasets.jl.
Optionally specify global attributes for the `.nc` file.

The metadata of the arrays in `Xs`, as well as their dimensions, are properly written
in the `.nc` file and any necessary type convertions happen automatically.

**WARNING**: We assume that any dimensions shared between the `Xs` are identical.

See also [`ncread`](@ref).
"""
function ncwrite(file::String, X::ClimArray; globalattr = Dict())
    ncwrite(file, (X,); globalattr)
end
function ncwrite(file::String, Xs; globalattr = Dict())

    # TODO: Fixing this is very easy. Simply make a `"ncells"` dimension, and then write
    # the `"lon"` and `"lat"` cfvariables to the nc file by decomposing the coordinates
    # into longitude and latitude.
    if any(X -> hasdim(X, Coord), Xs)
        error("""
        Outputing `UnstructuredGrid` coordinates to .nc files is not yet supported,
        but it is an easy fix, see source of `ncwrite`.
        """)
    end

    ds = NCDataset(file, "c"; attrib = globalattr)
    # NCDataset("file.nc", "c"; attrib = globalattr) do ds
        for (i, X) in enumerate(Xs)
            n = string(X.name)
            if n == ""
                n = "x$i"
                @warn "$i-th ClimArray has no name, naming it $(n) instead."
            end
            println("processing variable $(n)...")
            add_dims_to_ncfile!(ds, dims(X))
            attrib = X.attrib
            isnothing(attrib) && (attrib = Dict())
            dnames = dim_to_commonname.(dims(X))
            data = Array(X)
            defVar(ds, n, data, (dnames...,); attrib)
        end
        close(ds)
    # end
end

function add_dims_to_ncfile!(ds::NCDatasets.AbstractDataset, dimensions::Tuple)
    dnames = dim_to_commonname.(dimensions)
    for (i, d) ∈ enumerate(dnames)
        haskey(ds, d) && continue
        println("writing dimension $d...")
        v = dimensions[i].val
        # this conversion to DateTime is necessary because CFTime.jl doesn't support Date
        eltype(v) == Date && (v = DateTime.(v))
        l = length(v)
        defDim(ds, d, l) # add dimension entry
        attrib = dimensions[i].metadata
        if (isnothing(attrib) || attrib == DimensionalData.NoMetadata()) && haskey(DEFAULT_ATTRIBS, d)
            @warn "Dimension $d has no attributes, adding default attributes (mandatory)."
            attrib = DEFAULT_ATTRIBS[d]
        end
        # write dimension values as a variable as well (mandatory)
        defVar(ds, d, v, (d, ); attrib = attrib)
    end
end


text:
module Deform

using JuLIP
using JuLIP.Potentials
using Main.Fast_Hertz
using Main.Minimise_Energy
using Main.Granular

import JuLIP: energy, forces, hessian_pos

export deform_Config

function deform_Config(at::AbstractAtoms, deformType::AbstractString, deformValue::AbstractFloat)
    Lx = get_data(at, "Length")
    Ly = get_data(at, "Height")
    N = get_data(at, "Num")
    ratio = get_data(at, "ratio")
    x, y, _ = xyz(at)
    r = at.M
    if deformType == "phi"
        ϕ = deformValue
        set_data!(at, "Phi", ϕ)
        ρ = sqrt(2.0*ϕ*Lx*Ly/N/π/(1 + ratio^2))
        at.M = ρ*(0.4*at.Z + 0.6)
    elseif deformType == "grow"
        grow_ratio = deformValue
        at.M = at.M*grow_ratio
        ϕ = π*sum(at.M.^2)/Lx/Ly
        set_data!(at, "Phi", ϕ)
    elseif deformType == "strain"
        γ = deformValue
        strainType = get_data(at, "StrainType")
        strain = get_data(at, "Strain")
        strain = strain + γ
        set_data!(at, "Strain", strain)
        if strainType == "SimpleShear"
            x = x + γ*y
            mycell = [Lx strain*Ly 0.0; 0.0 Ly 0.0; 0.0 0.0 1.0]
        elseif strainType == "PureShear"
            γ_temp = γ/(1 + strain)
            Lx = Lx * (1 + γ_temp)
            Ly = Ly / (1 + γ_temp)
            set_data!(at, "Length", Lx)
            set_data!(at, "Height", Ly)
            x = x * (1 + γ_temp)
            y = y / (1 + γ_temp)
        end
        set_cell!(at, mycell)
        X = vecs([x; y; zeros(N)'])
        at.X = X
    end

    return Main.Minimise_Energy.min_Energy(at)
end

end

ext: jl


text:
module maagar

end #module

max_stars_repo_name: mosheduminer/maagar

max_issues_repo_name: mosheduminer/maagar

max_forks_repo_name: mosheduminer/maagar


text:
@testset "Fragment validation" begin
    str = """
    query myquery{field{...TopLevel}}

    query myquery2{field{... on User {...forInlineFragment}}}

    fragment TopLevel on Name {...SubTopLevel}

    fragment SubTopLevel on Name {field}

    fragment forInlineFragment on Name {field}
    """

    frags = GraphQLParser.find_fragments(GraphQLParser.parse(str))
    for name in ("forInlineFragment", "TopLevel", "SubTopLevel")
        @test haskey(frags, name)
    end

    str = """
    query myquery{field{...TopLevel}}
    """
    errors = GraphQLParser.validate(GraphQLParser.parse(str))
    @test only(errors) isa GraphQLParser.UnknownFragment
    @test only(errors[1].locations) == Loc(1, 24)

    str = """
    query myquery2{field{... on User {...TopLevel}}}
    """
    errors = GraphQLParser.validate(GraphQLParser.parse(str))
    @test only(errors) isa GraphQLParser.UnknownFragment
    @test errors[1].locations[1] == Loc(1, 38)

    # Counter Example 144
    str = """
    {
        dog {
            ...undefinedFragment
        }
    }
    """
    errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
    @test only(errors) isa GraphQLParser.UnknownFragment
    @test only(errors[1].locations) == Loc(3, 12)

    str = """
    query myquery1{field{...TopLevel}}

    query myquery2{field{...TopLevel2}}

    query myquery3{field{... on User {...forInlineFragment}}}

    fragment TopLevel2 on Name {...SubTopLevel}
    """
    errors = GraphQLParser.validate(GraphQLParser.parse(str))

    @test length(errors) == 3
    @test all(err -> isa(err, GraphQLParser.UnknownFragment), errors)

    # Example 137
    str = """
    {
        dog {
            ...fragmentOne
            ...fragmentTwo
        }
    }

    fragment fragmentOne on Dog {
        name
    }

    fragment fragmentTwo on Dog {
        owner {
            name
        }
    }
    """
    errors = GraphQLParser.validate(GraphQLParser.parse(str))
    @test isempty(errors)

    # Counter Example 138
    str = """
    {
        dog {
            ...fragmentOne
        }
    }

    fragment fragmentOne on Dog {
        name
    }

    fragment fragmentOne on Dog {
        owner {
            name
        }
    }
    """
    errors = GraphQLParser.validate(GraphQLParser.parse(str))
    @test length(errors) == 2
    @test errors[1] isa GraphQLParser.RepeatedFragmentDefinition
    @test errors[2] isa GraphQLParser.RepeatedFragmentDefinition

    # Counter Example 143
    str = """
    fragment nameFragment on Dog { # unused
        name
    }

    {
        dog {
            name
        }
    }
    """
    errors = GraphQLParser.validate(GraphQLParser.parse(str))
    @test only(errors) isa GraphQLParser.UnusedFragment
end

@testset "Operations" begin
    @testset "5.2.1.1 Operation Name Uniqueness" begin
        # Example 103
        str = """
        query getDogName {
            dog {
                name
            }
        }

        query getOwnerName {
            dog {
                owner {
                    name
                }
            }
        }
        """
        errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
        @test isempty(errors)

        # Counter Example 104
        str = """
        query getName {
            dog {
                name
            }
        }

        query getName {
            dog {
                owner {
                    name
                }
            }
        }
        """
        errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
        @test length(errors) == 2
        @test errors[1] isa GraphQLParser.RepeatedOperationDefinition
        @test errors[2] isa GraphQLParser.RepeatedOperationDefinition
        @test errors[1].locations == [Loc(1,1)]
        @test errors[2].locations == [Loc(7,1)]

        # Counter Example 105
        str = """
        query dogOperation {
            dog {
                name
            }
        }

        mutation dogOperation {
            mutateDog {
                id
            }
        }
        """
        errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
        @test length(errors) == 2
        @test errors[1] isa GraphQLParser.RepeatedOperationDefinition
        @test errors[2] isa GraphQLParser.RepeatedOperationDefinition
        @test errors[1].locations == [Loc(1,1)]
        @test errors[2].locations == [Loc(7,1)]
    end

    @testset "5.2.2.1 Lone Anonymous Operation" begin
        # Example 106
        str = """
        {
            dog {
                name
            }
        }
        """
        errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
        @test isempty(errors)

        # Counter Example 107
        str = """
        {
            dog {
                name
            }
        }

        query getName {
            dog {
                owner {
                    name
                }
            }
        }
        """
        errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
        @test only(errors) isa GraphQLParser.AnonymousOperationNotAlone
        @test errors[1].locations[1] == Loc(1,1)
    end
end

@testset "Variables" begin
    # Counter Example 165
    str = """
    query houseTrainedQuery(\$atOtherHomes: Boolean, \$atOtherHomes: Boolean) {
        dog {
            isHouseTrained(atOtherHomes: \$atOtherHomes)
        }
    }
    """
    errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
    @test length(errors) == 2
    @test errors[1] isa GraphQLParser.RepeatedVariableDefinition
    @test errors[2] isa GraphQLParser.RepeatedVariableDefinition

    # Example 166
    str = """
    query A(\$atOtherHomes: Boolean) {
        ...HouseTrainedFragment
    }

    query B(\$atOtherHomes: Boolean) {
        ...HouseTrainedFragment
    }

    fragment HouseTrainedFragment on Query {
        dog {
            isHouseTrained(atOtherHomes: \$atOtherHomes)
        }
    }
    """
    @test GraphQLParser.is_valid_executable_document(str)

    # Example 170
    str = """
    query variableIsDefined(\$atOtherHomes: Boolean) {
        dog {
            isHouseTrained(atOtherHomes: \$atOtherHomes)
        }
    }
    """
    @test GraphQLParser.is_valid_executable_document(str)

    # Counter Example 171
    str = """
    query variableIsNotDefined {
        dog {
            isHouseTrained(atOtherHomes: \$atOtherHomes)
        }
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnknownVariable

    # Example 172
    str = """
    query variableIsDefinedUsedInSingleFragment(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    @test GraphQLParser.is_valid_executable_document(str)

    # Counter Example 173
    str = """
    query variableIsNotDefinedUsedInSingleFragment {
        dog {
            ...isHouseTrainedFragment
        }
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnknownVariable

    # Counter Example 174
    str = """
    query variableIsNotDefinedUsedInNestedFragment {
        dog {
            ...outerHouseTrainedFragment
        }
    }

    fragment outerHouseTrainedFragment on Dog {
        ...isHouseTrainedFragment
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnknownVariable

    # Example 175
    str = """
    query houseTrainedQueryOne(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    query houseTrainedQueryTwo(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    @test GraphQLParser.is_valid_executable_document(str)

    # Counter Example 176
    str = """
    query houseTrainedQueryOne(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    query houseTrainedQueryTwoNotDefined {
        dog {
            ...isHouseTrainedFragment
        }
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnknownVariable

    # Counter Example 177
    str = """
        query variableUnused(\$atOtherHomes: Boolean) {
        dog {
            isHouseTrained
        }
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnusedVariable

    # Example 178
    str = """
    query variableUsedInFragment(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    @test GraphQLParser.is_valid_executable_document(str)

    # Counter Example 179
    str = """
    query variableNotUsedWithinFragment(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedWithoutVariableFragment
        }
    }

    fragment isHouseTrainedWithoutVariableFragment on Dog {
        isHouseTrained
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnusedVariable

    # Counter Example № 180
    str = """
    query queryWithUsedVar(\$atOtherHomes: Boolean) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    query queryWithExtraVar(\$atOtherHomes: Boolean, \$extra: Int) {
        dog {
            ...isHouseTrainedFragment
        }
    }

    fragment isHouseTrainedFragment on Dog {
        isHouseTrained(atOtherHomes: \$atOtherHomes)
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.UnusedVariable
end

@testset "Directives" begin
    # Counter Example 163
    str = """
    query (\$foo: Boolean = true, \$bar: Boolean = false) {
        field @skip(if: \$foo) @skip(if: \$bar)
    }
    """
    doc = GraphQLParser.parse(str)
    errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
    @test only(errors) isa GraphQLParser.RepeatedDirectiveName

    # Example 164
    str = """
    query (\$foo: Boolean = true, \$bar: Boolean = false) {
        field @skip(if: \$foo) {
            subfieldA
        }
        field @skip(if: \$bar) {
            subfieldB
        }
    }
    """
    doc = GraphQLParser.parse(str)
    errors = GraphQLParser.validate_operations!(GraphQLParser.ValidationError[], GraphQLParser.parse(str))
    @test isempty(errors)
end

@testset "Input Objects" begin
    # Counter Example № 161
    str = """ {
        field(arg: { field: true, field: false })
    }
    """
    errors = GraphQLParser.validate_executable_document(str)
    @test only(errors) isa GraphQLParser.RepeatedInputObjectField
end


text:
using WaterLily
using BenchmarkTools
using Profile

function body_sim(radius = 8; Re = 250, T=Float64)
    body = AutoBody((x,t)-> √sum(abs2,x .- 2radius) - radius)
    return Simulation((6radius+2,4radius+2),[1.,0.],radius; body, ν=radius/Re,T)
end
function vort_sim(L = 2^5; Re = 250, T=Float64)
    function uλ(i,vx)
        x,y,z = @. (vx-1.5)*π/L              # scaled coordinates
        i==1 && return -sin(x)*cos(y)*cos(z) # u_x
        i==2 && return  cos(x)*sin(y)*cos(z) # u_y
        return 0.                            # u_z
    end
    return Simulation((L+2,L+2,L+2),zeros(3),L;uλ,ν=L/Re,U=1.,T)
end

benchmark_step() = @benchmark mom_step!($sim.flow,$sim.pois)
benchmark_vort() = @benchmark sim_step!(sim,10) setup=(sim=vort_sim()) seconds=30
benchmark_body() = @benchmark sim_step!(sim,10) setup=(sim=body_sim())


text:
using Pkg
Pkg.activate("./GraphBLASJulia/")

include("../../src/GraphBLAS_bfs.jl")
include("../../src/GraphBLAS_sparse_square_matrix.jl")

@testset "GraphBLASJulia" begin

### algorithm correctness test
    GrB_init(GrB_NONBLOCKING);

    dim = 6;

    A = GrB_Matrix(Int64, dim, dim)
    GrB_Matrix_setElement(A, 1, ZeroBasedIndex(1), ZeroBasedIndex(0))
    GrB_Matrix_setElement(A, 2, ZeroBasedIndex(3), ZeroBasedIndex(1))
    GrB_Matrix_setElement(A, 3, ZeroBasedIndex(0), ZeroBasedIndex(2))
    GrB_Matrix_setElement(A, 4, ZeroBasedIndex(3), ZeroBasedIndex(2))
    GrB_Matrix_setElement(A, 5, ZeroBasedIndex(2), ZeroBasedIndex(3))
    GrB_Matrix_setElement(A, 6, ZeroBasedIndex(3), ZeroBasedIndex(4))

    v = GrB_Vector(Int64, dim)

    @test GraphBLAS_bfs!(A, v, ZeroBasedIndex(0)) == (GrB_SUCCESS)

    SuiteSparseGraphBLAS.dropzeros!(v);

    sol = GrB_Vector(Int64, dim)
    GrB_Vector_setElement(sol, 1, ZeroBasedIndex(0))
    GrB_Vector_setElement(sol, 4, ZeroBasedIndex(1))
    GrB_Vector_setElement(sol, 2, ZeroBasedIndex(2))
    GrB_Vector_setElement(sol, 3, ZeroBasedIndex(3))
    GrB_Vector_setElement(sol, 4, ZeroBasedIndex(4))
    GrB_wait()

    @test v == sol

### exception throw tests

    # Input matrix and vector have different dimensions
    w = GrB_Vector(Int64, dim+1)
    @test_throws DimensionMismatch GraphBLAS_bfs!(A, w, ZeroBasedIndex(0))

    # Input matrix is not square
    B = GrB_Matrix(Int64, dim, dim+1)
    @test_throws DimensionMismatch GraphBLAS_bfs!(B, v, ZeroBasedIndex(0))

    # Input vector elements' type U is not an Integer subtype
    z = GrB_Vector(Float64, dim)
    @test_throws MethodError GraphBLAS_bfs!(A, z, ZeroBasedIndex(0))

    # Input matrix is not a GrB_Matrix{T}
    @test_throws MethodError GraphBLAS_bfs!(v, v, ZeroBasedIndex(0))

    # Input vector is not a GrB_Vector{U}
    @test_throws MethodError GraphBLAS_bfs!(A, A, ZeroBasedIndex(0))

    # source is not a ZeroBasedIndex
    @test_throws MethodError GraphBLAS_bfs!(A, v, 0)

    GrB_Matrix_free(A)
    GrB_Matrix_free(B)

    GrB_Vector_free(v)
    GrB_Vector_free(sol)
    GrB_Vector_free(w)
    GrB_Vector_free(z)

    GrB_finalize()

end


text:
function YUMAParser(arg0::JString)
    return YUMAParser((JString,), arg0)
end

function YUMAParser(arg0::JString, arg1::DataProvidersManager, arg2::TimeScales)
    return YUMAParser((JString, DataProvidersManager, TimeScales), arg0, arg1, arg2)
end

function get_almanacs(obj::YUMAParser)
    return jcall(obj, "getAlmanacs", List, ())
end

function get_prn_numbers(obj::YUMAParser)
    return jcall(obj, "getPRNNumbers", List, ())
end

function get_supported_names(obj::YUMAParser)
    return jcall(obj, "getSupportedNames", JString, ())
end

function load_data(obj::YUMAParser)
    return jcall(obj, "loadData", void, ())
end

function load_data(obj::YUMAParser, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function still_accepts_data(obj::YUMAParser)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
using DiffableDistributions

@assert DNormal <: Distributions.AbstractNormal

m = ForwardDiff.GradientNumber(1)
s = ForwardDiff.GradientNumber(3.2)
@assert isa(DNormal(m, s), DNormal)
@assert isa(DNormal(m, 3.2), DNormal)
@assert isa(DNormal(1, s), DNormal)
@assert isa(DNormal(1, 2), DNormal)
@assert isa(Distributions.Normal(m, 2), DNormal)
@assert isa(Distributions.Normal(1., s), DNormal)
@assert isa(Distributions.Normal(m, s), DNormal)


text:
mutable struct LongstaffSchwartzExerciseStrategy{M <: MarketModelBasisSystem, E <: MarketModelExerciseValue, E2 <: MarketModelExerciseValue} <: ExerciseStrategy
  basisSystem::M
  basisCoefficients::Vector{Vector{Float64}}
  exercise::E
  control::E2
  numeraires::Vector{Int}
  currentIndex::Int
  principalInNumerairePortfolio::Float64
  newPrincipal::Float64
  exerciseTimes::Vector{Float64}
  relevantTimes::Vector{Float64}
  isBasisTime::BitArray{1}
  isRebateTime::BitArray{1}
  isControlTime::BitArray{1}
  isExerciseTime::BitArray{1}
  rebateDiscounters::Vector{MarketModelDiscounter}
  controlDiscounters::Vector{MarketModelDiscounter}
  basisValues::Vector{Vector{Float64}}
  exerciseIndex::Vector{Int}
end

function LongstaffSchwartzExerciseStrategy(basisSystem::MarketModelBasisSystem,
                                          basisCoefficients::Vector{Vector{Float64}},
                                          evolution::EvolutionDescription,
                                          numeraires::Vector{Int},
                                          exercise::MarketModelExerciseValue,
                                          control::MarketModelExerciseValue)

  basisSystem = clone(basisSystem)
  exercise = clone(exercise)
  control = clone(control)

  check_compatibility(evolution, numeraires)
  relevantTimes = evolution.evolutionTimes
  isBasisTime = is_in_subset(relevantTimes, basisSystem.evolution.evolutionTimes)
  isRebateTime = is_in_subset(relevantTimes, exercise.evolution.evolutionTimes)
  isControlTime = is_in_subset(relevantTimes, control.evolution.evolutionTimes)

  exerciseIndex = Vector{Int}(undef, length(relevantTimes))
  isExerciseTime = falses(length(exerciseIndex))
  exerciseTimes = Vector{Float64}()
  v = exercise.isExerciseTime
  exercises = idx = 1
  @inbounds @simd for i in eachindex(relevantTimes)
    exerciseIndex[i] = exercises
    if isRebateTime[i]
      isExerciseTime[i] = v[idx]
      idx += 1
      if isExerciseTime[i]
        push!(exerciseTimes, relevantTimes[i])
        exercises += 1
      end
    end
  end

  rateTimes = evolution.rateTimes
  rebateTimes = possible_cash_flow_times(exercise)
  @inbounds rebateDiscounters = MarketModelDiscounter[MarketModelDiscounter(rebateTimes[i], rateTimes) for i in eachindex(rebateTimes)]

  controlTimes = possible_cash_flow_times(control)
  @inbounds controlDiscounters = MarketModelDiscounter[MarketModelDiscounter(controlTimes[i], rateTimes) for i in eachindex(controlTimes)]

  basisSizes = number_of_functions(basisSystem)
  @inbounds basisValues = Vector{Float64}[zeros(basisSizes[i]) for i in eachindex(basisSizes)]

  return LongstaffSchwartzExerciseStrategy(basisSystem, basisCoefficients, exercise, control, numeraires, -1, 0.0, 0.0, exerciseTimes, relevantTimes,
                                          isBasisTime, isRebateTime, isControlTime, isExerciseTime, rebateDiscounters, controlDiscounters,
                                          basisValues, exerciseIndex)
end

relevant_times(ls::LongstaffSchwartzExerciseStrategy) = ls.relevantTimes

function reset!(ls::LongstaffSchwartzExerciseStrategy)
  reset!(ls.exercise)
  reset!(ls.control)
  reset!(ls.basisSystem)
  ls.currentIndex = 1
  ls.principalInNumerairePortfolio = ls.newPrincipal = 1.0

  return ls
end

function next_step!(ls::LongstaffSchwartzExerciseStrategy, currentState::CurveState)
  ls.principalInNumerairePortfolio = ls.newPrincipal

  if ls.isRebateTime[ls.currentIndex]
    next_step!(ls.exercise, currentState)
  end

  if ls.isControlTime[ls.currentIndex]
    next_step!(ls.control, currentState)
  end

  if ls.isBasisTime[ls.currentIndex]
    next_step!(ls.basisSystem, currentState)
  end

  if ls.currentIndex < length(ls.numeraires)
    numeraire = ls.numeraires[ls.currentIndex]
    nextNumeraire = ls.numeraires[ls.currentIndex + 1]
    ls.newPrincipal *= discount_ratio(currentState, numeraire, nextNumeraire)
  end

  ls.currentIndex += 1

  return ls
end

function get_exercise!(ls::LongstaffSchwartzExerciseStrategy, currentState::CurveState)
  exerciseIndex = ls.exerciseIndex[ls.currentIndex-1]
  # println("IDX: ", exerciseIndex)

  exerciseCF = get_value(ls.exercise, currentState)
  exerciseValue = exerciseCF.amount * numeraire_bonds(ls.rebateDiscounters[exerciseCF.timeIndex], currentState, ls.numeraires[ls.currentIndex - 1]) /
                  ls.principalInNumerairePortfolio

  controlCF = get_value(ls.control, currentState)
  controlValue = controlCF.amount * numeraire_bonds(ls.controlDiscounters[controlCF.timeIndex], currentState, ls.numeraires[ls.currentIndex - 1]) /
                  ls.principalInNumerairePortfolio

  set_values!(ls.basisSystem, currentState, ls.basisValues[exerciseIndex])

  alphas = ls.basisCoefficients[exerciseIndex]
  continuationValue = dot(alphas, ls.basisValues[exerciseIndex])

  # println("alphas: ", alphas)
  # println("basis: ", ls.basisValues[exerciseIndex])
  # println("exerciseVal: ", exerciseValue)
  # println("continuationValue: ", continuationValue)

  return exerciseValue >= continuationValue
end

clone(ls::LongstaffSchwartzExerciseStrategy) = LongstaffSchwartzExerciseStrategy(clone(ls.basisSystem), deepcopy(ls.basisCoefficients), clone(ls.exercise),
                                                clone(ls.control), copy(ls.numeraires), ls.currentIndex, ls.principalInNumerairePortfolio, ls.newPrincipal,
                                                copy(ls.exerciseTimes), copy(ls.relevantTimes), copy(ls.isBasisTime), copy(ls.isRebateTime), copy(ls.isControlTime),
                                                copy(ls.isExerciseTime), deepcopy(ls.rebateDiscounters), deepcopy(ls.controlDiscounters), deepcopy(ls.basisValues),
                                                copy(ls.exerciseIndex))


text:
function __init__()
    RLDatasets.d4rl_init()
    RLDatasets.d4rl_pybullet_init()
    RLDatasets.atari_init()
    RLDatasets.rl_unplugged_atari_init()
    RLDatasets.bsuite_init()
    RLDatasets.dm_init()
    RLDatasets.deep_ope_d4rl_init()
end


text:
# this package
using Projeto2Solvers

using Logging
using NLPModels
using Test
using CUTEst

include("unc.jl")


text:
# TODO: more tests!
using Test
using ChainRulesCore
using LinearAlgebra: Diagonal
using ChainRulesCore: extern, accumulate, accumulate!, store!, @scalar_rule,
    Wirtinger, wirtinger_primal, wirtinger_conjugate,
    Zero, One, Casted, cast,
    DNE, Thunk, Casted, DNERule, WirtingerRule
using Base.Broadcast: broadcastable

@testset "ChainRulesCore" begin
    include("differentials.jl")
    include("rules.jl")
    include("rule_types.jl")
end


text:
struct TestSummary{L, A}
    pvalues::Vector{Float64}
    keys::Vector{String}
    df::Vector{String}
    labels::Vector{String}
    statistic::Vector{String}
    value::Vector{String}
    locals::L
    accessor::A
end
struct Experiment{name, A, T}
    data_accessor::A
    target::T
    major_finding::String
    comments::String
    data::DataFrame
    tests::Dict{String,TestSummary}
    nbirds::Int
    foodtypes::Vector{Food}
end
function Experiment{name}(a::A, t::T, args...) where {name, A, T}
    Experiment{name, A, T}(a, t, args...)
end

name(::Experiment{N}) where N = N
nbirds(exp::Experiment) = exp.nbirds
foodtypes(exp::Experiment) = exp.foodtypes
target(exp::Experiment) = exp.target

default_keycols(x) = setdiff(names(x), ["μ", "sem", "firstinspection", "n"])
extract(data, da) = convert(Vector{Float64}, vcat([data[idx, c] for (c, idx) in da]...))
function results(da, data, tests, ::Any)
    res = extract(data, da)
    append!(res, significancecode.(tests.pvalues))
end
function results(da::Base.RefValue, data, tests, e::Experiment{K}) where K
    ks = default_keycols(data)
    targetdf = K ∈ (:Cheke11_specsat, :Cheke11_planning) ? summarize(e, e.data) : e.data
    if dropmissing(data[!,ks]) != dropmissing(targetdf[!,ks])
        @warn K data[!,ks] targetdf[!,ks]
    end
    da[] = [(c, findall(x -> x !== missing, targetdf[!, c]))
            for c in intersect(["μ", "sem", "firstinspection"], names(targetdf))]
    target = extract(targetdf, da[])
    for test in tests.tests
        targetdf = e.tests[test.id].df
        append!(target, significancecode.(e.tests[test.id].pvalues))
        a = e.tests[test.id].accessor
        isassigned(a) || continue
        computeddf = df(a[], test.content, nothing)
        targetdf == computeddf || @warn K test.id targetdf computeddf
    end
    e.target[] = target
    results(da[], data, tests, e)
end
function results(e, data)
    tests = statistical_tests(e, data)
    ag = summarize(e, data)
    results(e.data_accessor, ag, tests, e)
end
function results(e::Experiment{:Clayton0103}, data)
    res = vcat([results(EXPERIMENTS[key], res)
               for (key, res) in zip(CLAYTON0103_EXPERIMENTS, data)]...)
    if isa(e.target, Base.RefValue)
        e.target[] = vcat([EXPERIMENTS[key].target[] for key in CLAYTON0103_EXPERIMENTS]...)
        e.data_accessor[] = nothing
    end
    res
end
function asdataframe(e::Experiment{K}, x, k = Ref(0)) where K
    data = K ∈ (:Cheke11_specsat, :Cheke11_planning) ? summarize(e, e.data) : deepcopy(e.data)
    k0 = k[]
    for (col, idx) in e.data_accessor
        data[idx, col] .= x[k[]+1:k[]+length(idx)]
        k[] += length(idx)
    end
    tests = x[k[]+1:k0 + length(e.target)]
    k[] += length(tests)
    data, tests
end
function asdataframe(::Experiment{:Clayton0103}, x, k = Ref(0))
    [asdataframe(e, x, k) for e in CLAYTON0103_EXPERIMENTS]
end


text:
using Documenter
using ExternalDocstrings

makedocs(
    sitename = "ExternalDocstrings",
    format = Documenter.HTML(),
    modules = [ExternalDocstrings],
    # Ref:
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.makedocs
)

deploydocs(
    repo = "github.com/tkf/ExternalDocstrings.jl",
    devbranch = "main",
    push_preview = true,
    # Ref:
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
)


text:
"""
Kx Systems provides C API that can be used both on the server side (calling
functions defined in the q binary) or on the server side by linking with a
`c.o` object file.  With a few exceptions, the same functions are available
on the server and the client side and they are declared in a C header `k.h`.

This module is a thin wrapper for the kdb+ C API.
"""
module _k  # k.h wrappers
export k, b9, d9, okx, kclose
export ymd, dj
export r0, r1
export ktj, ka, kb, ku, kg, kh, ki, kj, ke, kf, sn, ss, ks, kc, kd, kz, kt
export ja, js, jk, jv
export ktn, knk, kp, xT, xD
export xa, xt, t, xr, r, xg, xh, xi, xj, xe, xf, xs, xn, xp, n, xk, xx, xy
export kG, kH, kI, kJ, kE, kF, kC, kS, kK, kX
export B_, C_, S_, G_, H_, I_, J_, E_, F_, V_, U_, K_
export KB, UU, KG, KH, KI, KJ, KE, KF, KC, KS, KP, KM, KD, KZ, KN, KU, KV, KT,
       XT, XD, KK, EE

include("startup.jl")

#########################################################################
# k.h
const B_ = Bool  # not in k.h
const C_ = Cchar
const S_ = Cstring
const G_ = Cuchar
const H_ = Cshort
const I_ = Cint
const J_ = Clonglong
const E_ = Cfloat
const F_ = Cdouble
const V_ = Void
const U_ = UInt128
struct k0
    m::C_
    a::C_
    t::C_  # type code
    u::C_
    r::I_  # reference count
end
const K_ = Ptr{k0}

TYPE_LETTERS = "kb  ghijefcspmdznuvt"
for (t, x) in enumerate(TYPE_LETTERS)
    isspace(x) || @eval const $(Symbol("K", uppercase(x))) = $(I_(t-1))
end
# guid
const UU = I_(2)
# table,dict
const XT = I_(98)   #   x->k is XD
const XD = I_(99)   #   kK(x)[0] is keys. kK(x)[1] is values.

const EE = I_(128)  #   error

Base.show(io::IO, ::Type{K_}) = write(io, "K_")

# reference management
"Increment the object's reference count"
r0(x::K_) = ccall((@k_sym :r0), K_, (K_,), x)
"Decrement the object's reference count"
r1(x::K_) = ccall((@k_sym :r1), K_, (K_,), x)

# head accessors
const xt = t(x::K_) = unsafe_load(x).t
const xr = r(x::K_) = unsafe_load(x).r

# scalar accessors
xg(x::K_) = unsafe_load(Ptr{G_}(x+8))
xh(x::K_) = unsafe_load(Ptr{H_}(x+8))
xi(x::K_) = unsafe_load(Ptr{I_}(x+8))
xj(x::K_) = unsafe_load(Ptr{J_}(x+8))
xe(x::K_) = unsafe_load(Ptr{E_}(x+8))
xf(x::K_) = unsafe_load(Ptr{F_}(x+8))
xs(x::K_) = (s = unsafe_load(Ptr{S_}(x+8)); s == C_NULL?"null":unsafe_string(s))

# vector accessors
const xn = n(x::K_) = unsafe_load(Ptr{J_}(x+8))
xp(x::K_) = unsafe_string(Ptr{UInt8}(x+16), xn(x))

kG(x::K_) = unsafe_wrap(Array, Ptr{G_}(x+16), (x|>n,))
kH(x::K_) = unsafe_wrap(Array, Ptr{H_}(x+16), (x|>n,))
kI(x::K_) = unsafe_wrap(Array, Ptr{I_}(x+16), (x|>n,))
kJ(x::K_) = unsafe_wrap(Array, Ptr{J_}(x+16), (x|>n,))
kE(x::K_) = unsafe_wrap(Array, Ptr{E_}(x+16), (x|>n,))
kF(x::K_) = unsafe_wrap(Array, Ptr{F_}(x+16), (x|>n,))
kC(x::K_) = unsafe_wrap(Array, Ptr{C_}(x+16), (x|>n,))
kS(x::K_) = unsafe_wrap(Array, Ptr{S_}(x+16), (x|>n,))
kK(x::K_) = unsafe_wrap(Array, Ptr{K_}(x+16), (x|>n,))
# Not in k.h, but useful
kX(::Type{C}, x::K_) where {C} = unsafe_wrap(Array, Ptr{C}(x+16), (x|>n,))

# table and dict accessors
xk(x::K_) = unsafe_load(Ptr{K_}(x+8))
xx(x::K_) = unsafe_load(Ptr{K_}(x+16), 1)
xy(x::K_) = unsafe_load(Ptr{K_}(x+16), 2)

# scalar constructors
"Create an atom of type and value"
ktj(t::Integer, x::Integer) = ccall((@k_sym :ktj), K_, (I_, J_), t, x)
"Create an atom of type"
ka(x::Integer) = ccall((@k_sym :ka), K_, (I_,), x)
"Create a boolean"
kb(x::Integer) = ccall((@k_sym :kb), K_, (I_,), x)
"Create a guid"
ku(x::U_) = (p = ka(-UU); unsafe_store!(Ptr{U_}(p+16), x); p)
ku(x::Integer) = ku(U_(x))
"Create a byte"
kg(x::Integer) = ccall((@k_sym :kg), K_, (I_,), x)
"Create a short"
kh(x::Integer) = ccall((@k_sym :kh), K_, (I_,), x)
"Create an int"
ki(x::Integer) = ccall((@k_sym :ki), K_, (I_,), x)
"Create a long"
kj(x::Integer) = ccall((@k_sym :kj), K_, (J_,), x)
"Create a real"
ke(x::Real) = ccall((@k_sym :ke), K_, (F_,), x)
"Create a float"
kf(x::Real) = ccall((@k_sym :kf), K_, (F_,), x)
"Create a char"
kc(x::Integer) = ccall((@k_sym :kc), K_, (I_,), x)
kc(x::Char) = kc(UInt8(x))
const _AnyString = Union{String, Symbol, Cstring}
"Intern n chars from a string"
sn(x::_AnyString, n::Integer) = ccall((@k_sym :sn), S_, (S_,I_), x, n)
"Intern a string"
ss(x::_AnyString) = ccall((@k_sym :ss), S_, (S_,), x)
"Create a symbol"
ks(x::_AnyString) = ccall((@k_sym :ks), K_, (S_,), x)
"Create a date"
kd(x::Integer) = ccall((@k_sym :kd), K_, (I_,), x)
"Create a datetime (deprecated)"
kz(x::AbstractFloat) = ccall((@k_sym :kz), K_, (F_,), x)
"Create a time"
kt(x::Integer) = ccall((@k_sym :kt), K_, (I_,), x)

# vector constructors
"Create a char array from string"
kp(x::AbstractString) = ccall((@k_sym :kp), K_, (S_,), x)
"Create a simple list of type and length"
ktn(t::Integer, n::Integer) = ccall((@k_sym :ktn), K_, (I_, J_), t, n)
#knk(n) = begin @assert n == 0; ccall((@k_sym :knk), K_, (I_,), 0) end
"Create a mixed list of length"
function knk(n::Integer, x::K_...)
    r = ktn(0, n)
    for i in 1:n
        unsafe_store!(Ptr{K_}(r+16), x[i], i)
    end
    return r
end
# table, dictionary
"Create a table from a dict"
xT(x::K_) = ccall((@k_sym :xT), K_, (K_, ), x)
"Create a dict"
xD(x::K_, y::K_) = ccall((@k_sym :xD), K_, (K_, K_), x, y)

# ja(K*,V*),js(K*,S),jk(K*,K),jv(K*k,K)
"Join an atom to a list"
ja(rx::Ref{K_}, y::Ref) = ccall((@k_sym :ja), K_, (Ref{K_}, Ptr{V_}), rx, y)
"Join a symbol to a list"
js(rx::Ref{K_}, y::S_) = ccall((@k_sym :js), K_, (Ref{K_}, S_), rx, y)
"Join a K_ object to a list"
jk(rx::Ref{K_}, y::K_) = ccall((@k_sym :jk), K_, (Ref{K_}, K_), rx, y)
"Join another K_ list to a list"
jv(rx::Ref{K_}, y::K_) = ccall((@k_sym :jv), K_, (Ref{K_}, K_), rx, y)

# K b9(I,K) and K d9(K)
b9(pe::Integer, x::K_) = ccall((@k_sym :b9), K_, (I_, K_), pe, x)
d9(x::K_) = ccall((@k_sym :d9), K_, (K_, ), x)
okx(x::K_) = ccall((@k_sym :okx), I_, (K_, ), x)
kclose(h::Integer) = ccall((@k_sym :kclose), V_, (I_, ), h)

# Dates
"Encode a year/month/day as q date"
ymd(y::Integer, m::Integer, d::Integer) = ccall((@k_sym :ymd),
    I_, (I_, I_, I_), y, m, d)
"Convert q date to yyyymmdd integer"
dj(j::Integer) = ccall((@k_sym :dj), I_, (I_, ), j)

if GOT_Q
    export dot_, #= avoid conflict with Base.dot. =# ee, dl, khp
    export krr, orr
    dot_(x::K_, y::K_) = ccall((@k_sym :dot), K_, (K_, K_), x, y)
    ee(x::K_) = ccall((@k_sym :ee), K_, (K_, ), x)
    # dl(V*f,I)
    dl(f::Ptr{V_}, n::Integer) = ccall((@k_sym :dl), K_, (Ptr{V_}, I_), f, n)
    # Simulate khp with hopen
    function khp(h::String, p::Integer)
        x = k(0, "hopen", ks(string(":", h, ":", p)))
        try
            return xi(x)
        finally
            r0(x)
        end
    end
    krr(x::AbstractString) = ccall((@k_sym :krr), K_, (S_, ), x)
    orr(x::AbstractString) = ccall((@k_sym :orr), K_, (S_, ), x)
else
    # communications (not included in q server)
    export khpun, khpu, khp, krr, orr
    # I khpun(const S,I,const S,I),khpu(const S,I,const S),khp(const S,I)
    khpun(h::String, p::Integer, u::String, n::Integer) = ccall((@k_sym :khpu),
        I_, (S_, I_, S_, I_), h, p, u, n)
    khpu(h::String, p::Integer, u::String) = ccall((@k_sym :khpu),
        I_, (S_, I_, S_), h, p, u)
    khp(h::String, p::Integer) = ccall((@k_sym :khp), I_, (S_, I_), h, p)
    function krr(x::AbstractString)
        e = ka(-128)
        unsafe_store!(Ptr{S_}(e+8), ss(x))
        e
    end
    orr(x::AbstractString) = krr(String(x, ": ", Libc.strerror()))
end

const K_NULL = K_(C_NULL)
# K k(I,const S,...)
# TODO: Use Julia metaprogramming to avoid repetition
k(h::Integer, m::String) = ccall((@k_sym :k), K_, (I_, S_, K_), h, m, K_NULL)
k(h::Integer, m::String, x1::K_) = ccall((@k_sym :k),
    K_, (I_, S_, K_, K_), h, m, x1, K_NULL)
k(h::Integer, m::String, x1::K_, x2::K_) = ccall((@k_sym :k),
    K_, (I_, S_, K_, K_, K_), h, m, x1, x2, K_NULL)
k(h::Integer, m::String, x1::K_, x2::K_, x3::K_) = ccall((@k_sym :k),
    K_, (I_, S_, K_, K_, K_, K_), h, m, x1, x2, x3, K_NULL)
k(h::Integer, m::String, x1::K_, x2::K_, x3::K_, x4::K_) = ccall((@k_sym :k),
    K_, (I_, S_, K_, K_, K_, K_, K_), h, m, x1, x2, x3, x4, K_NULL)
k(h::Integer, m::String,
    x1::K_, x2::K_, x3::K_, x4::K_, x5::K_) = ccall((@k_sym :k),
        K_, (I_, S_, K_, K_, K_, K_, K_, K_), h, m, x1, x2, x3, x4, x5, K_NULL)
k(h::Integer, m::String,
    x1::K_, x2::K_, x3::K_, x4::K_, x5::K_, x6::K_) = ccall((@k_sym :k),
        K_, (I_, S_, K_, K_, K_, K_, K_, K_, K_),
            h, m, x1, x2, x3, x4, x5, x6, K_NULL)
k(h::Integer, m::String,
    x1::K_, x2::K_, x3::K_, x4::K_, x5::K_, x6::K_, x7::K_) = ccall((@k_sym :k),
        K_, (I_, S_, K_, K_, K_, K_, K_, K_, K_, K_),
            h, m, x1, x2, x3, x4, x5, x6, x7, K_NULL)
k(h::Integer, m::String,
    x1::K_, x2::K_, x3::K_, x4::K_, x5::K_, x6::K_,
        x7::K_, x8::K_) = ccall((@k_sym :k),
            K_, (I_, S_, K_, K_, K_, K_, K_, K_, K_, K_, K_),
                h, m, x1, x2, x3, x4, x5, x6, x7, x8, K_NULL)
# Make sure we don't redefine basic pointer conversions.
using Core.Intrinsics.bitcast
Base.convert(::Type{K_}, p::K_) = p
Base.convert(::Type{K_}, p::Ptr) = bitcast(K_, p)
Base.convert(::Type{K_}, p::UInt) = bitcast(K_, p)
end # module k


text:
true # true
false # false

1 # 1
typeof(1) # Int64

1.0 # 1.0
typeof(1.0) # Float64

1 == 1.0 # true
1 === 1.0 # false

isa(1, Int64) # true

isa(1, Int) # true

typeof("a") # ASCIIString
typeof('a') # Char
typeof("λ") # UTF8String

# Vector (i.e. 1-d array)

a1 = ["a", "b", "c"]
typeof(a1) # Array{ASCIIString,1}

push!(a1, "d")
a1 # ["a", "b", "c", "d"]

a1[0] # BoundsError (Julia arrays are 1-indexed)
a1[1] # "a"
a1[end] # "d"
a1[endof(a1)] # "d"

a1[1:2] # First two elements
a1[end-1:end] # Last two elements

a2 = ["e", "f", "g", "h"]

vcat(a1, a2) # ["a", "b", ..., "g", "h"]

# Tuple

t1 = ("a", "b", "c")
typeof(t1) # Tuple{ASCIIString,ASCIIString,ASCIIString}
t2 = ("a", 1.0, :hello)
typeof(t2) # Tuple{ASCIIString,Float64,Symbol}

# Dictionary

d1 = Dict("a"=>1, "b"=>2)
typeof(d1) # Dict{Symbol,Int64}
d1["c"] = 3

# Set
s1 = Set([1, 2, 2, 1, 3, 2])
s1 # e.g. Set([2,3,1])


text:
module TestNChisq

using Base.Test
using Econometrics

expOut = 0.1441868
actOut = Econometrics.dnchisq(1.2, 3.2, 1.4)
@test_approx_eq_eps expOut actOut 1e-6

expOut = 0.1244554
actOut = Econometrics.pnchisq(1.2, 3.2, 1.4)
@test_approx_eq_eps expOut actOut 1e-6

expOut = 3.039873
actOut = Econometrics.qnchisq(0.4, 3.2, 1.4)
@test_approx_eq_eps expOut actOut 1e-6

p = 0.2
q = Econometrics.qnchisq(p, 3.2, 1.4)
p2 = Econometrics.pnchisq(q, 3.2, 1.4)
@test_approx_eq p p2

simVals = Econometrics.rnchisq(1000, 3.2, 1.4)
uVals = Econometrics.pnchisq(simVals, 3.2, 1.4)
simVals2 = Econometrics.qnchisq(uVals, 3.2, 1.4)
@test all(abs(simVals - simVals2) .< 1e-10)

end


text:
import Base:
    getindex, push!, close, flush

type DiskWriter
    offset::IOStream
    data::IOStream
end

function DiskWriter(filename)
    DiskWriter(open(filename * ".offset", "a"), open(filename, "a"))
end

function push!(db::DiskWriter, d::String)
    write(db.data, d)
    write(db.offset, position(db.data))
end

function flush(db::DiskWriter)
    flush(db.data)
    flush(db.offset)
end

function close(db::DiskWriter)
    close(db.data)
    close(db.offset)
end

type JsonDiskReader
    offset::Vector{Int}
    data::Vector{UInt8}
end

function DiskReader(filename)
    JsonDiskReader(Mmap.mmap(filename * ".offset", Vector{Int}), Mmap.mmap(filename, Vector{UInt8}))
end

function getindex(r::JsonDiskReader, i)::String
    off = i == 1 ? 1 : r.offset[i - 1]
    s = r.data[off:r.offset[i]]
    String(s)
end

function getindex(r::JsonDiskReader, indexes...)
    [getindex(r, i) for i in indexes]
end

max_stars_repo_name: INGEOTEC/ranix

max_issues_repo_name: INGEOTEC/ranix

max_forks_repo_name: INGEOTEC/ranix


text:
"""
    rainfarm(r,slope,nf,weight=1.;fglob=false, fsmooth=false, verbose=false)

Perform general RainFARM downscaling.

#Arguments
* `r`      : large-scale array to downscale
* `slope`  : spatial spectral slope
* `nf`     : refinement factor for spatial downscaling
* `weight` : weights for orographic downscaling
* `fglob`  : conserve global average over domain
* `fsmooth`: use smoothing instead of gp conservation
* `verbose`: provide some progress report
"""
function rainfarm(r, slope, nf, weight=1.; fglob=false, fsmooth=false, verbose=false)

(nax,nay,ntime)=size(r[:,:,:]);

if( nax==nay)
   nas=nax
else
   error("The spatial dimensions of r(nax,nay,nat) must be square (nax==nay)")
end

ns=nas*nf
nt=1; nat=1;
# This is the space_only version, downscaling nt is 1
#(ns, ns1)=size(lon); ns=max(ns,ns1);
#ns=size(lon); ns=ns[1];
#Recover ns (the fine scale res) from lon; we assume that RF works on squares, nlon=nlat
f=initmetagauss(slope,1,ns,1);
rd=zeros(ns,ns,ntime);

  for k=1:ntime
    r1=r[:,:,k];
    if verbose
@compat       @printf("Frame %d\r",k)
    end
    if Statistics.mean(r1)==0
	rd[:,:,k]=zeros(ns,ns);
    else
   	fm=downscale_spaceonly(r1,f,weight,fglob=fglob,fsmooth=fsmooth);
	rd[:,:,k]=fm;
    end
  end
  return rd
end


text:
"""
    testreport(nam::String,ext="")

Run the testreport script for one model config `nam` (or "all"),
with additional options (optional) speficied in `ext`

```
using MITgcmTools
testreport(MITgcm_config(configuration="front_relax"),"-norun")
#testreport(MITgcm_config(configuration="all"),"-norun")
```
"""
function testreport(config::MITgcm_config,ext="")
    nm=config.configuration
    try
        pth=pwd()
    catch e
        cd()
    end
    pth=pwd()
    cd(tempdir())
    println(pwd())
    if nm!=="all"
        lst=[nm]
    else
        exps=verification_experiments()
        lst=[exps[i].configuration for i in 1:length(exps)]
    end
    for nm in lst
        c=`$(MITgcm_path[1])/verification/testreport -t $(MITgcm_path[1])/verification/$(nm) $ext`
        isempty(ext) ? c=`$(MITgcm_path[1])/verification/testreport -t $(MITgcm_path[1])/verification/$(nm)` : nothing
        @suppress run(c)
    end
    cd(pth)
    return true
end

import ClimateModels: compile, build, setup, clean

"""
    clean(config::MITgcm_config)

Cancel any remaining task (config.channel) and clean up the run directory (via rm).

(part of the climate model interface as specialized for `MITgcm`)
"""
function clean(config::MITgcm_config)
    #cancel any remaining task
    while !isempty(config.channel)
        take!(config.channel)
    end
    #clean up run directory
    pp=joinpath(config.folder,string(config.ID),"run")
    isdir(pp) ? rm(pp,recursive=true) : nothing
    #
    return "no task left in pipeline"
end

"""
    build(config::MITgcm_config)

Build the model using `genmake2`, `make depend`, and `make`. The first two link all
code files, headers, etc  in the `build/` folder before compiling the model

(part of the climate model interface as specialized for `MITgcm`)
"""
function build(config::MITgcm_config)
    nam=config.configuration
    try
        pth=pwd()
    catch e
        cd()
    end
    pth=pwd()
    cd("$(MITgcm_path[1])/verification/$(nam)/build")
    try
        @suppress run(`../../../tools/genmake2 -mods=../code`) #$ext
        @suppress run(`make clean`)
        @suppress run(`make depend`)
        @suppress run(`make -j 4`)
    catch e
        println("model compilation may have failed")
    end
    cd(pth)
    return true
end

"""
    build(config::MITgcm_config,options::String)

Build the model using `genmake2`, `make depend`, and `make` unless otherwise
specified via `options`. The `genmake2` and `make depend` commands link all
code files, headers, etc  in the `build/` folder before `make` compiles the model.

(part of the climate model interface as specialized for `MITgcm`)
"""
function build(config::MITgcm_config,options::String)
    nam=config.configuration
    if options=="--allow-skip"
        tst=!isfile(joinpath(MITgcm_path[1],"verification",nam,"build","mitgcmuv"))
        tst ? build(config) : nothing
    else
        build(config)
    end
    return true
end

"""
    compile(config::MITgcm_config)

Compile the model using `make` in `build/` that has already been setup.

(part of the climate model interface as specialized for `MITgcm`)
"""
function compile(config::MITgcm_config)
    nam=config.configuration
    try
        pth=pwd()
    catch e
        cd()
    end
    pth=pwd()
    cd("$(MITgcm_path[1])/verification/$(nam)/build")
    try
        @suppress run(`make -j 4`)
    catch e
        println("model compilation may have failed")
    end
    cd(pth)
    return true
end

"""
    setup(config::MITgcm_config)

Create a `run/` folder and link everything there as needed to be ready to run model as
normally done for most-standard MITgcm configurations (incl. `prepare_run` and `mitgcmuv`).
Call `ClimateModels.git_log_init(config)` to setup git tracker and
`put!(config.channel,MITgcm_launch)` to be executed via `launch(config)` later.

(part of the climate model interface as specialized for `MITgcm`)
"""
function setup(config::MITgcm_config)
    !isdir(joinpath(config.folder)) ? mkdir(joinpath(config.folder)) : nothing
    !isdir(joinpath(config.folder,string(config.ID))) ? mkdir(joinpath(config.folder,string(config.ID))) : nothing

    pth_run=joinpath(config.folder,string(config.ID),"run")
    !isdir(pth_run) ? mkdir(pth_run) : nothing

    pth_log=joinpath(config.folder,string(config.ID),"log","tracked_parameters")
    pth_mv=joinpath(config.folder,string(config.ID),"original_parameters")

    if !isfile(joinpath(pth_run,"data"))&&isfile(joinpath(pth_log,"data"))
        p=pth_log
        f=readdir(p)
        [symlink(joinpath(p,f[i]),joinpath(pth_run,f[i])) for i in 1:length(f)]
    end

    p="$(MITgcm_path[1])/verification/$(config.configuration)/input"
    tmpA=readdir(p)
    f=tmpA[findall([!isfile(joinpath(pth_run,tmpA[i])) for i in 1:length(tmpA)])]
    [symlink(joinpath(p,f[i]),joinpath(pth_run,f[i])) for i in 1:length(f)]

    #replace relative paths with absolutes then exe prepare_run
    if isfile(joinpath(pth_run,"prepare_run"))
		try
			pth=pwd()
		catch e
			cd()
		end
        pth=pwd()
        cd(pth_run)
        #
        fil="prepare_run"
        meta = read(fil,String)
        meta = split(meta,"\n")
        ii=findall(occursin.("../../",meta))
        for i in ii
            meta[i]=replace(meta[i],"../../" => "$(MITgcm_path[1])/verification/")
        end
        ii=findall(occursin.("../",meta))
        for i in ii
            meta[i]=replace(meta[i],"../" => "$(MITgcm_path[1])/verification/$(config.configuration)")
        end
        #rm old file from run dir
        rm(fil)
        #write new file in run dir
        txt=["$(meta[i])\n" for i in 1:length(meta)]
        fid = open(fil, "w")
		[write(fid,txt[i]) for i in 1:length(txt)]
    	close(fid)
        #execute prepare_run
        chmod(fil,0o777)
        @suppress run(`./$(fil)`)
        #
        cd(pth)
    end

    if !islink(joinpath(pth_run,"mitgcmuv"))
        f="$(MITgcm_path[1])/verification/$(config.configuration)/build/mitgcmuv"
        symlink(f,joinpath(pth_run,"mitgcmuv"))
    end

    logdir=joinpath(config.folder,string(config.ID),"log")
    !isdir(logdir) ? ClimateModels.git_log_init(config) : nothing

    #Replace namelists with editeable versions in log/
    #
    #- read from run folder, rewrite to log/parameter_files
    #- mv all namelists to ../original_parameter_files
    #- link from log/parameter_files to here (run/)
    #(- add to git with message = original params)

    function list_namelist_files(pth_run)
            tmpA=readdir(pth_run)
            tmpA=tmpA[findall([length(tmpA[i])>3 for i in 1:length(tmpA)])]
            tmpA=tmpA[findall([tmpA[i][1:4]=="data"||tmpA[i]=="eedata"||
                    tmpA[i]=="prepare_run" for i in 1:length(tmpA)])]
    end
    nmlfiles=list_namelist_files(pth_run)

    if !isdir(pth_log)
        mkdir(pth_log)

        params=OrderedDict()
        for fil in nmlfiles
            nml=read(joinpath(pth_run,fil),MITgcm_namelist())
            write(joinpath(pth_log,fil),nml)
            #
            ni=length(nml.groups); tmp1=OrderedDict()
            [push!(tmp1,(nml.groups[i] => nml.params[i])) for i in 1:ni]
            tmp2=""
            fil=="data" ? tmp2="main" : nothing
            fil=="eedata" ? tmp2="eedata" : nothing
            occursin("data.",fil) ? tmp2=fil[6:end] : nothing
            if !isempty(tmp2)
                push!(params,(Symbol(tmp2) => tmp1))
            end
        end

        push!(config.inputs,params...)

        !isdir(pth_mv) ? mkdir(pth_mv) : nothing
        for fil in nmlfiles
            mv(joinpath(pth_run,fil),joinpath(pth_mv,fil))
            symlink(joinpath(pth_log,fil),joinpath(pth_run,fil))
        end

        ClimateModels.git_log_prm(config)
    end

    #add model run to scheduled tasks
    put!(config.channel,MITgcm_launch)

    return true
end

"""
    MITgcm_launch(config::MITgcm_config)

Go to `run/` folder and effectively call `mitgcmuv > output.txt`

(part of the climate model interface as specialized for `MITgcm`)
"""
function MITgcm_launch(config::MITgcm_config)
    try
        pth=pwd()
    catch e
        cd()
    end
    pth=pwd()
    cd(joinpath(config.folder,string(config.ID),"run"))
    tmp=["STOP NORMAL END"]
    try
        @suppress run(pipeline(`./mitgcmuv`,"output.txt"))
    catch e
        tmp[1]="model run may have failed"
    end
    cd(pth)
    return tmp[1]
end


text:
using Test

using DomainSets, DomainIntegrals, GridArrays, StaticArrays

using BasisFunctions
using BasisFunctions.Test: generic_test_measure

@testset "Measure" begin
    μ = BasisFunctions.GenericWeight(UnitInterval(),x->1.)
    generic_test_measure(μ)
    # @test_throws ErrorException !(isnormalized(μ))

    μ = LebesgueDomain(UnitInterval())
    generic_test_measure(μ)

    μ = LegendreWeight()
    generic_test_measure(μ)
    @test !isnormalized(μ)

    μ = FourierWeight()
    generic_test_measure(μ)
    @test isnormalized(μ)

    @test lebesguemeasure(UnitInterval()) isa FourierWeight
    @test lebesguemeasure(ChebyshevInterval()) isa LegendreWeight
    @test lebesguemeasure(0.3..0.4) isa LebesgueDomain

    μ = ChebyshevTWeight()
    generic_test_measure(μ)
    @test !isnormalized(μ)

    μ = ChebyshevUWeight()
    generic_test_measure(μ)
    @test !isnormalized(μ)

    μ = JacobiWeight(rand(),rand())
    generic_test_measure(μ)
    @test !isnormalized(μ)

    m = mapto(0..1, -1..1)
    μ = mappedmeasure(m,FourierWeight())
    generic_test_measure(μ)

    m = FourierWeight()
    μ = productmeasure(m,m)
    @test components(μ) == (component(μ,1),component(μ,2))
    io = IOBuffer()
    show(io, μ)
    @test length(take!(io))>0
    support(μ)
    x = SVector(rand(),rand())
    @test weightfun(μ,x)≈weightfunction(μ)(x)
end


text:
function _task_id(tname::String)
    _date = Dates.now()
    tstr = Dates.format(_date, "yyyy-mm-dd-HHMMSS")
    return isempty(tname) ? tstr : string(tstr, "-", tname)
end

max_stars_repo_name: josePereiro/GitWorkers

max_issues_repo_name: josePereiro/GitWorkers

max_forks_repo_name: josePereiro/GitWorkers


text:
# takes a topology, which for now is just an adjacency list; also list of neurons, but we
# should be able to just auto-detect all the neurons in the topology
function Network(neurons, topology; name = Base.gensym(:Network))

    all_neurons = Set(getproperty.(neurons, :sys))
    eqs = Equation[]
    params = Set()
    states = Set() # place holder
    defaultmap = Pair[]
    all_systems = AbstractSystem[]
    push!(all_systems, all_neurons...)
    post_neurons = Set()

    # what types of synapses do we have
    all_synapses = [synapse.second[2] for synapse in topology]
    all_synapse_types = [nameof(x.sys) for x in all_synapses]
    synapse_types = unique(all_synapse_types)

    # how many of each synapse type are there
    synapse_counts = Dict{Symbol,Int64}()

    for n in synapse_types
        c = count(x -> isequal(x, n), all_synapse_types)
        push!(synapse_counts, n => c)
    end

    # Extract reversal potentials
    reversals = unique([x.reversal for x in all_synapses])
    push!(params, reversals...)
    rev_meta = [getmetadata(x, ConductorEquilibriumCtx).val for x in reversals]

    for (i,j) in zip(reversals, rev_meta)
        push!(defaultmap, i => ustrip(Float64, mV, j))
    end

    voltage_fwds = Set()

    # create a unique synaptic current for each post-synaptic target
    for synapse in topology
        pre = synapse.first.sys # pre-synaptic neuron system
        post = synapse.second[1].sys # post-synaptic neuron system
        push!(post_neurons, post)
        Erev = synapse.second[2].reversal
        syntype = synapse.second[2].sys # synapse system
        synname = nameof(syntype)
        syn = @set syntype.name = Symbol(synname, synapse_counts[synname]) # each synapse is uniquely named
        synapse_counts[synname] -= 1
        push!(all_systems, syn)
        push!(voltage_fwds, syn.Vₘ ~ pre.Vₘ)

        if post.Isyn ∈ Set(vcat((get_variables(x.lhs) for x in eqs)...))
            idx = findfirst(x -> isequal([post.Isyn], get_variables(x.lhs)), eqs)
            eq = popat!(eqs, idx)
            eq = eq.lhs ~ eq.rhs + (syn.g * (post.Vₘ - Erev))
            push!(eqs, eq)
        else
            push!(eqs, post.Isyn ~ syn.g * (post.Vₘ - Erev))
        end
    end

    for nonpost in setdiff(all_neurons, post_neurons)
        push!(eqs, D(nonpost.Isyn) ~ 0)
    end

    append!(eqs, collect(voltage_fwds))
    network_system = ODESystem(eqs, t, states, params; defaults = defaultmap, name = name )
    return compose(network_system, all_systems)
end


text:
#
# generate_cwrapper.jl for tetgen copied from JuliaGeometry/tetgenbuilder
# In the moment this file is dormant.
#

# (c) Simon Danisch

# Permission is hereby granted, free  of charge, to any person obtaining
# a  copy  of this  software  and  associated documentation  files  (the
# "Software"), to  deal in  the Software without  restriction, including
# without limitation  the rights to  use, copy, modify,  merge, publish,
# distribute, sublicense,  and/or sell  copies of  the Software,  and to
# permit persons to whom the Software  is furnished to do so, subject to
# the following conditions:

# The  above  copyright  notice  and this  permission  notice  shall  be
# included in all copies or substantial portions of the Software.

# THE  SOFTWARE IS  PROVIDED  "AS  IS", WITHOUT  WARRANTY  OF ANY  KIND,
# EXPRESS OR  IMPLIED, INCLUDING  BUT NOT LIMITED  TO THE  WARRANTIES OF
# MERCHANTABILITY,    FITNESS    FOR    A   PARTICULAR    PURPOSE    AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING  FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

cd(@__DIR__)
open("cwrapper.cpp", "w") do io
    println(io, """
    #include "tetgen.h" // Defined tetgenio, tetrahedralize().
    """)

    for (REAL, ext) in ["double" => "f64"]
        println(io, """
        typedef struct {
          int firstnumber; // 0 or 1, default 0.
          int mesh_dim;    // must be 3.

          $REAL *pointlist;
          $REAL *pointattributelist;
          $REAL *pointmtrlist;
          int *pointmarkerlist;
          int numberofpoints;
          int numberofpointattributes;
          int numberofpointmtrs;

          int *tetrahedronlist;
          $REAL *tetrahedronattributelist;
          $REAL *tetrahedronvolumelist;
          int *neighborlist;
          int numberoftetrahedra;
          int numberofcorners;
          int numberoftetrahedronattributes;

          tetgenio::facet *facetlist;
          int *facetmarkerlist;
          int numberoffacets;

          $REAL *holelist;
          int numberofholes;

          $REAL *regionlist;
          int numberofregions;

          $REAL *facetconstraintlist;
          int numberoffacetconstraints;

          $REAL *segmentconstraintlist;
          int numberofsegmentconstraints;

          int *trifacelist;
          int *trifacemarkerlist;
          int numberoftrifaces;

          int *edgelist;
          int *edgemarkerlist;
          int numberofedges;
        } TetGenIO$ext;

        void copy_tetio(tetgenio* in, TetGenIO$ext* out){
            out->firstnumber = in->firstnumber;
            out->mesh_dim = in->mesh_dim;

            out->pointlist = in->pointlist;
            out->pointattributelist = in->pointattributelist;
            out->pointmtrlist = in->pointmtrlist;
            out->pointmarkerlist = in->pointmarkerlist;
            out->numberofpoints = in->numberofpoints;
            out->numberofpointattributes = in->numberofpointattributes;
            out->numberofpointmtrs = in->numberofpointmtrs;

            out->tetrahedronlist = in->tetrahedronlist;
            out->tetrahedronattributelist = in->tetrahedronattributelist;
            out->tetrahedronvolumelist = in->tetrahedronvolumelist;
            out->neighborlist = in->neighborlist;
            out->numberoftetrahedra = in->numberoftetrahedra;
            out->numberofcorners = in->numberofcorners;
            out->numberoftetrahedronattributes = in->numberoftetrahedronattributes;

            out->facetlist = in->facetlist;
            out->facetmarkerlist = in->facetmarkerlist;
            out->numberoffacets = in->numberoffacets;

            out->holelist = in->holelist;
            out->numberofholes = in->numberofholes;

            out->regionlist = in->regionlist;
            out->numberofregions = in->numberofregions;

            out->facetconstraintlist = in->facetconstraintlist;
            out->numberoffacetconstraints = in->numberoffacetconstraints;

            out->segmentconstraintlist = in->segmentconstraintlist;
            out->numberofsegmentconstraints = in->numberofsegmentconstraints;

            out->trifacelist = in->trifacelist;
            out->trifacemarkerlist = in->trifacemarkerlist;
            out->numberoftrifaces = in->numberoftrifaces;

            out->edgelist = in->edgelist;
            out->edgemarkerlist = in->edgemarkerlist;
            out->numberofedges = in->numberofedges;

        }
        void copy_tetio(TetGenIO$ext* in, tetgenio* out){
            out->firstnumber = in->firstnumber;
            out->mesh_dim = in->mesh_dim;

            out->pointlist = in->pointlist;
            out->pointattributelist = in->pointattributelist;
            out->pointmtrlist = in->pointmtrlist;
            out->pointmarkerlist = in->pointmarkerlist;
            out->numberofpoints = in->numberofpoints;
            out->numberofpointattributes = in->numberofpointattributes;
            out->numberofpointmtrs = in->numberofpointmtrs;

            out->tetrahedronlist = in->tetrahedronlist;
            out->tetrahedronattributelist = in->tetrahedronattributelist;
            out->tetrahedronvolumelist = in->tetrahedronvolumelist;
            out->neighborlist = in->neighborlist;
            out->numberoftetrahedra = in->numberoftetrahedra;
            out->numberofcorners = in->numberofcorners;
            out->numberoftetrahedronattributes = in->numberoftetrahedronattributes;

            out->facetlist = in->facetlist;
            out->facetmarkerlist = in->facetmarkerlist;
            out->numberoffacets = in->numberoffacets;

            out->holelist = in->holelist;
            out->numberofholes = in->numberofholes;

            out->regionlist = in->regionlist;
            out->numberofregions = in->numberofregions;

            out->facetconstraintlist = in->facetconstraintlist;
            out->numberoffacetconstraints = in->numberoffacetconstraints;

            out->segmentconstraintlist = in->segmentconstraintlist;
            out->numberofsegmentconstraints = in->numberofsegmentconstraints;

            out->trifacelist = in->trifacelist;
            out->trifacemarkerlist = in->trifacemarkerlist;
            out->numberoftrifaces = in->numberoftrifaces;

            out->edgelist = in->edgelist;
            out->edgemarkerlist = in->edgemarkerlist;
            out->numberofedges = in->numberofedges;

        }
        extern "C" {

          TetGenIO$ext tetrahedralize$ext(TetGenIO$ext jl_in, char* command){
            tetgenio in, out;
            copy_tetio(&jl_in, &in);
            tetrahedralize(command, &in, &out);
            TetGenIO$ext jl_out;
            copy_tetio(&out, &jl_out);
            in.initialize(); // don't free pointers from julia!
            out.initialize(); // don't free pointers for julia!
            return jl_out;
          }

        }

        """)
    end
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SSGraphBLAS_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SSGraphBLAS")
JLLWrappers.@generate_main_file("SSGraphBLAS", UUID("7ed9a814-9cab-54e9-8e9e-d9e95b4d61b1"))
end  # module SSGraphBLAS_jll


text:
# ---
# title: Loading Images
# description: Loading FITS images from files.
# author: "[William Thompson](https://github.com/sefffal)"
# cover: assets/loading-images.png
# ---

# We'll start by downloading a sample image. If you have an image stored locally,
# you would skip this step.
using AstroImages
fname = download(
    "http://www.astro.uvic.ca/~wthompson/astroimages/fits/656nmos.fits",
    "eagle-656nmos.fits"
);

# Load the image by filename.
# If unspecified, the image is loaded from the first image-HDU in the fits file.
img = AstroImage("eagle-656nmos.fits")

# --- save covers --- #src
mkpath("assets")  #src
save("assets/loading-images.png", imview(img)) #src


text:
module Uprating

"""
Semi-complete indexing routine using OBR quarterly data.
"""

include( "prices.jl" )

export uprate, UPRATE_MAPPINGS

end


text:
using ..Atmos
using ..Atmos: thermo_state, turbulence_tensors
using ..Mesh.Topologies
using ..Mesh.Grids
using ..MoistThermodynamics
using CLIMAParameters.Atmos.SubgridScale: inv_Pr_turb
using LinearAlgebra

Base.@kwdef mutable struct AtmosCollectedDiagnostics
    zvals::Union{Nothing, Array} = nothing
    repdvsr::Union{Nothing, Array} = nothing
end
const CollectedDiagnostics = AtmosCollectedDiagnostics()

include("diagnostic_vars.jl")

function atmos_default_init(dgngrp::DiagnosticsGroup, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    FT = eltype(Q)
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    if Array ∈ typeof(Q).parameters
        localvgeo = grid.vgeo
    else
        localvgeo = Array(grid.vgeo)
    end

    CollectedDiagnostics.zvals = zeros(FT, Nqk * nvertelem)
    CollectedDiagnostics.repdvsr = zeros(FT, Nqk * nvertelem)

    @visitQ nhorzelem nvertelem Nqk Nq begin
        z = localvgeo[ijk, grid.x3id, e]
        MH = localvgeo[ijk, grid.MHid, e]
        CollectedDiagnostics.zvals[Nqk * (ev - 1) + k] += MH * z
        CollectedDiagnostics.repdvsr[Nqk * (ev - 1) + k] += MH
    end

    # compute the full number of points on a slab
    MPI.Allreduce!(CollectedDiagnostics.repdvsr, +, mpicomm)

    CollectedDiagnostics.zvals ./= CollectedDiagnostics.repdvsr
end

# thermodynamic variables of interest
function vars_thermo(FT)
    @vars begin
        q_liq::FT
        q_ice::FT
        q_vap::FT
        T::FT
        θ_liq_ice::FT
        θ_dry::FT
        θ_v::FT
        e_int::FT
        h_m::FT
        h_t::FT
    end
end
num_thermo(FT) = varsize(vars_thermo(FT))
thermo_vars(array) = Vars{vars_thermo(eltype(array))}(array)

function compute_thermo!(bl, state, aux, ijk, e, thermoQ)
    e_tot = state.ρe / state.ρ
    ts = thermo_state(bl, state, aux)
    e_int = internal_energy(ts)
    Phpart = PhasePartition(ts)

    th = thermo_vars(thermoQ[ijk, e])
    th.q_liq = Phpart.liq
    th.q_ice = Phpart.ice
    th.q_vap = vapor_specific_humidity(ts)
    th.T = air_temperature(ts)
    th.θ_liq_ice = liquid_ice_pottemp(ts)
    th.θ_dry = dry_pottemp(ts)
    th.θ_v = virtual_pottemp(ts)
    th.e_int = e_int

    # Moist and total henthalpy
    R_m = gas_constant_air(ts)
    th.h_m = e_int + R_m * th.T
    th.h_t = e_tot + R_m * th.T

    return nothing
end

# horizontal averages
function vars_horzavg(FT)
    @vars begin
        ρ::FT
        ρu::FT
        ρv::FT
        ρw::FT
        e_tot::FT
        ρq_tot::FT
        q_liq::FT
        q_vap::FT
        θ_liq_ice::FT
        θ_dry::FT
        θ_v::FT
        e_int::FT
        h_m::FT
        h_t::FT
        qt_sgs::FT
        ht_sgs::FT
    end
end
num_horzavg(FT) = varsize(vars_horzavg(FT))
horzavg_vars(array) = Vars{vars_horzavg(eltype(array))}(array)

function compute_horzsums!(
    atmos::AtmosModel,
    state,
    diffusive_flux,
    aux,
    k,
    ijk,
    ev,
    e,
    Nqk,
    nvertelem,
    MH,
    localaux,
    thermoQ,
    horzsums,
    LWP,
    currtime,
)
    th = thermo_vars(thermoQ[ijk, e])
    hs = horzavg_vars(horzsums[Nqk * (ev - 1) + k])
    hs.ρ += MH * state.ρ
    hs.ρu += MH * state.ρu[1]
    hs.ρv += MH * state.ρu[2]
    hs.ρw += MH * state.ρu[3]
    hs.e_tot += MH * state.ρe
    hs.q_liq += MH * th.q_liq
    hs.q_vap += MH * th.q_vap
    hs.θ_liq_ice += MH * th.θ_liq_ice
    hs.θ_dry += MH * th.θ_dry
    hs.θ_v += MH * th.θ_v
    hs.e_int += MH * th.e_int
    hs.h_m += MH * th.h_m
    hs.h_t += MH * th.h_t

    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        hs.ρq_tot += MH * state.moisture.ρq_tot
    end

    ν, D_t, _ = turbulence_tensors(atmos, state, diffusive_flux, aux, currtime)

    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        d_q_tot = (-D_t) .* diffusive_flux.moisture.∇q_tot
        hs.qt_sgs += MH * state.ρ * d_q_tot[end]
    end

    d_h_tot = -D_t .* diffusive_flux.∇h_tot
    hs.ht_sgs += MH * state.ρ * d_h_tot[end]

    # liquid water path
    # this condition is also going to be used to get the number of points that
    # exist on a horizontal plane provided all planes have the same number of
    # points
    # TODO adjust for possibility of non equivalent horizontal slabs
    if ev == floor(nvertelem / 2) && k == floor(Nqk / 2)
        # TODO: uncomment the line below after rewriting the LWP assignment below using aux.∫dz...?
        # aux = extract_aux(dg, localaux, ijk, e)
        LWP[1] += MH * (localaux[ijk, 1, e] + localaux[ijk, 2, e])
    end

    return nothing
end

function compute_diagnosticsums!(
    atmos,
    state,
    k,
    ijk,
    ev,
    e,
    Nqk,
    nvertelem,
    MH,
    thermoQ,
    horzavgs,
    dsums,
)
    zvals = CollectedDiagnostics.zvals
    th = thermo_vars(thermoQ[ijk, e])
    ha = horzavg_vars(horzavgs[Nqk * (ev - 1) + k])
    ds = diagnostic_vars(dsums[Nqk * (ev - 1) + k])

    u = state.ρu[1] / state.ρ
    v = state.ρu[2] / state.ρ
    w = state.ρu[3] / state.ρ
    ũ = ha.ρu / ha.ρ
    ṽ = ha.ρv / ha.ρ
    w̃ = ha.ρw / ha.ρ
    ẽ = ha.e_tot / ha.ρ
    q̃_tot = ha.ρq_tot / ha.ρ
    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        q_tot = state.moisture.ρq_tot / state.ρ
    end

    # state and functions of state
    ds.u += MH * ũ
    ds.v += MH * ṽ
    ds.w += MH * w̃
    ds.e_tot += MH * ẽ
    ds.q_tot += MH * ha.ρq_tot / ha.ρ
    ds.q_liq += MH * ha.q_liq
    ds.thd += MH * ha.θ_dry
    ds.thl += MH * ha.θ_liq_ice
    ds.thv += MH * ha.θ_v
    ds.e_int += MH * ha.e_int
    ds.h_m += MH * ha.h_m
    ds.h_t += MH * ha.h_t
    ds.qt_sgs += MH * ha.qt_sgs
    ds.ht_sgs += MH * ha.ht_sgs

    # vertical fluxes
    ds.vert_eddy_mass_flux += MH * (w - w̃) * (state.ρ - ha.ρ)
    ds.vert_eddy_u_flux += MH * (w - w̃) * (u - ha.ρu / ha.ρ)
    ds.vert_eddy_v_flux += MH * (w - w̃) * (v - ha.ρv / ha.ρ)
    ds.vert_eddy_ql_flux += MH * (w - w̃) * (th.q_liq - ha.q_liq)
    ds.vert_eddy_qv_flux += MH * (w - w̃) * (th.q_vap - ha.q_vap)
    ds.vert_eddy_thd_flux += MH * (w - w̃) * (th.θ_dry - ha.θ_dry)
    ds.vert_eddy_thv_flux += MH * (w - w̃) * (th.θ_v - ha.θ_v)
    ds.vert_eddy_thl_flux += MH * (w - w̃) * (th.θ_liq_ice - ha.θ_liq_ice)
    # TODO: temporary fix
    if isa(atmos.moisture, EquilMoist)
        ds.vert_eddy_qt_flux += MH * (w - w̃) * (q_tot - q̃_tot)
        ds.vert_qt_flux += MH * w * q_tot
    end

    # variances
    ds.uvariance += MH * (u - ũ)^2
    ds.vvariance += MH * (v - ṽ)^2
    ds.wvariance += MH * (w - w̃)^2

    # skewness
    ds.wskew += MH * (w - w̃)^3

    # turbulent kinetic energy
    ds.TKE = 0.5 * (ds.uvariance + ds.vvariance + ds.wvariance)

    return nothing
end

"""
    atmos_default_collect(bl, currtime)

Perform a global grid traversal to compute various diagnostics.
"""
function atmos_default_collect(dgngrp::DiagnosticsGroup, currtime)
    mpicomm = Settings.mpicomm
    dg = Settings.dg
    Q = Settings.Q
    mpirank = MPI.Comm_rank(mpicomm)
    current_time = string(currtime)

    # make sure this time step is not already recorded
    dprefix = @sprintf(
        "%s_%s-%s-num%04d",
        dgngrp.out_prefix,
        dgngrp.name,
        Settings.starttime,
        dgngrp.num
    )
    dfilename = joinpath(Settings.output_dir, dprefix)
    docollect = [false]
    if mpirank == 0
        dfullname = full_name(dgngrp.writer, dfilename)
        if isfile(dfullname)
            @warn """
Diagnostics $(dgngrp.name) collection
    output file $dfullname exists
    skipping collection at $current_time"""
        else
            docollect[1] = true
        end
    end
    MPI.Bcast!(docollect, 0, mpicomm)
    if !docollect[1]
        return nothing
    end

    # extract grid information
    bl = dg.balancelaw
    grid = dg.grid
    topology = grid.topology
    N = polynomialorder(grid)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    npoints = Nq * Nq * Nqk
    nrealelem = length(topology.realelems)
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    # get the state, auxiliary and geo variables onto the host if needed
    if Array ∈ typeof(Q).parameters
        localQ = Q.realdata
        localaux = dg.auxstate.realdata
        localvgeo = grid.vgeo
        localdiff = dg.diffstate.realdata
    else
        localQ = Array(Q.realdata)
        localaux = Array(dg.auxstate.realdata)
        localvgeo = Array(grid.vgeo)
        localdiff = Array(dg.diffstate.realdata)
    end
    FT = eltype(localQ)

    nstate = num_state(bl, FT)
    nauxstate = num_aux(bl, FT)
    ndiff = num_diffusive(bl, FT)

    # thermo variables
    thermoQ = [zeros(FT, num_thermo(FT)) for _ in 1:npoints, _ in 1:nrealelem]

    # horizontal sums and the liquid water path
    l_LWP = zeros(FT, 1)
    horzsums = [zeros(FT, num_horzavg(FT)) for _ in 1:(Nqk * nvertelem)]

    # compute thermo variables and horizontal sums in a single pass
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, localQ, ijk, e)
        aux = extract_aux(dg, localaux, ijk, e)

        compute_thermo!(bl, state, aux, ijk, e, thermoQ)

        diffusive_flux = extract_diffusion(dg, localdiff, ijk, e)
        MH = localvgeo[ijk, grid.MHid, e]
        compute_horzsums!(
            bl,
            state,
            diffusive_flux,
            aux,
            k,
            ijk,
            ev,
            e,
            Nqk,
            nvertelem,
            MH,
            localaux,
            thermoQ,
            horzsums,
            l_LWP,
            currtime,
        )
    end

    # compute the horizontal and LWP averages
    repdvsr = CollectedDiagnostics.repdvsr
    horzavgs = [zeros(FT, num_horzavg(FT)) for _ in 1:(Nqk * nvertelem)]
    for ev in 1:nvertelem
        for k in 1:Nqk
            hsum = MPI.Allreduce(horzsums[Nqk * (ev - 1) + k], +, mpicomm)
            horzavgs[Nqk * (ev - 1) + k] .= hsum ./ repdvsr[Nqk * (ev - 1) + k]
        end
    end
    LWP = zero(FT)
    LWP = MPI.Reduce(l_LWP[1], +, 0, mpicomm)
    if mpirank == 0
        LWP /= repdvsr[1]
    end

    # compute the diagnostics using the previous computed values
    dsums = [zeros(FT, num_diagnostic(FT)) for _ in 1:(Nqk * nvertelem)]
    @visitQ nhorzelem nvertelem Nqk Nq begin
        state = extract_state(dg, localQ, ijk, e)
        MH = localvgeo[ijk, grid.MHid, e]
        compute_diagnosticsums!(
            bl,
            state,
            k,
            ijk,
            ev,
            e,
            Nqk,
            nvertelem,
            MH,
            thermoQ,
            horzavgs,
            dsums,
        )
    end
    varvals = OrderedDict()
    varnames = flattenednames(vars_diagnostic(FT))
    for vari in 1:length(varnames)
        davg = zeros(FT, Nqk * nvertelem)
        for evk in 1:(Nqk * nvertelem)
            dsum = MPI.Reduce(dsums[evk][vari], +, 0, mpicomm)
            if mpirank == 0
                davg[evk] = dsum / repdvsr[evk]
            end
        end
        if mpirank == 0
            varvals[varnames[vari]] = davg
        end
    end

    # write diagnostics
    if mpirank == 0
        write_data(
            dgngrp.writer,
            dfilename,
            OrderedDict("z" => CollectedDiagnostics.zvals),
            varvals,
            currtime,
        )
    end

    # write LWP
    if mpirank == 0
        jldopen(
            joinpath(
                Settings.output_dir,
                "liquid_water_path-$(Settings.starttime).jld2",
            ),
            "a+",
        ) do file
            file[current_time] = LWP
        end
    end

    return nothing
end # function collect

function atmos_default_fini(dgngrp::DiagnosticsGroup, currtime) end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
@testset "modifiers.jl" begin
    @testset "reparametrized" begin
        # Interval by x->x²
        interval = Cube(1)
        ri = reparametrized(interval, x->x^2)
        @test @capture_out(print(ri)) == @capture_err(print(stderr, ri))
        @test ri(0.0) == SVector(0.00)
        @test ri(0.5) == SVector(0.25)
        @test ri(1.0) == SVector(1.00)
        r = rand()
        @test offsetframe(ri, r)[1] ≈ offsetframe(interval, r^2)[1]
        @test offsetframe(ri, r)[2] ≈ offsetframe(interval, r^2)[2]

        # Sphere by 2
        circle = Sphere(1)
        rc = reparametrized(circle, 2)
        @test @capture_out(print(rc)) == @capture_err(print(stderr, rc))
        @test rc(0.0) == rc(0.5) == rc(1.0)
        r = rand()
        @test offsetframe(rc, r)[1] ≈ offsetframe(circle, 2r)[1]
        @test offsetframe(rc, r)[2] ≈ offsetframe(circle, 2r)[2]
        @test norm(rand(rc)) ≈ 1

        # Identity reprametrization does not change the space.
        for i in 1:5
            @test reparametrized(Sphere(i), (args...) -> args) ≅ Sphere(i)
        end
    end

    @testset "scaled" begin
        #
        interval = Cube(1)
        si = scaled(interval, 2)
        @test @capture_out(print(si)) == @capture_err(print(stderr, si))
        @test si(0.0) == SVector(0.0)
        @test si(0.5) == SVector(1.0)
        @test si(1.0) == SVector(2.0)
        @test 2interval == si
        r = rand()
        @test offsetframe(si, r) == offsetframe(interval, 2r)
        @test scaled(interval, 2) == 2interval == interval * 2

        # Radius of sphere should match scaling.
        for i in 1:5
            r = rand()
            sc = r * Sphere(i)
            @test norm(rand(sc)) ≈ r
        end

        # Mulitplying by 1 doesn't change the space.
        for i in 1:5
            @test 1 * Sphere(i) ≅ Sphere(i)
        end
    end

    @testset "translated" begin
    end

    @testset "flattened" begin
        fs = flattened(Sphere(1))
        @test @capture_out(print(fs)) == @capture_err(print(stderr, fs))
        @test offsetframe(fs, 0.0)[1] ≈ SVector(1.0, 0.0)
        @test offsetframe(fs, 0.0)[2] ≈ reshape([0.0, 0.0, 1.0], (3, 1))
        @test fs ≅ Sphere(1)

        # Noop for spaces higher dimensions.
        for i in 2:5
            fi = flattened(Sphere(i))
            r = rand(i)
            @test offsetframe(fi, r...)[1] ≈ offsetframe(Sphere(i), r...)[1]
            @test offsetframe(fi, r...)[2] ≈ offsetframe(Sphere(i), r...)[2]
            @test fi ≅ Sphere(i)
        end
    end
end