text:
function getattribute(tiffdata::EzXML.Node, ::Type{Int}, attr::String)
    try
        parse(Int, tiffdata[attr])+1
    catch e
        (!isa(e, KeyError)) && rethrow(e)
        -1
    end
end

function getattribute(tiffdata::EzXML.Node, ::Type{Float64}, attr::String)
    try
        parse(Float64, tiffdata[attr])
    catch e
        (!isa(e, KeyError)) && rethrow(e)
        NaN
    end
end

"""
    $(SIGNATURES)

OMEXML is [very
flexible](https://docs.openmicroscopy.org/ome-model/6.0.1/ome-tiff/specification.html#fragment-1)
with its representation of the IFDs in the TIFF image. This function attempts to
handle many of the exceptions and update the passed collections with the proper
mapping of which TiffData elements correspond to which IFDs (and which files
these IFDs are located in) inside the TIFF image.

**Arguments**
- `ifd_index::OrderedDict{Int, NTuple{4, Int}}`: A mapping from IFD number to
  dimensions
- `ifd_files::OrderedDict{Int, Tuple{String, String}}`: A mapping from IFD
  number to the filepath and UUID of the file it's located in
- `obs_filepaths::Dict{String, Int}`: A list of observed filepaths mapped to the offset of their IFDs
- `image::EzXML.Node`: The OMEXML rooted at the current position
- `dims::NamedTuple`: Sizes of each dimension with the names as keys
- `tifffile::TiffFile`: A pointer to the root file
- `posidx::Int`: The index of the current position

The first two parameters should be then pumped through
[`OMETIFF.get_ifds`](@ref)
"""
function ifdindex!(ifd_index::OrderedDict{Int, NTuple{4, Int}},
                   ifd_files::OrderedDict{Int, Tuple{UUID, String}},
                   obs_filepaths::Dict{String, Int},
                   image::EzXML.Node,
                   dims::NamedTuple,
                   tifffile::TiffFile,
                   posidx::Int)

    tiffdatas = findall(".//ns:TiffData", image, ["ns"=>namespace(image)])

    ifd = 1
    # to avoid overwriting previous IFD indices from other images, we make sure
    # to start counting from the maximum previously observed IFD index.
    #
    # TODO: This assumes a dense numbering of the indices, is this always true?
    prev_ifd = length(ifd_index) > 0 ? maximum(keys(ifd_index)) : 0
    for tiffdata in tiffdatas
        try # if this tiffdata specifies the corresponding IFD
            ifd = parse(Int, tiffdata["IFD"]) + 1
        catch
            ifd = 1
        end

        uuid_node = findfirst("./ns:UUID", tiffdata, ["ns"=>namespace(tiffdata)])
        if uuid_node != nothing
            internal_filename = uuid_node["FileName"]

            uuid = nodecontent(uuid_node)
            filepath = joinpath(dirname(tifffile.filepath), internal_filename)

            # if this file isn't one we've observed before, increment the offset
            if !in(filepath, keys(obs_filepaths))
                obs_filepaths[filepath] = prev_ifd
            end

            ifd += obs_filepaths[filepath]

            ifd_files[ifd] = (make_uuid(uuid), filepath)
        end

        # get Z, C, T indices (in order specified by `dims`)
        indices = Tuple(getattribute(tiffdata, Int, "First$x") for x in keys(dims)[3:5])
        # how many ifds does this tiffdata correspond to
        p = getattribute(tiffdata, Int, "PlaneCount") - 1

        # if none of the Z, C, T indices are specified then we'll assume the
        # indices starting with the inner dimension, etc
        if all(indices .< 0)
            # index in the master ifd list
            idx = ifd
            # reverse iterate since we cycle the inner dimension the most
            for k=1:dims[5], j=1:dims[4], i=1:dims[3]
                ifd_index[idx] = (i, j, k, posidx)

                # if this tiffdata applies to multiple ifds then check that we
                # don't exceed the specified number of ifds
                (p > 1 && idx >= p+ifd-1) && break
                idx += 1
            end
        # if any of the indices are specified in the tiffdata then use these
        else
            indices = (indices..., posidx) # add the position index
            # all the indices that are not specified, we assume the first index
            ifd_index[ifd] = Tuple(pos > 0 ? pos : 1 for pos in indices)
        end
        prev_ifd = ifd
    end
end

"""
    get_unitful_axis(image, dimsize, stepsize, units) -> Range

Attempts to return a unitful range with a length of `dimsize`. Parameters
`stepsize` and `units` should be the XML tags in `image`.
"""
function get_unitful_axis(image::EzXML.Node, dimsize::Int, stepsize::String, units::String)
    try
        # OME-XML stores the step size
        increment = parse(Float64, image[stepsize])

        # parse units using Unitful
        unittype = uparse(image[units])

        # Create a unitful range
        return 0*unittype:increment*unittype:increment*(dimsize-1)*unittype
    catch
        # alternatively, just return the bare unitless range
        return 1:dimsize
    end
end

const axis_name_mapping = (X = :x, Y = :y, Z=:z, T=:time, C=:channel, P=:position)

"""
    build_axes(image) -> Tuple, Vector

Extracts the dimensions and axis information from the OMEXML data.

**Output**
- `NamedTuple{order, NTuple{6, Int}}`: the labeled 6 dimensions in the `order`
  that they are specified in the OMEXML.
- `Vector{AxisArray.Axis}`: List of `AxisArray.Axis` objects with units (if
  possible) in the same `order` as above

!!! warning
    There's no guarantee that the dimension sizes extracted here are correct if
    the acquisition was cancelled during a multiposition session. See
    [#38](https://github.com/tlnagy/OMETIFF.jl/issues/38). Downstream functions
    should be flexible and handle these cases.
"""
function build_axes(image::EzXML.Node)
    order = "YX"*join(replace(split(image["DimensionOrder"], ""), "X"=>"", "Y"=>""))
    order = Tuple(Symbol(dim) for dim in order)
    dims = NamedTuple{order, NTuple{5, Int}}(Tuple(parse(Int, image["Size$(x)"]) for x in order))
    dims = merge(dims, [:P=>1])

    # extract channel names
    channel_names = nodecontent.(findall("ns:Channel/@Name", image, ["ns"=>namespace(image)]))
    if isempty(channel_names)
        channel_names = ["C$x" for x in 1:dims[:C]]
    end

    vals = (
        X=get_unitful_axis(image, dims[:X], "PhysicalSizeX", "PhysicalSizeXUnit"),
        Y=get_unitful_axis(image, dims[:Y], "PhysicalSizeY", "PhysicalSizeYUnit"),
        Z=get_unitful_axis(image, dims[:Z], "PhysicalSizeZ", "PhysicalSizeZUnit"),
        C=to_symbol.(channel_names),
        T=get_unitful_axis(image, dims[:T], "TimeIncrement", "TimeIncrementUnit"),
        P=1
    )

    axes = [Axis{axis_name_mapping[key]}(vals[key]) for key in keys(dims)]

    dims, axes
end

"""
    $(SIGNATURES) -> AxisArray

Extracts the actual acquisition times from the OME-XML data. Takes `containers`, a vector of the XML nodes
corresponding to the root of each image.
"""
function get_elapsed_times(containers::Vector{EzXML.Node}, master_dims::NamedTuple, masteraxis::Vector{Axis}; default_unit=Unitful.s)

    # get the used dims that aren't x or y (since each xy plane has one elapsed time)
    nonxydims = NamedTuple{Tuple(k for (k, v) in pairs(master_dims) if k != :X && k != :Y && v > 1)}(master_dims)
    elapsed_times = fill(NaN*default_unit, values(nonxydims))

    # OMETIFF stores the Z, T, C info for each plane in TheZ, TheT, and TheC attributes
    attrnames = ["The$k" for k in keys(nonxydims) if k != :P]
    didx = fill(1, length(nonxydims))

    # we have to track position independently since it isn't stored in the plane tags
    # get which index corresponds to position (usually it's last)
    posdim = findfirst(x->x==:P, keys(nonxydims))
    unitstr = string(default_unit)

    for (pos, container) in enumerate(containers)
        planes = findall("ns:Plane[@DeltaT]", container, ["ns"=>namespace(container)])

        for plane in planes
            for d in 1:length(didx)
                # if this is the position dimension, get it from the current Pixel object since it's
                # not stored in the plane itself
                if d == posdim
                    didx[d] = pos
                else
                    didx[d] = getattribute(plane, Int, attrnames[d])
                end
            end
            try
                unitstr = plane["DeltaTUnit"]
            catch
            end
            unittype = uparse(unitstr)
            elapsed_times[didx...] = getattribute(plane, Float64, "DeltaT")*unittype
        end
    end
    used_axes = (masteraxis[i] for i in findall(x->x in keys(nonxydims), keys(master_dims)))
    AxisArray(elapsed_times, used_axes...)
end


text:
mutable struct Disjunct
    constraints::Vector{AbstractConstraint}
    indicator::VariableRef
end

mutable struct DisjunctionConstraint <: AbstractConstraint
    disjuncts::Vector{Disjunct}
end

mutable struct GDPdata
    gdp_variable_refs::Vector
    gdp_variable_names::Vector
    disjunctions::Dict
end

function GDPModel()
    model = Model()
    model.ext[:GDPdata] = GDPdata()
    # model.optimize_hook =

    return model
end

ext: jl


text:
const Vec = SVector
struct FontExtent{T}
    vertical_bearing::Vec{2, T}
    horizontal_bearing::Vec{2, T}

    advance::Vec{2, T}
    scale::Vec{2, T}
end

import Base: /, *, ==

import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, DefaultArrayStyle, materialize!, flatten, Style, broadcasted
BroadcastStyle(::Type{<: FontExtent}) = Style{FontExtent}()
BroadcastStyle(::Style{FontExtent}, x) = Style{FontExtent}()
BroadcastStyle(x, ::Style{FontExtent}) = Style{FontExtent}()

function broadcasted(op::Function, f::FontExtent, scaling::StaticVector)
    FontExtent(
        op.(f.vertical_bearing, scaling[1]),
        op.(f.horizontal_bearing, scaling[2]),
        op.(f.advance, scaling),
        op.(f.scale, scaling),
    )
end
function broadcasted(op::Function, ::Type{T}, f::FontExtent) where T
    FontExtent(
        map(x-> op(T, x), f.vertical_bearing),
        map(x-> op(T, x), f.horizontal_bearing),
        map(x-> op(T, x), f.advance),
        map(x-> op(T, x), f.scale),
    )
end

function FontExtent(fontmetric::FreeType.FT_Glyph_Metrics, scale::T = 64.0) where T <: AbstractFloat
    FontExtent(
        Vec{2, T}(fontmetric.vertBearingX, fontmetric.vertBearingY) ./ scale,
        Vec{2, T}(fontmetric.horiBearingX, fontmetric.horiBearingY) ./ scale,
        Vec{2, T}(fontmetric.horiAdvance, fontmetric.vertAdvance) ./ scale,
        Vec{2, T}(fontmetric.width, fontmetric.height) ./ scale
    )
end
function ==(x::FontExtent, y::FontExtent)
    x.vertical_bearing == y.vertical_bearing &&
    x.horizontal_bearing == y.horizontal_bearing &&
    x.advance == y.advance &&
    x.scale == y.scale
end
function FontExtent(fontmetric::FreeType.FT_Glyph_Metrics, scale::Integer)
    FontExtent(
        div.(Vec{2, Int}(fontmetric.vertBearingX, fontmetric.vertBearingY), scale),
        div.(Vec{2, Int}(fontmetric.horiBearingX, fontmetric.horiBearingY), scale),
        div.(Vec{2, Int}(fontmetric.horiAdvance, fontmetric.vertAdvance), scale),
        div.(Vec{2, Int}(fontmetric.width, fontmetric.height), scale)
    )
end

const FREE_FONT_LIBRARY = FT_Library[C_NULL]

function ft_init()
    global FREE_FONT_LIBRARY
    FREE_FONT_LIBRARY[1] != C_NULL && error("Freetype already initalized. init() called two times?")
    err = FT_Init_FreeType(FREE_FONT_LIBRARY)
    return err == 0
end

function ft_done()
    global FREE_FONT_LIBRARY
    FREE_FONT_LIBRARY[1] == C_NULL && error("Library == CNULL. FreeTypeAbstraction.done() called before init(), or done called two times?")
    err = FT_Done_FreeType(FREE_FONT_LIBRARY[1])
    FREE_FONT_LIBRARY[1] = C_NULL
    return err == 0
end

function newface(facename, faceindex::Real = 0, ftlib = FREE_FONT_LIBRARY)
    face = FT_Face[C_NULL]
    err = FT_New_Face(ftlib[1], facename, Int32(faceindex), face)
    if err != 0
        error("Couldn't load font $facename with error $err")
    end
    face
end

mutable struct FTFont
    ft_ptr::FreeType.FT_FaceRec
    function FTFont(ft_ptr::FreeType.FT_FaceRec)
        obj = new(ft_ptr)
        finalizer(FT_Done_Face, obj)
        return obj
    end
end

function FTFont(path::String)
    face = Ref{FT_Face}(C_NULL)
    err = FT_New_Face(ftlib[1], path, Int32(faceindex), face)
    if err != 0
        error("Couldn't load font $facename with error $err")
    end
    return FTFont(face[])
end

setpixelsize(face, x, y) = setpixelsize(face, (x, y))

function setpixelsize(face, size)
    err = FT_Set_Pixel_Sizes(face[1], UInt32(size[1]), UInt32(size[2]))
    if err != 0
    error("Couldn't set the pixel size for font with error $err")
    end
end

function kerning(c1::Char, c2::Char, face::Array{Ptr{FreeType.FT_FaceRec},1}, divisor::Float32)
    i1 = FT_Get_Char_Index(face[], c1)
    i2 = FT_Get_Char_Index(face[], c2)
    kernVec = Vector{FreeType.FT_Vector}(undef, 1)
    err = FT_Get_Kerning(face[], i1, i2, FreeType.FT_KERNING_DEFAULT, pointer(kernVec))
    err != 0 && return zero(Vec{2, Float32})
    return Vec{2, Float32}(kernVec[1].x / divisor, kernVec[1].y / divisor)
end

function loadchar(face, c::Char)
    err = FT_Load_Char(face[1], c, FT_LOAD_RENDER)
    @assert err == 0
end

function renderface(face, c::Char, pixelsize = (32,32))
    setpixelsize(face, pixelsize)
    faceRec = unsafe_load(face[1])
    loadchar(face, c)
    glyphRec    = unsafe_load(faceRec.glyph)
    @assert glyphRec.format == FreeType.FT_GLYPH_FORMAT_BITMAP
    return glyphbitmap(glyphRec.bitmap), FontExtent(glyphRec.metrics)
end

function getextent(face, c::Char, pixelsize)
    setpixelsize(face, pixelsize)
    faceRec = unsafe_load(face[1])
    loadchar(face, c)
    glyphRec = unsafe_load(faceRec.glyph)
    FontExtent(glyphRec.metrics)
end

function glyphbitmap(bmpRec::FreeType.FT_Bitmap)
    @assert bmpRec.pixel_mode == FreeType.FT_PIXEL_MODE_GRAY
    bmp = Matrix{UInt8}(undef, bmpRec.width, bmpRec.rows)
    row = bmpRec.buffer
    if bmpRec.pitch < 0
        row -= bmpRec.pitch * (rbmpRec.rows - 1)
    end

    for r = 1:bmpRec.rows
        srcArray = unsafe_wrap(Array, row, bmpRec.width)
        bmp[:, r] = srcArray
        row += bmpRec.pitch
    end
    return bmp
end

one_or_typemax(::Type{T}) where {T<:Union{Real,Colorant}} = T<:Integer ? typemax(T) : oneunit(T)

"""
    renderstring!(img::AbstractMatrix, str::String, face, pixelsize, y0, x0;
    fcolor=one_or_typemax(T), bcolor=zero(T), halign=:hleft, valign=:vbaseline) -> Matrix

Render `str` into `img` using the font `face` of size `pixelsize` at coordinates `y0,x0`.

# Arguments
* `y0,x0`: origin is in upper left with positive `y` going down
* `fcolor`: foreground color; typemax(T) for T<:Integer, otherwise one(T)
* `bcolor`: background color; set to `nothing` for transparent
* `halign`: :hleft, :hcenter, or :hright
* `valign`: :vtop, :vcenter, :vbaseline, or :vbottom
"""
function renderstring!(
        img::AbstractMatrix{T}, str::String, face, pixelsize, y0, x0;
        fcolor::T = one_or_typemax(T), bcolor::Union{T,Nothing} = zero(T),
        halign::Symbol = :hleft, valign::Symbol = :vbaseline
    ) where T<:Union{Real,Colorant}
    bitmaps = Vector{Matrix{UInt8}}(undef, lastindex(str))
    metrics = Vector{FontExtent{Int}}(undef, lastindex(str))
    ymin = ymax = sumadvancex = 0
    for (istr, char) = enumerate(str)
        bitmap, metric_float = renderface(face, char, pixelsize)
        metric = round.(Int, metric_float)
        bitmaps[istr] = bitmap
        metrics[istr] = metric
        w, h = metric.scale
        bx, by = metric.horizontal_bearing
        ymin = min(ymin, by - h)
        ymax = max(ymax, by)
        sumadvancex += metric.advance[1]
    end

    px = x0 - (halign == :hright ? sumadvancex : halign == :hcenter ? sumadvancex >> 1 : 0)
    py = y0 + (
        valign == :vtop ? ymax : valign == :vbottom ? ymin :
        valign == :vcenter ? (ymax - ymin) >> 1 + ymin : 0
    )
    bitmapmax = typemax(eltype(bitmaps[1]))

    imgh, imgw = size(img)
    if bcolor != nothing
        img[
            clamp(py-ymax+1, 1, imgh) : clamp(py-ymin-1, 1, imgh),
            clamp(px-1, 1, imgw) : clamp(px+sumadvancex-1, 1, imgw)
        ] .= bcolor
    end

    local prev_char::Char
    for (istr, char) = enumerate(str)
        w, h = metrics[istr].scale
        bx, by = metrics[istr].horizontal_bearing
        if istr == 1
            prev_char = char
        else
            kx, ky = map(x-> round(Int, x), kerning(prev_char, char, face, 64.0f0))
            px += kx
        end
        cliprowlo, cliprowhi = max(0, by-py), max(0, h+py-by-imgh)
        clipcollo, clipcolhi = max(0, bx-px), max(0, w+px-bx-imgw)
        if bcolor == nothing
            for row = 1+cliprowlo : h-cliprowhi, col = 1+clipcollo : w-clipcolhi
                bitmaps[istr][col,row]==0 && continue
                c1 = bitmaps[istr][col,row] / bitmapmax * fcolor
                img[row+py-by, col+px-bx] = T <: Integer ? round(T, c1) : T(c1)
            end
        else
            for row = 1+cliprowlo : h-cliprowhi, col = 1+clipcollo : w-clipcolhi
                bitmaps[istr][col, row] == 0 && continue
                w1 = bitmaps[istr][col, row] / bitmapmax
                c1 = w1 * fcolor
                c0 = (1.0 - w1) * bcolor
                img[row + py - by, col + px - bx] = T <: Integer ? round(T, c1 + c0) : T(c1 + c0)
            end
        end
        px += metrics[istr].advance[1]
    end
    img
end


text:
# Author: Soma S Dhavala

module domainUtils
# Functions defined in this file can be imported as modules

using JSON
using LightGraphs
using Requests
using DataStructures
using Distances
using JuMP

# given the file/url containing domain map extract all the concepts
# in the domain model - in particular read the relationships and  definitions
# create various graphs based on relationships and node attributes in the Domain Model

# Domain Model V2, 4th March, 2016
# Ref: https://github.com/ekstep/Learning-Platform/tree/master/docs/domain_model_v2

# read the domain model (concepts, relationships)
# returns status and the domain model object (Julia dictonary)

function readDomainModel(strname, src="url")
  status=false
  domainModelData = []
  try
    if(src=="url")
      resp = get(strname)
      if(statuscode(resp)==200)
        domainModelData = readall(resp)
      else
        error("domain model url response failed")
      end
    elseif(src=="file")
      if(!isfile(filename))
        error("File name specifying Domain Model is not a valid file")
      elseif (!isreadable(filename))
        error("File containing Domain Model is not readable")
      else
        # read the file containing the domain model
        f=open(filename)
        domainModelData=readall(f)
        close(f)
      end
    else
      error("Neither an url nor a file name is supplied to read Domain Model")
    end

    domainModelData=JSON.parse(domainModelData)

    # all the data is stored in the "result" key
    # concept definitions on are in "concepts" nested under "results"

    if(!haskey(domainModelData,"result"))
      error("invalid schema. \"result\" attribute not found.")
    elseif (!haskey(domainModelData["result",],"concepts"))
      error("invalid schema. \"concepts\" attribute not found.")
    end

    #  check concept definition
    dC = domainModelData["result",]["concepts",]

    if(!haskey(dC[1],"identifier"))
      error("invalid schema. \"result\" attribute not found. invalid schema")
    end

    #  check relation definition
    if (!haskey(domainModelData["result",],"relations"))
      error("invalid schema. \"relations\" attribute not found.")
    end

    # check relation schema
    dR = domainModelData["result",]["relations",]
    if (!haskey(dR[1],"startNodeId"))
      error("invalid schema. \"startNodeId\" attribute not found.")
    elseif (!haskey(dR[1],"endNodeId"))
      error("invalid schema. \"endNodeId\" attribute not found.")
    elseif (!haskey(dR[1],"relationType"))
      error("invalid schema. \"relationType\" attribute not found.")
    end

    status=true

  finally
    return status,domainModelData
  end
end

function getConceptsFromConceptMap(domainModelData)
  # read sorted unique list of concepts

  #  read all concepts
  dC = domainModelData["result",]["concepts",]
  concepts = []

  nC = length(dC)
  for i in 1:nC
    x = dC[i]
    ln = x["identifier",]
    #println(ln)
    push!(concepts,ln)
  end
  sort!(unique(concepts))

  # create a look-up from concept-to-index
  conceptHash = Dict()
  ind=1
  for c in concepts
    conceptHash[c]=ind
    ind=ind+1
  end

#  in-case s are preesent in the relationships, but not in the  definitions
#   # relationship dictionary
#   dR = domainModelData["result",]["relations",]
#   nR = length(dR)
#   for i in 1:nR
#     x = dR[i]
#     ln = x["startNodeId",]
#     rn = x["endNodeId",]
#     val = x["relationType",]
#     if val == "isParentOf"
#       push!(s,ln)
#       push!(s,rn)
#     end
#   end
#   s = union(s)
  return concepts, conceptHash
end

function getConceptGraph(domainModelData,relationType="isParentOf")

  # create a graph where nodes are concepts and edges satisfy "relationType"
  # property as defined in the Domain Model

  concepts, conceptHash = getConceptsFromConceptMap(domainModelData)
  n=length(concepts)
  cGraph=Array{Int64,2}(zeros(n,n))

  #  read relation data
  cDict = Dict()
  dR = domainModelData["result",]["relations",]
  nR = length(dR)

  for i in 1:nR
    x = dR[i]
    ln = x["startNodeId",]
    rn = x["endNodeId",]
    # interpreation of the relatiobship
    # startNodeId is parentOf of endNodeId
    val = x["relationType",]
    if val == relationType
      cDict[ln,rn] = 1
    end
  end

  for ci in concepts
    i=conceptHash[ci]
    for cj in concepts
      j=conceptHash[cj]
      if (haskey(cDict,(ci,cj)))
        cGraph[i,j]=1
      end
    end
  end
  return cGraph
end

function getGraphFromAdj(cGraph,directed=true)
  # create an LighGraph graph object from adj matrix
  n,m=size(cGraph)
  if(n!=m)
    error("graph has to be a square matrix")
  end

  if(directed)
    g = DiGraph(n)
  else
    g = Graph(n)
    # make the adj symmetric for undirected graph
    cGraph = (cGraph+cGraph')
    cGraph[cGraph.>1]=1
  end
  for i in 1:n
    for j in 1:n
      if (cGraph[i,j]==1)
        add_edge!(g,i,j)
      end
    end
  end
  return g
end

function getConceptSimilarityByRelation(cGraph,directed=true,dist="geodesic",weight=1,maxDepth=6)

  # given an graph (adj matrix), compute a concept similarity matrix (between concepts)
  g = getGraphFromAdj(cGraph*weight,directed)
  n = size(cGraph)[1]

  B=Array{Float64,2}(zeros(n,n))

  # compute shortes-distance between every pair of concepts (Dijkstra's shortest path/ Geodesic distance)
  for i in 1:n
    d = dijkstra_shortest_paths(g, i).dists
    d[d.>=maxDepth]=maxDepth
      B[i,:]=d
  end
  # turn "normalized" distance into a simialrity matrix
  B = 1 - (B/maximum(B))
  return B
end

function getConceptGradeCoverageVec(domainModelData,impute=true,directed=false,gradeTag="gradeLevel",weight=1)

  # grade is an attribute at concept level. each concept can have multple grades
  # prepare, for each concept, count of each grade (for concepts, the count can be either 0 or 1)
  # but at dimension level, they can vary

  concepts, conceptHash = getConceptsFromConceptMap(domainModelData)

  n=length(concepts)

  dC = domainModelData["result",]["concepts",]
  # get all the unique grade levels
  grades = []

  # in the 1st pass, look-up the grades to create unique grades levels
  for i in 1:n
    x = dC[i]
    if(haskey(x,"gradeLevel"))
      for val in x["gradeLevel"]
        push!(grades,val)
      end
    end
  end
  # need to figure out how to sort strings (crazy)
  #sort!(unique(grades))
  grades = (unique(grades))

  # create a hash for grades (grade-to-index)
  gradeHash = Dict()
  ind=1
  for grade in grades
    gradeHash[grade]=ind
    ind=ind+1
  end

  # 2nd pass
  # create  by grade count matrix
  # rows are s, cols are grade levels
  cGraph=Array{Int64,2}(zeros(n,ind-1))

  # default assume it is grade 1
  # debug: set grade count to zero
  cGraph[:,1] = 0

  for i in 1:n
    x = dC[i]
    c = x["identifier",]
    if(haskey(x,"gradeLevel"))
      vals = x["gradeLevel",]
        for val in vals
          cGraph[conceptHash[c],gradeHash[val]]+=1
        end
    end
  end

  # if grades need to be imputed at dimension and domain level (as of Domain Model v2)
  # using the parent-child relationship to aggregate the children's grade-count-vecs

  if(impute)
    rGraph = getConceptGraph(domainModelData)
    for c in concepts
      # node under consideration
      i = conceptHash[c]
      # get the parents of this list
      # parenet is any node from which a "i" (node of interest) can be reached
      # see if node i can be reached from node j, if yes, then "i" is a parent of "j"
      parList = find(rGraph[i,:])
      m = length(parList)
      if(m>0)
        for parent in parList
          cGraph[parent,:] += cGraph[i,:]
        end
      end
    end

  end
  return cGraph,grades,gradeHash,concepts,conceptHash
end

# obsolete
function getGradeVecSim(g1,g2,dist=Euclidean(),method="ordinal")
  # computes distance between two grade freq vectors
  n = length(g1)
  w = exp(0:n-1)
  w=w/sum(w)
  if(length(g2)!=n)
    error("vecs of not same size")
  end
  g1=vec(g1/sum(g1))
  g2=vec(g2/sum(g2))

  d = zeros(n,1)
  for i in 1:n
    d[i]=1*evaluate(CosineDist(),g1,circshift(g2,i-1))
  end
  return d
end

# obsolete
function getGradeDiffWeights(grades,method="ordinal")
  # park - need to develop
  m=length(grades)
  W=Array{Int64,2}(zeros(m,m))
  return W
end

function convertFloatToInt64(x,digits=2)
  # given a vec, normalize it, convert to an interger
  x = x/sum(x);
  M = 10^digits;
  x = round(x*10^digits);
  x = vec(map(tmp -> convert(Int64,round(tmp)),x));
  return x,M;
end

function getMinCostDistance(x,y,w,wext,verbose=true)
  # lables
  # x: source or supply
  # y: sink or demand
  n = length(x);
  if(n!=length(y))
    error("two vectors must of same size")
  end

  if((size(w)[1]!=size(w)[2]) | size(w)[1]!=n)
    error("cost matrix must be square")
  end
  # convert the vecs to integers
  x,M = convertFloatToInt64(x,2);
  y,M = convertFloatToInt64(y,2);

  # create exraneous source and sink nodes
  # in the case where supply is not equal to demand
  dif = sum(x)-sum(y);
  if(dif>0)
    # if supply is more, create extraneous demand node
    push!(x,0);
    push!(y,dif);
  else
    # if demand is more, create extraneous supply node
    push!(x,-dif);
    push!(y,0);
  end
  dif = sum(x)-sum(y);
  println("x: ",x)
  println("y: ",y)
  @assert dif == 0;

  N = n+1;
  M = sum(x);
  cost = Array{Int64,2}(zeros(N,N));
  # cost metric for given vec
  cost[1:n,1:n] = w;
  # cost to reach from ext.source to any sink node
  cost[N,1:n] = wext;
  # cost to reach ext.sink to any source node
  cost[1:n,N] = wext;

  mdl = Model();

  # flow is non-negative
  @defVar(mdl, flow[i=1:N, j=1:N] >= 0);
  # objective function:
  # minimize the overall flow from source to sink
  @setObjective(mdl, Min, sum{cost[i,j] * flow[i,j], i=1:N, j=1:N});
  # constraint: net out flow at source node should be equal its specified supply
  @addConstraint(mdl, xyconstr[i=1:1:N], sum{flow[i,j], j=1:N} == x[i]);
    # constraint: net out flow at source node should be equal its specified demand
  @addConstraint(mdl, xyconstr[j = 1:N], sum{flow[i,j], i=1:N} == y[j]);

  if(verbose)
    println("Solving original problem...")
  end
  status = solve(mdl);

  if status == :Optimal
    d = getObjectiveValue(mdl)/M/maximum(cost);
    flow = getValue(flow)
    if(verbose)
      @printf("Optimal!\n");
      @printf("Objective value: %d\n", d);
      @printf("Source:\n")
      println(x);
      @printf("Target:\n")
      println(y);
      @printf("Flow:\n");
      for j = 1:n
        @printf("\t%s", string("d",j));
      end
      @printf("\t%s", string("d.ext"));
      @printf("\n");
      for i = 1:n
        @printf("%s", string("s",i));
        for j = 1:N
          @printf("\t%d", (flow[i,j]));
        end
        @printf("\n");
      end
      @printf("%s", string("s.ext"));
      for j = 1:N
        @printf("\t%d", (flow[N,j]));
      end
      @printf("\n");
    end
  else
    d = [];
    if(verbose)
      @printf("No solution\n");
    end
  end
  return status, d, mdl, flow,cost;
end

function getConceptSimilarityByGrade(cGraph)
  # compute similarity between each pair of concepts (nodes) based on grade coverage
  # for future use: possibly consider that distances need to be weighed by grade levels
  # eg: grade-1 is farther from grade-5 than grade-1 is from grade-3
  d = pairwise(Euclidean(), cGraph)
  return(d)
end

end

max_stars_repo_name: harshavardhanc/sunbird-analytics

max_issues_repo_name: harshavardhanc/sunbird-analytics

max_forks_repo_name: harshavardhanc/sunbird-analytics


text:
### Optimize
# Needs https://github.com/JuliaLang/julia/pull/38287
import Core.Compiler: optimize
function optimize(interp::KernelInterpreter, opt::OptimizationState, params::OptimizationParams, @nospecialize(result))
    nargs = Int(opt.nargs) - 1
    ir = run_passes(opt.src, nargs, opt)
    Core.Compiler.finish(opt, params, ir, result)
end

import Core.Compiler: CodeInfo, convert_to_ircode, copy_exprargs, slot2reg, compact!, coverage_enabled, adce_pass!
import Core.Compiler: ssa_inlining_pass!, getfield_elim_pass!, type_lift_pass!, verify_ir, verify_linetable
function run_passes(ci::CodeInfo, nargs::Int, sv::OptimizationState)
    preserve_coverage = coverage_enabled(sv.mod)
    ir = convert_to_ircode(ci, copy_exprargs(ci.code), preserve_coverage, nargs, sv)
    ir = slot2reg(ir, ci, nargs, sv)
    # TODO: Domsorting can produce an updated domtree - no need to recompute here
    ir = compact!(ir)
    ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, ci.propagate_inbounds)
    ir = compact!(ir)
    ir = getfield_elim_pass!(ir) # SROA
    ir = adce_pass!(ir)
    ir = type_lift_pass!(ir)
    ir = compact!(ir)
    verify_ir(ir)
    verify_linetable(ir.linetable)
    return ir
end


text:
@testset "fexpr" begin

    fexpr = :( function bar(a, b::T1, c::T2) where {T1<:AbstractArray} where {T2<:Int64}
           print("vector method")
     end)
    f, args, wherestack, body = unwrap_fun(fexpr, true, true)
    fexpr1 = wrap_fun(f, args, wherestack, body)

    eval(fexpr1)

    methods(bar)

    # @test fexp1 == fexpr

end


text:
abstract type Model end
abstract type TerraceOnly <:Model end
abstract type TerraceKnickpoint <:Model end
abstract type TerraceDiffusion <:Model end
abstract type TerraceKnickpointDiffusion <:Model end

struct ModelOpts{T} end

function solve(solver::Type{ModelOpts{TerraceOnly}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer1, buffer2,
    ScratchTerrace_FEM, Scratch_FEM, to)

    # -- TERRACE PROFILE
    zt = Terrace.z
    # @timeit to "newton terrace" terracenewton!(
    #     TerracePhysics,
    #     (Δt * yr),
    #     zt,
    #     h_sea,
    #     id_shore.terrace,
    #     buffer1,
    #     buffer2)

    @timeit to "newton terrace" terracenewton2!(
        TerracePhysics,
        (Δt * yr),
        zt,
        h_sea,
        id_shore.terrace)

    return Terrace, River,to
end

function solve(solver::Type{ModelOpts{TerraceKnickpoint}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer1, buffer2,
    ScratchTerrace_FEM, Scratch_FEM, to)

    # -- RIVER PROFILE
    @timeit to "newton river" terracenewton!(
        TerracePhysics,
        Δt * yr,
        River.z,
        h_sea,
        id_shore.river + 1,
        buffer1,
        buffer2)

    # -- TERRACE PROFILE
    @timeit to "newton terrace" terracenewton!(
        TerracePhysics,
        Δt * yr,
        Terrace.z,
        h_sea,
        id_shore.terrace,
        buffer1,
        buffer2)

    # KNICKPOINT MIGRATION ADVECTION -> FLUX LIMITER TVD
    @timeit to "update r" updater!(RiverPhysics, RiverArrays, River)
    @timeit to "TVD" TVD!(River, buffer1, RiverPhysics, RiverArrays, Δt, id_shore.river)

    return Terrace, River,to
end

function solve(solver::Type{ModelOpts{TerraceKnickpointDiffusion}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer1, buffer2,
    ScratchTerrace_FEM, Scratch_FEM, to)

   # -- RIVER PROFILE
    @timeit to "newton river" terracenewton!(
        TerracePhysics,
        Δt * yr,
        River.z,
        h_sea,
        id_shore.river + 1,
        buffer1,
        buffer2)

    # -- TERRACE PROFILE
    @timeit to "newton terrace" terracenewton!(
        TerracePhysics,
        Δt * yr,
        Terrace.z,
        h_sea,
        id_shore.terrace,
        buffer1,
        buffer2)

    # KNICKPOINT MIGRATION ADVECTION -> FLUX LIMITER TVD
    @timeit to "update r" updater!(RiverPhysics, RiverArrays, River)
    @timeit to "TVD" TVD!(River, buffer1, RiverPhysics, RiverArrays, Δt, id_shore.river)

    # Diffusion over Δt (FEM)
    @timeit to "FEM River" River.z = femsolver(Δt, Δx, nn, River, h_sea,
                                               Scratch_FEM,
                                               HillSlopePhysics,
                                               SubmarinePhysics)

    # TERRACE
    # Diffusion over Δt
    @timeit to "FEM Terrace" Terrace.z = femsolver(
            Δt,
            Terrace,
            h_sea,
            ScratchTerrace_FEM,
            Scratch_FEM,
            HillSlopePhysics
    )

    return Terrace, River,to
end

function solve(solver::Type{ModelOpts{TerraceDiffusion}},
    TerracePhysics, RiverPhysics, RiverArrays, HillSlopePhysics, SubmarinePhysics,
    Δt, River, Terrace, h_sea, id_shore, buffer1, buffer2,
    ScratchTerrace_FEM, Scratch_FEM, to)

    # -- TERRACE PROFILE
    @timeit to "newton terrace" terracenewton!(
        TerracePhysics,
        (Δt * yr),
        Terrace.z,
        h_sea,
        (id_shore.terrace),
        buffer1,
        buffer2)

    # TERRACE
    @timeit to "FEM Terrace" Terrace.z = femsolver(
            Δt,
            Terrace,
            h_sea,
            ScratchTerrace_FEM,
            Scratch_FEM,
            HillSlopePhysics
    )

    return Terrace, River,to
end

function old_solver!(Terrace, TerracePhysics, h_sea, id_shore, dt)
    P, h_wb, βz = TerracePhysics.P0, TerracePhysics.h_wb, TerracePhysics.βz
    for i in id_shore.terrace.i1: id_shore.terrace.i2
        @inbounds Terrace.z[i] += - βz .* P .* exp(-4*(h_sea-Terrace.z[i])/h_wb) * dt
    end
end


text:
using Vlasiator, LaTeXStrings, SHA, LazyArtifacts
using Suppressor: @capture_out, @capture_err, @suppress_out
using REPL.TerminalMenus
using Test

group = get(ENV, "TEST_GROUP", :all) |> Symbol

function nanmaximum(x::AbstractArray{T}) where T<:AbstractFloat
   result = convert(eltype(x), NaN)
   for i in x
      if !isnan(i)
         if isnan(result) || i > result
            result = i
         end
      end
   end
   result
end

function simulateInput(keys...)

   keydict = Dict(:up => "\e[A",
                  :down => "\e[B",
                  :enter => "\r")

   for key in keys
      if isa(key, Symbol)
         write(Base.stdin.buffer, keydict[key])
      else
         write(Base.stdin.buffer, "$key")
      end
   end
end

@testset "Vlasiator.jl" begin
   rootpath = artifact"testdata"

   files = joinpath.(rootpath, ("bulk.1d.vlsv", "bulk.2d.vlsv", "bulk.amr.vlsv"))
   meta1 = load(files[1])
   meta2 = load(files[2])
   meta3 = load(files[3])

   if group in (:read, :all)
      @testset "Reading files" begin
         @test_throws ArgumentError load("data")
         meta = meta1
         @test ndims(meta) == 1
         @test startswith(repr(meta), "File: bulk.1d.vlsv")
         @test size(meta) == 529201
         # Variable strings reading
         @test meta.variable[end] == "vg_boundarytype"
         # Variable info reading
         varinfo = readvariablemeta(meta, "proton/vg_rho")
         @test startswith(repr(varinfo), "Variable in LaTeX")
         @test varinfo.unit == "1/m^3"
         # Velocity mesh display
         @test startswith(repr(meta.meshes["proton"]), "vblocks")
         # Parameter checking
         @test hasparameter(meta, "dt")
         # Parameter reading
         t = readparameter(meta, "time")
         @test t == 10.0
         @test_throws ArgumentError meta["nonsense"]
         # Do-Block syntax
         t = load(files[1]) do meta
            readparameter(meta, "time")
         end
         @test t == 10.0
         # unsorted ID
         @test readvariable(meta, "CellID", false) == 10:-1:1
         indexRef = [10, 9, 8, 7, 6, 5, 4, 3, 2, 1]
         @test meta.cellindex == indexRef
         # sorted var by default
         @test meta["vg_boundarytype"] == [4, 4, 1, 1, 1, 1, 1, 1, 3, 3]
         # ID finding (noAMR)
         loc = [2.0, 0.0, 0.0]
         id = getcell(meta, loc)
         coords = getcellcoordinates(meta, id)
         @test coords == (3.0, 0.0, 0.0)
         @test readvariable(meta, "proton/vg_rho", id)[1] == 1.2288102f0
         # ID in a line
         point1 = [-4.0, 0.0, 0.0]
         point2 = [4.0, 0.0, 0.0]
         cellids, _, _ = getcellinline(meta, point1, point2)
         @test cellids == 4:7
         point1 = [-10.1, 0.0, 0.0]
         @test_throws DomainError getcellinline(meta, point1, point2)
         point2 = [10.1, 0.0, 0.0]
         @test_throws DomainError getcellinline(meta, point1, point2)

         # Nearest ID with VDF stored
         @test getnearestcellwithvdf(meta, id) == 5

         # velocity space reading
         vcellids, vcellf = readvcells(meta, 5; species="proton")
         V = getvcellcoordinates(meta, vcellids; species="proton")
         @test V[end] == [2.45f0, 1.95f0, 1.95f0]
         @test_throws ArgumentError readvcells(meta, 2)
         vcids = Vlasiator.reorder(meta.meshes["proton"], vcellids)
         @test vcids[5] == 0x00000029
         f = Vlasiator.reconstruct(meta.meshes["proton"], vcellids, vcellf)
         @test f[CartesianIndex(26, 20, 20)] == 85.41775f0
         @test getdensity(meta, f) ≈ 1.8255334f0
         @test getdensity(meta, vcellf) ≈ 1.8255334f0
         @test getvelocity(meta, f)[1] ≈ 1.0f0 rtol=3e-3
         @test getvelocity(meta, vcellids, vcellf)[1] ≈ 1.0f0 rtol=3e-3
         @test getpressure(meta, f) ≈ zeros(Float32, 6) atol=1e-16
         @test getpressure(meta, vcellids, vcellf) ≈ zeros(Float32, 6) atol=1e-16
         @test getmaxwellianity(meta, f) ≈ 5.741325243685855 rtol=1e-4
         @test getmaxwellianity(meta, vcellids, vcellf) ≈ 5.741325243685855 rtol=1e-4

         # AMR data reading, DCCRG grid
         metaAMR = meta3
         sliceoffset = abs(metaAMR.coordmin[2])
         idlist, indexlist = getslicecell(metaAMR, sliceoffset, 2,
            metaAMR.coordmin[2], metaAMR.coordmax[2])

         # ID finding (AMR)
         loc = (2.5e6, 2.5e6, 2.5e6) # exact cell center
         id = getcell(metaAMR, loc)
         @test getcellcoordinates(metaAMR, id) == loc

         data = readvariable(metaAMR, "proton/vg_rho")
         dataslice = refineslice(metaAMR, idlist, data[indexlist], :y)
         @test sum(dataslice) ≈ 7.690352275026747e8
         let err = nothing
            try
               getslicecell(metaAMR, sliceoffset, 1, -2., -1.)
            catch err
            end
            @test err isa Exception
         end
         # AMR ID finding
         loc = [12*Vlasiator.RE, 0.0, 0.0]
         @test getcell(metaAMR, loc) == 0x00000000000002d0

         # AMR level
         @test metaAMR.maxamr == 2
         @test getlevel(metaAMR, idlist[1]) == 1

         # DCCRG utilities
         @test isparent(metaAMR, 1)
         @test !isparent(metaAMR, 1080)
         @test getchildren(metaAMR, 1) == getsiblings(metaAMR, 129)
         @test getparent(metaAMR, 129) == 1
         @test_throws ArgumentError getparent(metaAMR, 5)
         @test_throws ArgumentError getsiblings(metaAMR, 1)

         # FS grid
         data = readvariable(metaAMR, "fg_e")
         ncells, namr = metaAMR.ncells, metaAMR.maxamr
         @test size(data) == (3, ncells[1]*namr^2, ncells[2]*namr^2, ncells[3]*namr^2)
         @test data[:,16,8,8] == [-1.3758785f-7, 3.2213068f-4, -3.1518404f-4]

         # Compare two VLSV files
         @test issame(files[1], files[1])
      end
   end

   if group in (:derive, :all)
      @testset "Derived variables" begin
         meta = meta1
         @test meta["Vmag"] |> sortperm == [3, 4, 1, 2, 8, 9, 10, 5, 7, 6]
         @test meta["Protated"][2,2,3] == 1.7952461f-29
         @test meta["Panisotropy"][4] == 1.0147696f0
         @test meta["Tanisotropy"][4] == 1.0147695404826917
         @test meta["Agyrotropy"][1] |> isnan

         meta = meta2
         @test meta["Bmag"][4] == 3.0052159f-9

         @test meta["Emag"][1,10,99] ≈ 2.6120074f-6

         @test meta["VS"] |> nanmaximum == 1.3726345956957596e6

         @test meta["VA"] |> nanmaximum == 2.3202628822256166e8

         @test readvariable(meta, "VA", [1,2])[2] == 65507.75496283364

         @test meta["MA"][end] == 10.700530839822328

         @test meta["MS"][end] == 16.9375888861409

         @test meta["Vpar"][1] == 698735.3f0

         @test meta["Vperp"][1] == 40982.48f0

         @test_throws ArgumentError meta["Epar"]

         @test_throws ArgumentError meta["Eperp"]

         @test meta["T"][1] == 347619.9817319378

         @test meta["Pram"][1] == 8.204415428337215e-10

         @test meta["Pb"][1] == 3.59502523919876e-12

         @test meta["Beta"][1] == 1.3359065984817116

         @test meta["BetaStar"][1] == 229.55170154977864

         @test meta["Poynting"][:,10,10] == [-3.677613f-11, 8.859047f-9, 2.4681486f-9]

         @test meta["IonInertial"][1] == 5.389771470423157e8

         @test readvariable(meta, "Larmor", UInt64[1])[1] == 322324.70603759587

         @test meta["Gyroperiod"][1] == 21.834297799454554

         @test Vlasiator.getdata2d(meta, "Gyroperiod") |> ndims == 2

         @test meta["Gyrofrequency"][1] == 0.04579950356933307

         @test meta["J"][1,1000] == 1.4322914301940016e-15

         @test meta["Omegap"][1] == 209.5467447842415

         @test meta["Plasmaperiod"][1] == 0.0047722048893178645

         @test meta["MagneticTension"][3, 40, 50] == -3.6856588f-19
      end
      @testset "VLSV writing" begin
         meta = meta1
         # Obtain unsorted derived variables, workaround #59
         vmag = readvariable(meta, "Vmag", meta.cellid)
         pa = readvariable(meta, "Panisotropy", meta.cellid)
         vars = Vector{Tuple{VecOrMat, String, VarInfo}}(undef, 0)
         push!(vars, (vmag, "vmag", VarInfo("m/s", L"$\mathrm{m}/mathrm{s}$", L"$V$", "")))
         push!(vars, (pa, "panisotropy", VarInfo("", "", "", "")))

         write_vlsv(files[1], "bulk_new.vlsv", vars)
         sha_str = bytes2hex(open(sha1, "bulk_new.vlsv"))
         @test sha_str == "2b209fb022a2db3b013e01606a60c1fac360a79d"

         rm("bulk_new.vlsv", force=true)
      end
   end

   if group in (:utility, :all)
      @testset "Rotation" begin
         e1 = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
         e2 = [0.0 1.0 0.0; 1.0 0.0 0.0; 0.0 0.0 1.0]
         R = Vlasiator.getRotationMatrix(e1, e2)
         @test R == e2

         v = fill(1/√3, 3)
         θ = π / 4
         R = Vlasiator.getRotationMatrix(v, θ)
         @test inv(R) ≈ R'

         Rᵀ = Vlasiator.rotateTensorToVectorZ(R, v)
         @test Rᵀ ≈ [1/√2 -1/√2 0.0; 1/√2 1/√2 0.0; 0.0 0.0 1.0]
      end
      @testset "Curvature, Divergence" begin
         let dx = ones(Float32, 3)
            # 2D
            A = ones(Float32, 3,3,1,3)
            @test sum(Vlasiator.curl(A, dx)) == 0.0
            # 3D
            A = ones(Float32, 3,3,3,3)
            @test sum(Vlasiator.curl(A, dx)) == 0.0
            @test sum(Vlasiator.divergence(A, dx)) == 0.0
         end
      end
      @testset "Gradient" begin
         let A = ones(Float32, 3,3,3), dx = ones(Float32, 3)
            @test sum(Vlasiator.gradient(A, dx)) == 0.0
         end
         let A = ones(Float32, 3,3)
            @test sum(Vlasiator.gradient(A)) == 0.0
         end
         let A = ones(Float32, 3)
            @test sum(Vlasiator.gradient(A)) == 0.0
         end
      end
      @testset "FluxFunction" begin
         b = reshape(1:75, 3, 5, 5)
         dx = [1,2]
         flux = compute_flux_function(b, dx, 1)
         @test flux[3,3] == 29
         # saddle point test func
         flux = [x^2 - y^2 for x in -10:1.0:10, y in -10:1.0:10]
         xi_, oi_ = find_reconnection_points(flux; method=1)
         @test xi_ == [11; 11;;]
         xi_, oi_ = find_reconnection_points(flux; method=2)
         @test xi_ == [11; 11;;]
      end
   end

   if group in (:vtk, :all)
      @testset "VTK" begin
         meta = meta2 # no amr
         data, ghostType = Vlasiator.fillmesh(meta, ["proton/vg_rho"])
         @test size(data[1][1]) == (1, 63, 100, 1)

         meta = meta3 # amr
         write_vtk(meta, vars=["proton/vg_rho", "fg_b", "proton/vg_v"])
         sha_str = bytes2hex(open(sha1, "bulk.amr_3.vti"))
         @test sha_str == "8a2bb0a15c5dcc329f88821036df840a86eef9d5"

         filesaved = ["bulk.amr.vthb", "bulk.amr_1.vti", "bulk.amr_2.vti", "bulk.amr_3.vti"]
         rm.(filesaved, force=true)
      end
   end

   if group in (:log, :all)
      @testset "Log" begin
         file = joinpath(rootpath, "logfile.txt")
         timestamps, speed = readlog(file)
         @test length(speed) == 50 && speed[end] == 631.2511f0
      end
   end

   if group in (:monitor, :all)
      @testset "Monitor" begin
         output = @capture_out begin
            n = 2e6    # [amu/m³]
            v = 6e5    # [m/s]
            T = 5e5    # [K]
            B = 5e-9   # [T]
            check_plasma_characteristics(n, v, T, B)
         end
         @test startswith(output, "---")
      end
   end

   if group in (:plot, :all)
      @testset "PyPlot" begin
         using PyPlot
         ENV["MPLBACKEND"]="agg" # no GUI
         # 1D
         meta = meta1
         line = plot(meta, "proton/vg_rho")[1]
         @test line.get_ydata() == meta["proton/vg_rho"]
         line = plot(meta, "proton/vg_v", comp=:3)[1]
         @test line.get_ydata() == meta["proton/vg_v"][3,:]
         centers = plotmesh(meta, projection="y")
         points = centers.get_offsets()
         @test size(points) == (10, 2)
         fig = plt.figure()
         ax = fig.add_subplot(projection="3d")
         centers = plotmesh(meta, projection="3d")
         points = centers.get_offsets() # only 2D from 3D coords, might be improved
         @test size(points) == (10, 2)
         close(fig)

         @test_throws ArgumentError pcolormesh(meta, "proton/vg_rho")

         loc = [2.0, 0.0, 0.0]
         var = vdfslice(meta, loc).get_array()
         @test var[786] == 238.24398578141802
         @test_throws ArgumentError vdfslice(meta, loc, species="helium")
         output = @capture_err begin
            var = vdfslice(meta, loc; slicetype=:bperp).get_array()
         end
         @test var[786] == 4.02741885708042e-10

         output = @capture_err begin
            vdfslice(meta, loc; verbose=true)
         end
         @test startswith(output, "[ Info:")

         # 2D
         meta = meta2
         var = pcolormesh(meta, "proton/vg_rho").get_array()
         @test var[end-2] == 999535.8f0 && length(var) == 6300
         var = pcolormesh(meta, "proton/vg_rho", extent=[0,1,0,2]).get_array()
         @test var[end] == 0.0 && length(var) == 6
         var = pcolormesh(meta, "proton/vg_rho";
            extent=[-2e7, 2e7, -2e7, 2e7], axisunit=SI, colorscale=Log).get_array()
         @test var[end] == 1.00022675f6 && length(var) == 100
         var = pcolormesh(meta, "fg_b").get_array()
         @test var[1] == 3.0058909f-9
         var = pcolormesh(meta, "proton/vg_v", comp=:x, colorscale=SymLog).get_array()
         @test var[2] == -699935.2f0
         p = streamplot(meta, "proton/vg_v", comp="xy")
         @test typeof(p) == PyPlot.PyObject
         var = quiver(meta, "proton/vg_v", axisunit=SI, stride=1).get_offsets()
         @test size(var) == (6300, 2)

         # 3D AMR
         meta = meta3
         var = pcolormesh(meta, "proton/vg_rho").get_array()
         @test var[255] == 1.0483886f6 && length(var) == 512
         var = pcolormesh(meta, "proton/vg_v").get_array()
         @test var[255] == 99992.586f0
         @test_throws ArgumentError pcolormesh(meta, "fg_b")
      end

      @testset "Plots" begin
         isdefined(Main, :PyPlot) && import Vlasiator.vdfslice
         include("../src/plot/plots.jl")
         RecipesBase.is_key_supported(k::Symbol) = true
         # 1D
         meta = meta1
         rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), meta, "proton/vg_rho")
         @test getfield(rec[1], 1)[:seriestype] == :line &&
            rec[1].args[1] isa LinRange

         rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(),
            VDFSlice((meta, [0.0,0.0,0.0])))
         @test getfield(rec[1], 1)[:seriestype] == :histogram2d

         # 2D
         meta = meta2
         rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), meta, "proton/vg_rho")
         @test getfield(rec[1], 1)[:seriestype] == :heatmap &&
            rec[1].args[1] isa LinRange
      end

      @testset "Plot UI" begin
         meta = meta1
         simulateInput(:enter)
         @suppress_out pui(meta; suppress_output=true)
         @test true # if it reaches here, then pass
         meta = meta2
         simulateInput(fill(:down,5), :enter, :enter);
         @suppress_out pui(meta; suppress_output=true)
         @test true # if it reaches here, then pass
      end
   end
   for meta in (meta1, meta2, meta3)
      close(meta.fid) # required for Windows?
   end
end


text:
using Koala
using KoalaFlux
using Flux
using Test
import DataFrames

X, y = load_ames()
@test KoalaFlux.ordinal_categorical_features(X) ==
    (Symbol[:OverallQual, :GrLivArea, :x1stFlrSF, :TotalBsmtSF,
            :BsmtFinSF1, :LotArea, :GarageCars, :GarageArea,
            :YearRemodAdd, :YearBuilt],
     Symbol[:Neighborhood, :MSSubClass])

t = KoalaFlux.FrameToFluxInputTransformer()
tM = Machine(t, X)
XX = transform(tM, X[2:6,:])
@test XX.categorical[2,1] ==  XX.categorical[2,3]
@test XX.categorical[1,1] != XX.categorical[1,2]

for ftr in tM.scheme.categorical_features
    delete!(X, ftr)
end
t = KoalaFlux.FrameToFluxInputTransformer()
tM = Machine(t, X)
XX = transform(tM, X[2:6,:])

# 3D embedding for a two-class categorical:
A1 = Float64[1 2;
             3 4;
             0 5]

# 2D embedding for a three-class categorical:
A2 = Float64[1 2 3;
             0 1 5]

embedding = [A1, A2]
@test KoalaFlux.x(embedding, [1.2, 2.3, 4.3, 5.6], [2, 1]) ==
    [1.2, 2.3, 4.3, 5.6, 2.0, 4.0, 5.0, 1.0, 0.0]

X, y = load_ames()
train, test = partition(eachindex(y), 0.7)
flux = FluxRegressor()
flux.n = 3
flux.learning_rate = 0.002
fluxM = Machine(flux, X, y, train)
fit!(fluxM, train)

function get_creator(p)
    function creator(n_inputs)
        n_hidden = round(Int, sqrt(n_inputs))
        return Flux.Chain(Flux.Dense(n_inputs, n_hidden, Flux.σ),
                          Flux.Dropout(p), Flux.Dense(n_hidden, 1))
    end
    return creator
end

flux.network_creator = get_creator(0.5)
learning_curve(fluxM, train, test, [2,4,6])

flux.embedding0 = fluxM.report[:embedding]
flux.chain0 = fluxM.report[:chain]
fit!(fluxM)

# tests for CategoricalEmbedder:

letters = ['a', 'b', 'a']
weirds = ["goo", "goo", "prickles"]
numbers = [1.0, 2.0, 1.0]
yy = [-1.0, 0.0, 1.0]

XX = DataFrames.DataFrame(letters=letters, numbers=numbers, weirds=weirds)

flux = FluxRegressor()
flux.n = 1
flux.dimension_formula = x -> x
flux.learning_rate = 0.0
fluxM = Machine(flux, XX, yy, 1:3)
fit!(fluxM, 1:3)
afloat, bfloat = fluxM.report[:embedding][1].data[1,:]
t = CategoricalEmbedder(fluxM)
tM = Machine(t, XX)
s = tM.scheme
@test transform(s, :letters, 1, letters) ≈ [afloat, bfloat, afloat]
transform(tM, XX)
@test transform(tM, XX[[:letters]])[:letters__1] == transform(tM, XX)[:letters__1]


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

script = """
# Issue #11948
f(x) = x+1
workspace()
@assert !isdefined(:f)
LastMain.f(2)

# PR #12990
io = IOBuffer()
show(io, Pair)
@assert String(take!(io)) == "Pair"
@assert !Base.inbase(LastMain)
"""
exename = Base.julia_cmd()
run(`$exename --startup-file=no -e $script`)

# issue #17764
script2 = """
mutable struct Foo end
workspace()
mutable struct Foo end
@assert Tuple{Type{LastMain.Foo}} !== Tuple{Type{Main.Foo}}
"""
run(`$exename --startup-file=no -e $script2`)


text:
"""
Bridges `CP.Sort` to MILP constraints by adding O(n²) binary variables, with a
transportation-like model.

## Detailed model

Let `x` be the array to sort and `y` its sorted copy, both vectors having
length `n`. This bridge handles the constraint `[y..., x...]`-in-`CP.Sort(n)`.

Two sets of variables:
* `f[i, j]`: real variable, `i` from 1 to `n`; the "flow" from the array to
  sort to the sorted copy, equal to the value of `x[i]` and `y[j]`
* `a[i, j]`: binary variable, `i` from 1 to `n`; `a[i, j]` indicates whether
  the flow `f[i, j]` is nonzero

Constraints:
* Flow coming from the array to sort `x`:
    ``x_i = \\sum_{j=1}^n f_{i,j} \\qquad \\forall i \\in \\{1, 2\\dots n\\}``
* Flow going to the sorted array:
    ``y_j = \\sum_{i=1}^n f_{i,j} \\qquad \\forall j \\in \\{1, 2\\dots n\\}``
* The flow from one value of the array to sort can only go to one element of
  the sorted array:
    ``\\sum_{i=1}^n a_{i,j} = 1 \\qquad \\forall j \\in \\{1, 2\\dots n\\}``
* The flow to one value of the sorted array can only come from one element of
  the array to sort:
    ``\\sum_{j=1}^n a_{i,j} = 1 \\qquad \\forall i \\in \\{1, 2\\dots n\\}``
* The flow `f[i, j]` is related to the binary variable `a[i, j]`:
    ``U a_{i,j} \\leq f_{i,j} \\leq L a_{i,j} \\qquad \\forall (i,j) \\in \\{1, 2\\dots n\\}^2``
* The array `y` must be sorted:
    ``y_{j-1} \\leq y_{j} \\qquad \\forall j \\in \\{2, 3\\dots n\\}``
"""
struct Sort2MILPBridge{T} <: MOIBC.AbstractBridge
    vars_flow::Matrix{MOI.VariableIndex}
    vars_unicity::Matrix{MOI.VariableIndex}
    vars_unicity_bin::Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}
    cons_transportation_x::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_transportation_y::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_unicity_x::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_unicity_y::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}
    cons_flow_gt::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}
    cons_flow_lt::Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    cons_sort::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}
end

function MOIBC.bridge_constraint(
    ::Type{Sort2MILPBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.Sort,
) where {T}
    return MOIBC.bridge_constraint(
        Sort2MILPBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{Sort2MILPBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.Sort,
) where {T}
    f_scalars = MOIU.scalarize(f)
    dim = s.dimension

    f_sorted = f_scalars[1:dim]
    f_to_sort = f_scalars[(dim + 1):(2 * dim)]

    # For this formulation work, both lower and upper bounds are required on
    # the values of the array to sort.
    for i in 1:dim
        @assert CP.has_lower_bound(model, f_to_sort[i])
        @assert CP.has_upper_bound(model, f_to_sort[i])
    end

    # Create the new variables.
    vars_flow = Matrix{MOI.VariableIndex}(undef, dim, dim)
    vars_unicity = Matrix{MOI.VariableIndex}(undef, dim, dim)
    vars_unicity_bin = Matrix{MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}}(undef, dim, dim)
    for i in 1:dim
        vars_flow[i, :] = MOI.add_variables(model, dim)
        vars_unicity[i, :], vars_unicity_bin[i, :] = MOI.add_constrained_variables(model, [MOI.ZeroOne() for _ in 1:dim])
    end

    # Write the transportation equations.
    cons_transportation_x = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_sorted[i] - sum(one(T) .* vars_flow[i, :]),
            MOI.EqualTo(zero(T))
        )
        for i in 1:dim
    ]
    cons_transportation_y = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            f_to_sort[j] - sum(one(T) .* vars_flow[:, j]),
            MOI.EqualTo(zero(T))
        )
        for j in 1:dim
    ]

    # Unicity constraints, so that only one input and one output contribute to each flow.
    cons_unicity_x = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            sum(one(T) .* vars_unicity[i, :]),
            MOI.EqualTo(one(T))
        )
        for i in 1:dim
    ]
    cons_unicity_y = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}[
        MOI.add_constraint(
            model,
            sum(one(T) .* vars_unicity[:, j]),
            MOI.EqualTo(one(T))
        )
        for j in 1:dim
    ]

    # Relate flows and unicity.
    cons_flow_gt = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}}(undef, dim, dim)
    cons_flow_lt = Matrix{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}(undef, dim, dim)
    for i in 1:dim
        L = T(CP.get_lower_bound(model, f_to_sort[i]))
        U = T(CP.get_upper_bound(model, f_to_sort[i]))

        for j in 1:dim
            cons_flow_gt[i, j] = MOI.add_constraint(
                model,
                vars_flow[i, j] - U * vars_unicity[i, j],
                MOI.GreaterThan(zero(T))
            )
            cons_flow_lt[i, j] = MOI.add_constraint(
                model,
                vars_flow[i, j] - L * vars_unicity[i, j],
                MOI.LessThan(zero(T))
            )
        end
    end

    # Ensure that the end array is sorted.
    cons_sort = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}[
        MOI.add_constraint(
            model,
            f_sorted[i] - f_sorted[i + 1],
            MOI.GreaterThan(zero(T))
        )
        for i in 1:(dim - 1)
    ]

    return Sort2MILPBridge(
        vars_flow,
        vars_unicity,
        vars_unicity_bin,
        cons_transportation_x,
        cons_transportation_y,
        cons_unicity_x,
        cons_unicity_y,
        cons_flow_gt,
        cons_flow_lt,
        cons_sort,
    )
end

function MOI.supports_constraint(
    ::Type{Sort2MILPBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.Sort},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{Sort2MILPBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{Sort2MILPBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
    ]
end

function MOI.get(b::Sort2MILPBridge, ::MOI.NumberOfVariables)
    return length(b.vars_flow) + length(b.vars_unicity)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return length(b.vars_unicity)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return length(b.cons_flow_gt) + length(b.cons_sort)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return length(b.cons_flow_lt)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return length(b.cons_transportation_x) + length(b.cons_transportation_y) + length(b.cons_unicity_x) + length(b.cons_unicity_y)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return vcat(vec(b.vars_flow), vec(b.vars_unicity))
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return copy(b.vars_unicity_bin)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return vcat(vec(b.cons_flow_gt), b.cons_sort)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.LessThan{T},
    },
) where {T}
    return vec(b.cons_flow_lt)
end

function MOI.get(
    b::Sort2MILPBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.EqualTo{T},
    },
) where {T}
    return vcat(b.cons_transportation_x, b.cons_transportation_y, b.cons_unicity_x, b.cons_unicity_y)
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using PyCall

macro python_call(expr)
    quote
        try
            return $(esc(expr))
        catch e
            if isa(e, PyCall.PyError)
                printstyled("Uncaught Python exception:\n", bold = true, color = :red)
                traceback.print_exception(e.T, e.val, e.traceback)
            end
            rethrow()
        end
    end
end


text:
# граница xэмминга про упаковку шарами
function hammingND(n, d)
  t = div(d - 1, 2)
  binoms = [binomial(n, i) for i in 0:t]
  println("binoms = $binoms")
  z = sum(binoms)
  println("z = $z")
  M = 2^n / z
  return (M, log2(M))
end

function hammingNK(n, k)
  t = 0
  b = 2^(n - k)
  while sum([binomial(n, i) for i in 0:t]) ≤ b
    t += 1
  end
  return t - 1
end

function varshamovGuilbertND(n, d)
  k = 0
  binomSum = sum([binomial(n - 1, i) for i in 0:d-2])
  println("binomSum = $binomSum")

  while 2^(n - k) > binomSum
    k += 1
  end

  return k - 1
end

function varshamovGuilbertNK(n, k)
  d = 2
  b = 2^(n - k)
  while b > sum([binomial(n - 1, i) for i in 0:d-2])
    d += 1
  end
  return d - 1
end

function graismer(k, d)
  n = 0
  for i = 0:k-1
    n += ceil(Int64, d / 2^i)
  end
  return n
end

function graismerD(n, k)
  # n > sum(d/2^i) from 0 to k
  d = 0
  while true
    sums = sum([ceil(Int64, d / 2^i) for i in 0:k-1])
    if n > sums
      d += 1
    else
      break
    end
  end

  return d - 1
end

n, k, d = 21, 13, 14

# наxождение k по n и d
#
# M, klog = hammingND(n, d)
# println("граница Хэмминга: M ≤ $M, k ≈ $(klog)")

# println("==============================")

# mink = varshamovGuilbertND(n, d)
# println("граница Варшамова-Гилберта: k ≥ $mink")

# println("==============================")
# println(graismer(9, d))
# println("###################################")

# наxождение d по n и k
#
t = hammingNK(n, k)
println("граница Хэмминга: t ≤ $t")

println("==============================")

mind = varshamovGuilbertNK(n, k)
println("граница Варшамова-Гилберта: d ≥ $mind")

println("==============================")
println("граница Грайсмера: $(graismerD(n, k))")
println("###################################")

max_stars_repo_name: ChShersh/university-courses

max_issues_repo_name: ChShersh/university-courses

max_forks_repo_name: ChShersh/university-courses


text:
# Generates the files exp_generated/exp_X.jl using the GraphMatFun package
using GraphMatFun, LinearAlgebra

rhov=[0; 0.015; 0.25; 0.95; 2.1; 5.4];
for s=1:8
    push!(rhov, rhov[end]*2);
end

for i=1:(size(rhov,1)-1)
    r=(rhov[i]+rhov[i+1])/2;
    # Get the exp_native_jl = Higham 2008 version of scaling and squaring
    (graph,_)=graph_exp_native_jl(r)
    compress_graph!(graph);
    A=randn(3,3); A=r*A/opnorm(A,1)
    E1=eval_graph(graph,A)
    E2=exp(A);
    ee=norm(E1-E2)/norm(E1);
    println("Generating graph $i for nA ∈ ($(rhov[i]), $(rhov[i+1]))");
    println("    Check that the graph gives same solution: $ee");
    alloc_function=k-> "getmem(cache,$k)"; # Special type of mem allocation
    lang=LangJulia(true,true,true,false,alloc_function,true)

    fname="exp_$(i).jl";
    priohelp=Dict();
    if i==4 # The topological sorter needs some help. cf https://github.com/matrixfunctions/GraphMatFun.jl/issues/61
        priohelp[:V]=-100000;
    end
    gen_code(fname,graph,
             lang=lang,funname="exp_$(i)",precomputed_nodes=[:A],
             priohelp=priohelp);

    # Post-processing
    lines=[];
    open(fname, "r") do infile
        lines=readlines(infile);
    end

    open(fname,"w") do outfile
        for (j,line) in enumerate(lines)
            # Function name with Val-trick
            line=replace(line,"exp_$(i)!(A" => "exp_gen!(cache,A,::Val{$(i)}");
            if (contains(line,"ValueOne")) # Not needed functionality
                continue;
            end
            # Make a LAPACK call instead of backslash
            #line=replace(line, r"memslots(\d+)\s*.?=\s*memslots(\d+)..?memslots(\d+)" => s"LAPACK.gesv!(memslots\2, memslots\3); memslots\1=memslots\3")
            line=replace(line, r"memslots(\d+)\s*.?=\s*memslots(\d+)..?memslots(\d+)" => s"ldiv_for_generated!(memslots\1, memslots\2, memslots\3)")
            # Make sure the output is in A
            line=replace(line, r"return memslots(\d+)" => s"copyto!(A,memslots\1)")

            # Support for Julia 1.6.2 with manual inplace add for uniformscaling
            line=replace(line,
                         r"mul!\(memslots(\d+),true,I.coeff(\d+),true,true\)"
                         => s"inplace_add!(memslots\1,I*coeff\2)");

            println(outfile,line);
        end
    end

end

max_stars_repo_name: JuliaDiffEq/ExponentialUtilities

max_issues_repo_name: JuliaDiffEq/ExponentialUtilities

max_forks_repo_name: JuliaDiffEq/ExponentialUtilities


text:
using Oceananigans.Fields: compute_at!

import Oceananigans.OutputWriters:
                        fetch_output,
                        construct_output,
                        serializeproperty!

# This is working just fine at the moment?
# But it will be veeeeery slow, as reconstruct_global_field is not
# a performant operation
function fetch_output(mrf::MultiRegionField, model)
    field = reconstruct_global_field(mrf)
    compute_at!(field, time(model))
    return parent(field)
end

function construct_output(mrf::MultiRegionField, grid, user_indices, with_halos)
    grid = user_output.grid

    # TODO: support non-default indices I guess
    # for that we have to figure out how to partition indices, eg user_indices is "global"
    # indices = output_indices(user_output, grid, user_indices, with_halos)

    indices = (:, :, user_indices[3]) # sorry user

    return construct_output(mrf, indices)
end

function serializeproperty!(file, location, mrf::MultiRegionField{LX, LY, LZ}) where {LX, LY, LZ}
    p = reconstruct_global_field(mrf)
    serializeproperty!(file, location * "/location", (LX(), LY(), LZ()))
    serializeproperty!(file, location * "/data", parent(p))
    serializeproperty!(file, location * "/boundary_conditions", p.boundary_conditions)

    return nothing
end

function serializeproperty!(file, location, mrg::MultiRegionGrid)
    file[location] = on_architecture(CPU(), reconstruct_global_grid(mrg))
end


text:
@testset "Cycles" begin
    completedg = CompleteDiGraph(4)
    pathdg = PathDiGraph(5)
    triangle = random_regular_graph(3, 2)
    quadrangle = random_regular_graph(4, 2)
    pentagon = random_regular_graph(5, 2)

    @test maxsimplecycles(pathdg) == 0
    @test maxsimplecycles(completedg) == 20
    @test maxsimplecycles(4) == 20
    @test maxsimplecycles(pathdg, false) == 84

    @test length(simplecycles(completedg)) == 20
    @test simplecycles(completedg) == @inferred(simplecycles_iter(completedg))
    @test simplecyclescount(completedg) == 20
    @test simplecycleslength(completedg) == ([0, 6, 8, 6], 20)

    @test simplecyclescount(pathdg) == 0
    @test length(simplecycles(pathdg)) == 0
    @test isempty(simplecycles(pathdg)) == true
    @test isempty(simplecycles_iter(pathdg)) == true
    @test simplecycleslength(pathdg) == (zeros(5), 0)

    @test simplecyclescount(completedg, 10) == 10
    @test simplecycleslength(completedg, 10)[2] == 10

    @test simplecyclescount(pathdg, 10) == 0
    @test isempty(simplecycles_iter(pathdg, 10)) == true
    @test simplecycleslength(pathdg, 10) == (zeros(5), 0)

    trianglelengths, triangletotal = simplecycleslength(DiGraph(triangle))
    @test sum(trianglelengths) == triangletotal

    quadranglelengths, quadrangletotal = simplecycleslength(DiGraph(quadrangle))
    @test sum(quadranglelengths) == quadrangletotal
    @test simplecycles(DiGraph(quadrangle)) == @inferred(simplecycles_iter(DiGraph(quadrangle)))

    pentagonlengths, pentagontotal = simplecycleslength(DiGraph(pentagon))
    @test sum(pentagonlengths) == pentagontotal

end


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

#=
In Julia the if condition is simple, Unlike python, matlab it won't use vectorized code.

    `if condition.....end`

    `if condition....else.....end`

=#

println("***********************1st Example*************************")

a = rand()
println("lets check the value of a $a")
if a>0.5
    println("wow!! great")

end

println("***********************2nd Example*************************")

a = rand()
println("lets check the value of a $a")
if a>0.5
    println("wow!! great")
else
    println("normal")

end

println("***********************3rd Example*************************")
a = rand()
b = rand()

c=((a>b) ? println("a is greater") : ("b is greater"))
print(c)

println("***********************4th Example*************************")
a = rand()
b = rand()

c=((a>b) && println("a is greater"))
print(c)


text:
# ------------ Com-pirate fixes ------------ #

# In the future, I hope to remove some of these.

import Zygote.literal_getproperty

# Fix for: https://github.com/FluxML/Zygote.jl/issues/717
Zygote.@adjoint function literal_getproperty(x, ::Val{f}) where f
    val = getproperty(x, f)
    function back(Δ)
        Zygote.accum_param(__context__, val, Δ) # === nothing && return
        if isimmutable(x)
            ((;Zygote.nt_nothing(x)..., Zygote.pair(Val(f), Δ)...), nothing)
        else
            dx = Zygote.grad_mut(__context__, x)
            dx[] = (;dx[]...,Zygote.pair(Val(f), Zygote.accum(getfield(dx[], f), Δ))...)
            return (dx, nothing)
        end
    end
    unwrap(val), back
end

# Whitelist includes vectorized calls. (Required for performance of tracing).
whitelist = [
             # Jaynes custom indicator call.
             :trace,

             # Base.
             :map, :filter, :_apply_iterate, :collect,

             # Interactions with the context.
             :learnable, :fillable, :factor,
            ]

# Fix for specialized tracing. (Not a big deal).
function recur(ir, to = self)
    pr = Pipe(ir)
    for (x, st) in pr
        isexpr(st.expr, :call) && begin
            ref = unwrap(st.expr.args[1])
            ref in whitelist || continue
            pr[x] = Expr(:call, to, st.expr.args...)
        end
    end
    finish(pr)
end

# Fix for _apply_iterate (used in contexts). (More annoying).
function f_push!(arr::Array, t::Tuple{}) end
f_push!(arr::Array, t::Array) = append!(arr, t)
f_push!(arr::Array, t::Tuple) = append!(arr, t)
f_push!(arr, t) = push!(arr, t)
function flatten(t::Tuple)
    arr = Any[]
    for sub in t
        f_push!(arr, sub)
    end
    return arr
end


text:
module HonlyHelperBouncySpikes

using ..Ahorn, Maple

@mapdef Entity "HonlyHelper/BouncySpikesUp" BouncySpikesUp(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, texture::String="objects/HonlyHelper/BouncySpikes/bouncer", FreezeFrameEnable::Bool=false)
@mapdef Entity "HonlyHelper/BouncySpikesDown" BouncySpikesDown(x::Integer, y::Integer, width::Integer=Maple.defaultSpikeWidth, texture::String="objects/HonlyHelper/BouncySpikes/bouncer", FreezeFrameEnable::Bool=false)
@mapdef Entity "HonlyHelper/BouncySpikesLeft" BouncySpikesLeft(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, texture::String="objects/HonlyHelper/BouncySpikes/bouncer", FreezeFrameEnable::Bool=false)
@mapdef Entity "HonlyHelper/BouncySpikesRight" BouncySpikesRight(x::Integer, y::Integer, height::Integer=Maple.defaultSpikeHeight, texture::String="objects/HonlyHelper/BouncySpikes/bouncer", FreezeFrameEnable::Bool=false)

const placements = Ahorn.PlacementDict()

BouncySpikes = Dict{String, Type}(
    "up" => BouncySpikesUp,
    "down" => BouncySpikesDown,
    "left" => BouncySpikesLeft,
    "right" => BouncySpikesRight
)

BouncySpikesUnion = Union{BouncySpikesUp, BouncySpikesDown, BouncySpikesLeft, BouncySpikesRight}

        for (dir, entity) in BouncySpikes
            key = "Bouncy Spikes ($(uppercasefirst(dir))) (HonlyHelper)"
            placements[key] = Ahorn.EntityPlacement(
                entity,
                "rectangle",
            )
        end

directions = Dict{String, String}(
    "HonlyHelper/BouncySpikesUp" => "up",
    "HonlyHelper/BouncySpikesDown" => "down",
    "HonlyHelper/BouncySpikesLeft" => "left",
    "HonlyHelper/BouncySpikesRight" => "right",
)

offsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (4, -4),
    "down" => (4, 4),
    "left" => (-4, 4),
    "right" => (4, 4),
)

BouncySpikesOffsets = Dict{String, Tuple{Integer, Integer}}(
    "up" => (0, 0),
    "down" => (0, -0),
    "left" => (0, 0),
    "right" => (-0, 0),
)

rotations = Dict{String, Number}(
    "up" => 0,
    "right" => pi / 2,
    "down" => pi,
    "left" => pi * 3 / 2
)

resizeDirections = Dict{String, Tuple{Bool, Bool}}(
    "up" => (true, false),
    "down" => (true, false),
    "left" => (false, true),
    "right" => (false, true),
)

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::BouncySpikesUnion)
    direction = get(directions, entity.name, "up")
    theta = rotations[direction] - pi / 2

    width = Int(get(entity.data, "width", 0))
    height = Int(get(entity.data, "height", 0))

    x, y = Ahorn.position(entity)
    cx, cy = x + floor(Int, width / 2) - 8 * (direction == "left"), y + floor(Int, height / 2) - 8 * (direction == "up")

    Ahorn.drawArrow(ctx, cx, cy, cx + cos(theta) * 24, cy + sin(theta) * 24, Ahorn.colors.selection_selected_fc, headLength=6)
end

function Ahorn.selection(entity::BouncySpikesUnion)
    if haskey(directions, entity.name)
        x, y = Ahorn.position(entity)

        width = Int(get(entity.data, "width", 8))
        height = Int(get(entity.data, "height", 8))

        direction = get(directions, entity.name, "up")

        ox, oy = offsets[direction]

        return Ahorn.Rectangle(x + ox - 4, y + oy - 4, width, height)
    end
end

Ahorn.minimumSize(entity::BouncySpikesUnion) = 8, 8

function Ahorn.resizable(entity::BouncySpikesUnion)
    if haskey(directions, entity.name)
        direction = get(directions, entity.name, "up")

        return resizeDirections[direction]
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::BouncySpikesUnion)
    if haskey(directions, entity.name)
        direction = get(directions, entity.name, "up")
	texture_used = get(entity.data, "texture", "objects/HonlyHelper/BouncySpikes/bouncer")
        BouncySpikesOffset = BouncySpikesOffsets[direction]

        width = get(entity.data, "width", 8)
        height = get(entity.data, "height", 8)

        for ox in 0:8:width - 8, oy in 0:8:height - 8
            drawX, drawY = (ox, oy) .+ offsets[direction] .+ BouncySpikesOffset
            Ahorn.drawSprite(ctx, "$(texture_used)_$(direction)00", drawX, drawY)
        end
    end
end

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
using UnitfulOfu, Test

# write your own tests here

println("Running tests for UnitfulOfu:")

tests = [
  "ofu_test.jl"
]

for test in tests
    println("\n  * $(test) *\n")
    include(test)
    println()
end


text:
using BinningAnalysis
using Test, Statistics, Random

@testset "All Tests" begin

    @testset "Logarithmic Binning" begin
        include("logbinning.jl")
    end

    @testset "Full Binning" begin
        include("fullbinning.jl")
    end

    @testset "Generic functions" begin
        x = [0.0561823, 0.846613, 0.813439, 0.357134, 0.157445, 0.103298, 0.948842, 0.629425, 0.290206, 0.00695332, 0.869828, 0.949165, 0.897995, 0.916239, 0.457564, 0.349827, 0.398683, 0.264218, 0.72754, 0.934315, 0.666448, 0.134813, 0.364933, 0.829088, 0.256443, 0.595029, 0.172097, 0.241686, 0.489935, 0.239663, 0.391291, 0.00751015, 0.138935, 0.0569876, 0.571786, 0.694996, 0.798602, 0.923308, 0.73978, 0.414774, 0.835145, 0.731303, 0.271647, 0.707796, 0.00348624, 0.0905812, 0.316176, 0.921054, 0.131037, 0.599667, 0.805071, 0.440813, 0.086516, 0.363658, 0.476161, 0.931257, 0.28974, 0.78717, 0.60822, 0.144024, 0.214432, 0.061922, 0.626495, 0.512072, 0.758078, 0.840485, 0.242576, 0.147441, 0.599222, 0.993569, 0.0365044, 0.0983033, 0.713144, 0.422394, 0.480044, 0.968745, 0.518475, 0.431319, 0.4432, 0.526007, 0.612975, 0.468387, 0.262145, 0.888011, 0.105744, 0.325821, 0.769525, 0.289073, 0.336083, 0.443037, 0.489698, 0.141654, 0.915284, 0.319068, 0.341001, 0.704346, 0.0794996, 0.0412352, 0.70016, 0.0195158]
        @test std_error(x) ≈ std_error(x; method=:log)
        @test std_error(x; method=:log) ≈ 0.03384198827323159
        @test std_error(x; method=:full) ≈ 0.030761063167290947
        @test std_error(x; method=:jackknife) ≈ 0.029551526503551463
        @test tau(x) ≈ tau(x; method=:log)
        @test tau(x; method=:log) ≈ 0.15572524869061377
        @test tau(x; method=:full) ≈ 0.04176737474771308

        x = x .+ 1im
        @test std_error(x)≈ std_error(x; method=:log)
        @test std_error(x; method=:log) ≈ 0.0338419882732316
        @test std_error(x; method=:full) ≈ 0.03076106316729094
        @test std_error(x; method=:jackknife) ≈ 0.029551526503551445
        @test tau(x) ≈ tau(x; method=:log)
        @test tau(x; method=:log) ≈ 0.15572524869061377
        @test tau(x; method=:full) ≈ 0.04176737474771308

        x = Array{Float64,2}[[0.311954 0.476706 0.314101; 0.420978 0.478085 0.0194284], [0.491386 0.493583 0.295477; 0.136896 0.1634 0.84641], [0.318536 0.705903 0.121377; 0.764174 0.240484 0.25894], [0.286923 0.172458 0.392881; 0.124348 0.140628 0.730131], [0.287451 0.221914 0.382938; 0.29568 0.249575 0.87685]]
        @test std_error(x) ≈ std_error(x; method=:log)
        @test std_error(x, method=:log) ≈ [0.03856212859659066 0.09765048982846934 0.04879645809318544; 0.11744109653881814 0.05978167652935139 0.1722342286233907]
        @test std_error(x, method=:full) ≈ [0.03856212859659072 0.09765048982846936 0.048796458093185405; 0.11744109653881814 0.05978167652935137 0.17223422862339066]
        @test tau(x) ≈ tau(x; method=:log)
        @test tau(x; method=:log) ≈ [0.0 0.0 0.0; 0.0 0.0 0.0]
        @test isapprox(tau(x; method=:full), [0.0 -5.55112e-17 -5.55112e-17; -1.11022e-16 0.0 0.0], atol=1e-16)

        x = ["this", "should", "error"]
        @test_throws ErrorException std_error(x, method=:jackknife)
        @test_throws ArgumentError std_error(x, method=:whatever)
    end

    @testset "Jackknife" begin
        g(x1, x2) = x1^2 - x2

        # Real
        ts = [0.00124803, 0.643089, 0.183268, 0.799899, 0.0857666, 0.955348, 0.165763, 0.765998, 0.63942, 0.308818]
        ts2 = [0.606857, 0.0227746, 0.805997, 0.978731, 0.0853112, 0.311463, 0.628918, 0.0190664, 0.515998, 0.0223728]
        @test isapprox.(
            jackknife_full(identity, ts),
            (0.45486176300000025, 0.0, 0.10834619757501414),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(identity, 1 ./ ts),
            (83.4370988407286, 0.0, 79.76537738034833),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(g, ts, ts2.^2),
            (-0.07916794438503649, 0.011738898528964592, 0.14501699232741938),
            atol = 1e-12
        ) |> all

        # check consistency with Julia's Statistics
        @test jackknife(identity, ts)[2] ≈ std(ts)/sqrt(length(ts))

        # Complex
        ts = Complex{Float64}[0.0259924+0.674798im, 0.329853+0.558688im, 0.821612+0.142805im, 0.0501703+0.801068im, 0.0309707+0.877745im, 0.937856+0.852463im, 0.669084+0.606286im, 0.887004+0.431615im, 0.763452+0.210563im, 0.678384+0.428294im]
        ts2 = Complex{Float64}[0.138147+0.11007im, 0.484956+0.127761im, 0.986078+0.702827im, 0.45161+0.878256im, 0.768398+0.954537im, 0.228518+0.260732im, 0.256892+0.437918im, 0.647672+0.749172im, 0.0587658+0.408715im, 0.0792651+0.288578im]
        @test isapprox.(
            jackknife_full(identity, ts),
            (
                0.519437840000001 + 0.5584324999999994im,
                -9.992007221626409e-16 + 0.0im,
                0.1428607353199917
            ),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(identity, 1 ./ ts),
            (
                0.6727428408881817 - 0.8112743686359876im,
                -9.992007221626409e-16 + 1.9984014443252818e-15im,
                0.20473472232587642
            ),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(g, ts, ts2.^2),
            (
                0.021446672721215643 + 0.06979705636190503im,
                0.007300401672734998 - 0.01055311354806504im,
                0.2752976586383889
            ),
            atol = 1e-12
        ) |> all

        # check consistency with Julia's Statistics
        @test jackknife(identity, ts)[2] ≈ std(ts)/sqrt(length(ts))
    end
end


text: println(STDERR, "Goodbye, World!")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#
#
# TODO: Fourier spectral, SpectralElements will require different functionality
#   - move broadcasting/indexing to AbstractField,
#   - if Fourier can be handled with vanilla Field, then rename it
#     SpectralField, else rename to LagrangePolyField
#
# TODO: create SpectralElementField and overload inner product by
#   - overload *(Adjoint{Field}, Field), norm(::Field, 2)
#
""" Scalar function field in D-dimensional space over a spectral basis"""
struct Field{T,D,Tarr <: AbstractArray{T,D}} <: AbstractField{T,D}
    array::Tarr
end

# display
function Base.summary(io::IO, u::Field{T,D}) where{T,D}
    println(io, "$(D)D scalar field of type $T")
    Base.show(io, typeof(u))
end
function Base.show(io::IO, ::MIME"text/plain", u::Field)
    iocontext = IOContext(io, :compact => true, :limit => true)
    Base.summary(iocontext, u)
    Base.show(iocontext, MIME"text/plain"(), u.array)
    println()
end

# vector indexing
Base.IndexStyle(::Field) = IndexLinear()
Base.getindex(u::Field, i::Int) = getindex(u.array, i)
Base.setindex!(u::Field, v, i::Int) = setindex!(u.array, v, i)
Base.size(u::Field) = (length(u.array),)

# allocation
Base.similar(u::Field, ::Type{T} = eltype(u), dims::Dims = size(u.array)) where{T} = Field(similar(u.array, T, dims))
Base.zero(u::Field, dims::Dims) = zero(u) # ignore dims since <: AbstractVector

# broadcast
Base.Broadcast.BroadcastStyle(::Type{<:Field}) = Broadcast.ArrayStyle{Field}()
function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{Field}},
                      ::Type{ElType}) where ElType
  u = find_fld(bc)
  Field(similar(Array{ElType}, axes(u.array)))
end

find_fld(bc::Base.Broadcast.Broadcasted) = find_fld(bc.args)
find_fld(args::Tuple) = find_fld(find_fld(args[1]), Base.tail(args))
find_fld(x) = x
find_fld(::Tuple{}) = nothing
find_fld(a::Field, rest) = a
find_fld(::Any, rest) = find_fld(rest)
#


text:
module Escape

using Base.Threads
using LinearAlgebra, Statistics, DelimitedFiles, Distributed, Serialization
using CSV, StanInterface, Graphs, MetaGraphs, Suppressor, DataFrames,
      HypothesisTests, MultipleTesting, BioSequences, FASTX,
      StatsBase, StatsFuns, Distributions, Loo, RecipesBase, Base.Threads, Conda,
      StaticArrays, ProgressMeter, JLD2, HAMdetector_model_binaries_jll, raxml_ng_jll

include("alleles.jl")
include("abstract_hla_data.jl")
include("replacement.jl")
include("abstract_hla_model.jl")
include("hla_dataset.jl")
include("dataset_rousseau.jl")
include("phylogenetic_tree.jl")
include("hla_data.jl")
include("split_hla_data.jl")
include("phylogenetic_background.jl")
include("epitope_prediction.jl")
include("epitope_map.jl")
include("hla_model.jl")
include("fisher_test.jl")
include("hla_model_result.jl")
include("loo.jl")
include("posterior_predictive_checks.jl")
include("split_hla_result.jl")

end


text:
#=
    NEMO: Next-generation Energy Modeling system for Optimization.
    https://github.com/sei-international/NemoMod.jl

    Copyright © 2018: Stockholm Environment Institute U.S.

    File description: Utilities for ahead-of-time compilation of |nemo. Typical usage in a REPL session:
        using NemoMod
        include(normpath(joinpath(pathof(NemoMod), "..", "..", "utils", "compilation.jl")))
        compilenemo()
=#

# BEGIN: Ensure all packages needed for stock |nemo compilation are installed.
using Pkg

if length(setdiff(["PackageCompiler", "Libdl", "JuMP", "SQLite", "DataFrames", "Distributed", "Dates", "ConfParser", "GLPKMathProgInterface", "Cbc"],
    collect(keys(Pkg.installed())))) > 0

    Pkg.add(setdiff(["PackageCompiler", "Libdl", "JuMP", "SQLite", "DataFrames", "Distributed", "Dates", "ConfParser", "GLPKMathProgInterface", "Cbc"],
        collect(keys(Pkg.installed()))))
end
# END: Ensure all packages needed for stock |nemo compilation are installed.

using PackageCompiler, Libdl, NemoMod

"""Clone of PackageCompiler.snoop_userimg(). Implemented here in order to allow a custom blacklist to be passed to PackageCompiler.snoop().
    This should be a temporary measure eventually obviated by upgrades to PackageCompiler (i.e., for compatability with Julia 1.0)."""
function snoop_userimg_nemo(userimg, packages::Tuple{String, String}...)
    snooped_precompiles = map(packages) do package_snoopfile
        package, snoopfile = package_snoopfile
        abs_package_path = if ispath(package)
            normpath(abspath(package))
        else
            normpath(Base.find_package(package), "..", "..")
        end
        file2snoop = normpath(abspath(joinpath(abs_package_path, snoopfile)))
        package = PackageCompiler.package_folder(PackageCompiler.get_root_dir(abs_package_path))
        isdir(package) || mkpath(package)
        precompile_file = joinpath(package, "precompile.jl")
        PackageCompiler.snoop(file2snoop, precompile_file, joinpath(package, "snooped.csv"),
            Vector{Pair{String, String}}(), String["Main", ", _}"])
        return precompile_file
    end
    # merge all of the temporary files into a single output
    open(userimg, "w") do output
        println(output, """
            # Prevent this from being put into the Main namespace
            Core.eval(Module(), quote
            """)
        for (pkg, _) in packages
            println(output, """
                import $pkg
                """)
        end
        println(output, """
            for m in Base.loaded_modules_array()
                Core.isdefined(@__MODULE__, nameof(m)) || Core.eval(@__MODULE__, Expr(:(=), nameof(m), m))
            end
            """)
        for path in snooped_precompiles
            open(input -> write(output, input), path)
            println(output)
        end
        println(output, """
            end) # eval
            """)
    end
    nothing
end  # snoop_userimg_nemo(userimg, packages::Tuple{String, String}...)

"""Replaces the Julia system image file with the image whose full path is specified in src.
    The logic in this function is necessary to avoid Julia locks on the system image file."""
function copysysimage(src::String)
    local sysimagepath::String = joinpath(PackageCompiler.default_sysimg_path(false), "sys.$(Libdl.dlext)")
        # Path to Julia system image file
    local sysimagedir::String = dirname(sysimagepath)
        # Directory containing Julia system image file

    # BEGIN: Try to delete any existing sys.dll.nemo_compile_backup files.
    for f in readdir(sysimagedir)
        try
            startswith(f, "sys.dll.nemo_compile_backup") && rm(joinpath(sysimagedir, f))
        catch
            # Just continue
        end
    end
    # END: Try to delete any existing sys.dll.nemo_compile_backup files.

    # BEGIN: Rename Julia system image file.
    if isfile(sysimagepath)
        local sysimagemoved::Bool = false  # Indicates whether Julia system image file has been renamed from "sys.dll"
        local counter::Int = 0  # Successive integer to append to new name for Julia system image file

        while !sysimagemoved
            counter += 1

            try
                mv(sysimagepath, sysimagepath * ".nemo_compile_backup" * string(counter); force = true)
                sysimagemoved = true
            catch
                # Just continue
            end
        end
    end
    # END: Rename Julia system image file.

    # Copy src to system image
    cp(src, sysimagepath)
end  # copysysimage(src::String)

"""Generates a new Julia system image including specified additional packages. Arguments:
    • packagesandusecases - One or more tuples of the following form:
        > First item - Name of the additional package to be included in the new system image.
        > Second item - Use cases that determine what functionality of the package is included in the new system image.
Example invocations:
    • compilenemo() - New system image includes NemoMod with compiled functionality based on standard NemoMod test cases. Support for GLPK and Cbc solvers is included
        in new system image since they're |nemo's default solvers.
    • compilenemo(("NemoMod", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "runtests.jl"))),
        ("CPLEX", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "test_cplex.jl")))) - New system image includes NemoMod and CPLEX; compiled functionality
        is based on standard NemoMod test cases plus a supplemental test case for CPLEX. Support for GLPK, Cbc, and CPLEX solvers is included in new system image."""
function compilenemo(packagesandusecases::Tuple{String, String}...=
    ("NemoMod", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "runtests.jl"))))
    # BEGIN: Restore stock system image.
    #= Background:
    |nemo all-in-one installer backs up original system image as lib\julia\sys.dll.nemo_backup.
    Calling PackageCompiler.compile_package() with force = true pushes current system image to sys.dll.packagecompiler_backup in lib\julia
        (if this file doesn't exist already).
    PackageCompiler.revert() copies any sys.dll from CPU target-specific subdirectory under PackageCompiler sysimg folder to lib\julia.
        Such sys.dll's are clean system images (no extra packages added) created by PackageCompiler (revert creates one).

    Logic:
        1) Restore lib\julia\sys.dll.nemo_backup if available.
        2) Otherwise, restore lib\julia\sys.dll.packagecompiler_backup if available.
        3) Otherwise, assume sys.dll is stock. Copy sys.dll to lib\julia\sys.dll.nemo_backup for future reference.
    =#
    local sysimagedir::String = PackageCompiler.default_sysimg_path(false)  # Directory containing Julia system image file

    if isfile(joinpath(sysimagedir, "sys.$(Libdl.dlext).nemo_backup"))
        copysysimage(joinpath(sysimagedir, "sys.$(Libdl.dlext).nemo_backup"))
    elseif isfile(joinpath(sysimagedir, "sys.$(Libdl.dlext).packagecompiler_backup"))
        copysysimage(joinpath(sysimagedir, "sys.$(Libdl.dlext).packagecompiler_backup"))
    else
        cp(joinpath(sysimagedir, "sys.$(Libdl.dlext)"), joinpath(sysimagedir, "sys.$(Libdl.dlext).nemo_backup"))
    end
    # END: Restore stock system image.

    # Create snooping file
    snoop_userimg_nemo(PackageCompiler.sysimg_folder("precompile.jl"), packagesandusecases...)

    # Compile new system image
    local imgfile::String = PackageCompiler.compile_package("NemoMod"; reuse = true, force = false, cpu_target = "native")

    # Replace system image
    copysysimage(imgfile)
    @info "compilenemo() replaced system image. This supersedes prior messages from PackageCompiler."
end  # compilenemo(packagesandusecases::Tuple{String, String}...= ("NemoMod", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "runtests.jl"))))

"""A convenience function for compiling |nemo with support for Gurobi, CPLEX, and/or Mosek
    included (corresponding parameter = true) or excluded (corresponding parameter = false).
    Support for GLPK and Cbc is included in all cases, even if all parameters are false."""
function compilenemo(; gurobi = false, cplex = false, mosek = false)
    # BEGIN: Ensure Julia packages for Gurobi, CPLEX, and Mosek, as needed, are installed.
    local pkgs::Array{String,1} = Array{String,1}()
        # Array of names of packages whose installation will be verified

    gurobi && push!(pkgs, "Gurobi")
    cplex && push!(pkgs, "CPLEX")
    mosek && push!(pkgs, "Mosek")

    pkgs = setdiff(pkgs, collect(keys(Pkg.installed())))

    length(pkgs) > 0 && Pkg.add(pkgs)
    # END: Ensure Julia packages for Gurobi, CPLEX, and Mosek, as needed, are installed.

    # BEGIN: Call compilenemo() with appropriate packagesandusecases.
    local puc::Array{Tuple{String,String},1} = Array{Tuple{String,String},1}()
        # Array of Tuples for packagesandusecases

    push!(puc, ("NemoMod", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "runtests.jl"))))
    gurobi && push!(puc, ("Gurobi", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "test_gurobi.jl"))))
    cplex && push!(puc, ("CPLEX", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "test_cplex.jl"))))
    mosek && push!(puc, ("Mosek", normpath(joinpath(pathof(NemoMod), "..", "..", "test", "test_mosek.jl"))))

    compilenemo(puc...)
    # END: Call compilenemo() with appropriate packagesandusecases.
end  # compilenemo(; gurobi = false, cplex = false, mosek = false)


text:
function foo(ch; xi=0, dx=1, xf=1)
    for _x in xi : dx : xf
        put!(ch, _x)
    end
end

ch = Channel(0)
task = @async foo(ch, xf = 10)
bind(ch, task)

for i in 1 : 9
    @show take!(ch)
end
task


text:
##############################################################################
##
## Parse FixedEffect
##
##
##############################################################################

function parse_fixedeffect(df::AbstractDataFrame, feformula)
    fe = FixedEffect[]
    id = Symbol[]
    for term in feformula.terms
        result = parse_fixedeffect(df, term)
        if result != nothing
            push!(fe, result[1])
            push!(id, result[2])
        end
    end
    return fe, id
end

# Constructor
parse_fixedeffect(x::Nothing) = nothing

# Constructors from dataframe + symbol
function parse_fixedeffect(df::AbstractDataFrame, a::Symbol)
    v = df[a]
    if isa(v, CategoricalVector)
        # x from x*id -> x + id + x&id
        return FixedEffect(v), a
    end
end

# Constructors from dataframe + expression
function parse_fixedeffect(df::AbstractDataFrame, a::Expr)
    _check(a) || throw("Expression $a should only contain & and variable names")
    factorvars, interactionvars = _split(df, allvars(a))
    if !isempty(factorvars)
        # x1&x2 from (x1&x2)*id
        fe = FixedEffect((df[v] for v in factorvars)...; interaction = _multiply(df, interactionvars))
        id = _name(allvars(a))
        return fe, id
    end
end

function _check(a::Expr)
    a.args[1] == :& && check(a.args[2]) && check(a.args[3])
end
check(a::Symbol) = true

function _split(df::AbstractDataFrame, ss::Vector{Symbol})
    catvars, contvars = Symbol[], Symbol[]
    for s in ss
        isa(df[s], CategoricalVector) ? push!(catvars, s) : push!(contvars, s)
    end
    return catvars, contvars
end

function _multiply(df, ss::Vector{Symbol})
    if isempty(ss)
        out = Ones(size(df, 1))
    else
        out = ones(size(df, 1))
        for j in 1:length(ss)
            _multiply!(out, df[ss[j]])
        end
    end
    return out
end

function _multiply!(out, v)
    for i in 1:length(out)
        if v[i] === missing
            # may be missing when I remove singletons
            out[i] = 0.0
        else
            out[i] = out[i] * v[i]
        end
    end
end

function _name(s::Vector{Symbol})
    if isempty(s)
        out = nothing
    else
        out = Symbol(reduce((x1, x2) -> string(x1)*"x"*string(x2), s))
    end
    return out
end


text:
# ## Generating dummy data
#
# Let's start by generating some dummy data with both numerical values and categorical values:

using MLJ, PrettyPrinting
MLJ.color_off() # hide
@load KNNRegressor
## input
X = (age    = [23, 45, 34, 25, 67],
     gender = categorical(['m', 'm', 'f', 'm', 'f']))
## target
height = [178, 194, 165, 173, 168];

# Note that the scientific type of `age` is `Count` here:

scitype(X.age)

# We will want to coerce that to `Continuous` so that it can be given to a regressor that expects such values.
#
# ## Declaring a pipeline
#
# A typical workflow for such data is to one-hot-encode the categorical data and then apply some regression model on the data.
# Let's say that we want to apply the following steps:
# 1. standardize the target variable (`:height`)
# 1. one hot encode the categorical data
# 1. train a KNN regression model
#
# The `@pipeline` macro helps you define such a simple (non-branching) pipeline of steps to be applied in order:

pipe = @pipeline MyPipe(X -> coerce(X, :age=>Continuous),
                       hot = OneHotEncoder(),
                       knn = KNNRegressor(K=3),
                       target = UnivariateStandardizer());

# Note the coercion of the `:age` variable to Continuous since `KNNRegressor` expects `Continuous` input.
# Note also the `target` keyword where you can specify a transformation of the target variable.
#
# Hyperparameters of this pipeline can be accessed (and set) using dot syntax:

pipe.knn.K = 2
pipe.hot.drop_last = true;

# Evaluation for a pipe can be done with the `evaluate!` method; implicitly it will construct machines that will contain the fitted parameters etc:

evaluate(pipe, X, height, resampling=Holdout(),
         measure=rms) |> pprint

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials