text:
AbstractAn

"""
    Temporal Superior Longitudinal Fasciculus
"""
const TemporalSuperiorLongitudinalFasciculus = BilateralStructure("temporal superior longitudinal fasciculus")

"""
    ParietalSuperiorLongitudinalFasciculus
"""
const ParietalSuperiorLongitudinalFasciculus = BilateralStructure("parietal superior longitudinal fasiculus")

"""
Fornix (F)	(Catani et al., 2002; Catani and Thiebaut de Schotten, 2008; Thiebaut de Schotten et al., 2011b; Wakana et al., 2004)
"""
const Fornix = BilateralStructure("fornix")

const ForcepsMajor = AnatomicalStructure("forceps major")
const PosteriorForceps = ForcepsMajor

const ForcepsMinor = AnatomicalStructure("forceps minor")
const AnteriorForceps = ForcepsMinor

const AnteriorThalamicRadiation = BilateralStructure("anterior thalamic radiation")

const ExternalCapsule = AnatomicalStructure("external capsule")

const InternalCapsule = AnatomicalStructure("external capsule")

# corpus callosum

###
AnteriorWhiteCommissure = AnatomicalStructure("anterior white commissure")
VentralWhiteCommissure = AnteriorWhiteCommissure

### Projection Pathways

"""
    AcousticRadiation (AR)
"""
const AcousticRadiation = AnatomicalStructure("acoustic radiation")

"""
    CorticospinalTract (CST)
"""
const CorticospinalTract = BilateralStructure("corticospinal tract")

"""
Corticobulbar Tract (CBT)	(Thiebaut de Schotten et al., 2011b)
"""
const CorticobulbarTract = AnatomicalStructure("corticobulbar Tract")

"""
Corticostriatal Pathway (CS)	(Fernandez-Miranda et al., 2008; Torgerson et al., 2015)
"""
const CorticostriatalPathway = AnatomicalStructure("corticostriatal pathway")

# TODO striatal inferior frontal cortex pathway
"""
    StriatalInferiorFrontalPathway
"""
const StriatalInferiorFrontalPathway = BilateralStructure("striatal inferior frontal cortex pathway")

# TODO superior cortic-striato-frontal pathway
"""
    SuperiorCorticoStriatoFrontalPathway
"""
const SuperiorCorticoStriatoFrontalPathway = BilateralStructure("superior cortic-striato-frontal pathway")

# TODO superior corticostriate-parietal cortex only
"""
    SuperiorCorticoStriatoParietalPathway
"""
const SuperiorCorticoStriatoParietalPathway = BilateralStructure("superior corticostriate-parietal cortex")

"""
    SuperiorCorticostriatePathway
"""
const SuperiorCorticostriatePathway = BilateralStructure(
    name = "superior corticostriate pathway",
    children = (
        SuperiorCorticoStriatoFrontalPathway,
    )
)

"""
    Corticothalamic Pathway (CT)	(Behrens et al., 2003; Wakana et al., 2004)
"""
const CorticothalamicPathway = BilateralStructure("corticothalamic pathway")

const SuperiorOpticRadiation = BilateralStructure("superior optic radiation")

const InferiorOpticRadiation = BilateralStructure("inferior optic radiation")

"""
Optic Radiation (OR)	(Rubino et al., 2005; Thiebaut de Schotten et al., 2011b)
"""
const OpticRadiation = BilateralStructure(
    name = "optic radiation",
    children = (
        SuperiorOpticRadiation,
        InferiorOpticRadiation
     )
   )
const GeniculocalcarineTract = OpticRadiation
const GeniculostriatePathway = OpticRadiation
const VisualRadiation = OpticRadiation

"""
Frontopontine Tract(FPT)	(Meola et al., 2016c)
"""
const FrontopontineTract = BilateralStructure("frontopontine tract")
const FrontopontineFibers = FrontopontineTract

"""
Occipitopontine Tract (OPT)	(Meola et al., 2016c)
"""
const OccipitopontineTract = BilateralStructure("occipitopontine tract")

"""
Parietopontine Tract (PPT)	(Meola et al., 2016c)
"""
const ParietopontineTract = BilateralStructure("parietopontine tract")

"""
Temporopontine Tract (TPT)	(Meola et al., 2016c)
"""
const TemporopontineTract = BilateralStructure("temporopontine tract")

### Association pathways
"""
Arcuate Fasciculus (AF)	(Axer et al., 2013; Fernandez-Miranda et al., 2015; Rilling et al., 2008; Saur et al., 2008; Thiebaut de Schotten et al., 2012)
"""
const ArcuateFasciculus = BilateralStructure("arcuate fasciculus")

"""
Cingulum (C)	(Catani et al., 2002; Catani and Thiebaut de Schotten, 2008; Thiebaut de Schotten et al., 2012; Thiebaut de Schotten et al., 2011b; Wakana et al., 2004)
"""
const Cingulum = BilateralStructure("cingulum")

"""
    ExtremeCapsule (EMC)

Extreme Capsule (EMC)	(Makris and Pandya, 2009)
"""
const ExtremeCapsule = BilateralStructure("extreme capsule")

# TODO Frontal Aslant Tract?
"""
Frontal Aslant Tract (FAT)

## Citations
* Catani et al., 2012
* Thiebaut de Schotten et al., 2012
"""
const FrontalAslantTract = AnatomicalStructure("frontal aslant tract")

"""
Inferior Fronto Occipital Fasciculus (IFOF)

* Catani et al., 2002
* Forkel et al., 2014
* Makris et al., 2007
* Meola et al., 2015
* Schmahmann and Pandya, 2007
* Thiebaut de Schotten et al., 2012
* Thiebaut de Schotten et al., 2011b
* Wakana et al., 2004)
"""
const InferiorFrontalOccipitalFasciculus = BilateralStructure("inferior frontal occipital fasciculus")

"""
InferiorLongitudinalFasciculus (ILF)	(Amemiya and Naito, 2016; Catani et al., 2002; Catani and Thiebaut de Schotten, 2008; Grossi et al., 2014; Koldewyn et al., 2014; Latini, 2015; Mandonnet et al., 2007; Thiebaut de Schotten et al., 2011b; Tusa and Ungerleider, 1985; Wakana et al., 2004)
"""
InferiorLongitudinalFasciculus = BilateralStructure("inferior longitudinal fasciculus")

"""
Middle Longitudinal Fasciculus (MdLF)

Makris et al., 2009
Makris et al., 2013
Wang et al., 2013
"""
const MiddleLongitudinalFasciculus = BilateralStructure("middle longitudinal fasciculus")

"""
Superior Longitudinal Fasciculus (SLF)	(Budisavljevic et al., 2016; Karlsgodt et al., 2008; Makris et al., 2005; Thiebaut de Schotten et al., 2011a; Wakana et al., 2004; Wang et al., 2016)
"""
const SuperiorLongitudinalFasciculus = BilateralStructure("superior longitudinal fasciculus")

# TODO U-fibers are pretty ambiguous for an ROI
"""

U-fiber (U)
"""
const Ufiber = AnatomicalStructure("U-fiber")

"""
Uncinate Fasciculus (UF)

* Catani and Thiebaut de Schotten, 2008
* Leng et al., 2016
* Thiebaut de Schotten et al., 2012
* Thiebaut de Schotten et al., 2011b
* Wakana et al., 2004
"""
const Uncinate = BilateralStructure("uncinate fasciculus")
const UncinateFasciculus = Uncinate

"""
Vertical Occipital Fasciculus (VOF)

* Takemura et al., 2016
* Yeatman et al., 2014
"""
const VerticalOccipitalFasciculus = BilateralStructure("vertical occipital fasciculus")

### Commisures

"""
Anterior Commissure (AC)

Catani et al., 2002
Catani and Thiebaut de Schotten, 2008
Thiebaut de Schotten et al., 2011b
"""
const AnteriorCommissure = AnatomicalStructure("anterior commissure")

"""
    Rostrum of Coprus Callosum
"""
const RostrumOfCorpusCallosum = AnatomicalStructure("rostrum of corpus callosum")

const GenuOfCorpusCallosum = AnatomicalStructure("genu of corpus callosum")

const BodyOfCorpusCallosum = AnatomicalStructure("body of corpus callosum")
const TrunkOfCorpusCallosum = BodyOfCorpusCallosum

const SpleniumOfCorpusCallosum = AnatomicalStructure("splenium of corpus callosum")

const IsthmusOfCorpusCallosum = AnatomicalStructure("isthmus of corpus callosum")

"""
Corpus Callosum (CC)

(Catani and Thiebaut de Schotten, 2008; Huang et al., 2005; Thiebaut de Schotten et al., 2011b; Wakana et al., 2004)
"""

const CorpusCallosum = AnatomicalStructure(
    name = "corpus callosum",
    children = (
        SpleniumOfCorpusCallosum,
        IsthmusOfCorpusCallosum,
        GenuOfCorpusCallosum,
        RostrumOfCorpusCallosum
    )
)

"""
Posterior Commissure (PC)	(Keene, 1938)
"""
const PosteriorCommissure = AnatomicalStructure("posterior commissure")

### Cerebellum
const SuperiorCerebellarPeduncle = AnatomicalStructure("superior cerebellar peduncle")
const MiddleCerebellarPeduncle = AnatomicalStructure("middle cerebellar peduncle")
const InferiorCerebellarPeduncle = AnatomicalStructure("inferior cerebellar peduncle")

### Brain stem
"""
    Central Tegmental Tract (CTT)
"""
const CentralTegmentalTract = BilateralStructure("central tegmental tract")

const DorsalLongitudinalFasciculus = BilateralStructure("dorsal longitudinal fasciculus")

const LateralLemniscus = BilateralStructure("lateral lemniscus")

const MedialLemniscus = BilateralStructure("medial lemniscus")

const MedialLongitudinalFasciculus = BilateralStructure("medial longitudinal fasciculus")

const RubrospinalTract = BilateralStructure("rubrospinal tract")

"""
Spinothalamic tract
Meola et al., 2016c
"""
const SpinothalamicTract = BilateralStructure("spinothalamic tract")


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "jemalloc"
version = v"5.3.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/jemalloc/jemalloc.git", "54eaed1d8b56b1aa528be3bdd1877e59c56fa90c")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/jemalloc/
autoconf

FLAGS=(--disable-initial-exec-tls)
if [[ "${target}" == *-freebsd* ]]; then
     FLAGS+=(--with-jemalloc-prefix)
fi
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} "${FLAGS[@]}"

make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libjemalloc", "jemalloc"], :libjemalloc)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

ext: jl


text:
VERSION < v"0.7.0-rc1" && __precompile__()
module RangeArrays

using Compat

include("matrix.jl")
include("repeatedrange.jl")

export RangeMatrix, RepeatedRangeMatrix

end # module


text:
@testset "book" begin

  @testset "level 1" begin
    res = book("BTC-USD", 1)
    @test isa(res, DataFrame)
    @test nrow(res) == 1
  end

  @testset "level 2" begin
    res = book("BTC-USD", 2)
    @test isa(res, DataFrame)
    @test nrow(res) == 50
  end

  @testset "level 3" begin
    res = book("BTC-USD", 3)
    @test isa(res, DataFrame)
  end

end


text:
module Geometry

#using LinearAlgebra
#using Statistics
using StaticArrays

import Statistics.mean, Statistics.normalize
import LinearAlgebra.cross, LinearAlgebra.norm, LinearAlgebra.dot, LinearAlgebra.I
import LinearAlgebra.det

import Base.+, Base.-, Base.*, Base./, Base.isnan, Base.inv
import Base.getindex, Base.min, Base.max
import Base.convert, Base.isapprox
import Base.size, Base.length, Base.iterate

export VectorLike, Vector3, Normal3, Point3
export Transformation, translation, scaling, rotation, look_at, swaps_handedness
export rotate_z_to
export X_AXIS, Y_AXIS, Z_AXIS, NEG_X_AXIS, NEG_Y_AXIS, NEG_Z_AXIS, IDENTITY_TRANSFORM

abstract type VectorLike end

# Three kinds of vectors are represented with different types: points in space, vectors between
# points, and normal vectors of surfaces. The reason is that these behave differently under
# transformations, and efficient implementation is easier if they are separated.
# The most important difference is that normals are associated with a surface, and when the
# surface is transformed, the normal must be transformed in a way that maintains its normality
# to the surface.

struct Vector3 <: VectorLike
    x::Float64
    y::Float64
    z::Float64
end

struct Point3 <: VectorLike
    x::Float64
    y::Float64
    z::Float64
end

struct Normal3 <: VectorLike
    x::Float64
    y::Float64
    z::Float64
end

const X_AXIS = Vector3(1,0,0)
const Y_AXIS = Vector3(0,1,0)
const Z_AXIS = Vector3(0,0,1)
const NEG_X_AXIS = Vector3(-1,0,0)
const NEG_Y_AXIS = Vector3(0,-1,0)
const NEG_Z_AXIS = Vector3(0,0,-1)

Vector3(x::Real) = Vector3(x,x,x)
Point3(x::Real) = Point3(x,x,x)
Normal3(x::Real) = Normal3(x,x,x)

Vector3(v::Vector3) = Vector3(v.x, v.y, v.z)
Normal3(v::Vector3) = Normal3(v.x, v.y, v.z)
Point3(v::Vector3) = Point3(v.x, v.y, v.z)

size(A::VectorLike) = 3
length(A::VectorLike) = 3

mean(A::Point3, B::Point3) = 0.5 * Point3(A.x+B.x, A.y+B.y, A.z+B.z)
+(A::Point3, B::Normal3) = Point3(A.x+B.x, A.y+B.y, A.z+B.z)
+(A::Point3, B::Vector3) = Point3(A.x+B.x, A.y+B.y, A.z+B.z)
+(B::Vector3, A::Point3) = A+B
+(v::Vector3, u::Vector3) = Vector3(v.x+u.x, v.y+u.y, v.z+u.z)
-(v::Vector3, u::Vector3) = Vector3(v.x-u.x, v.y-u.y, v.z-u.z)
-(v::Point3, u::Point3) = Vector3(v.x-u.x, v.y-u.y, v.z-u.z)
-(A::T) where {T<:VectorLike} = T(-A.x, -A.y, -A.z)

dot(A::VectorLike, B::VectorLike) = A.x*B.x + A.y*B.y + A.z*B.z
cross(A::VectorLike, B::VectorLike) = Vector3(A.y*B.z - A.z*B.y, A.z*B.x - A.x*B.z, A.x*B.y - A.y*B.x)
norm(A::VectorLike) = sqrt(A.x^2 + A.y^2 + A.z^2)
normalize(A::VectorLike) = A / norm(A)

min(A::T, B::T) where {T<:VectorLike} = T(min(A.x,B.x), min(A.y,B.y), min(A.z,B.z))
max(A::T, B::T) where {T<:VectorLike} = T(max(A.x,B.x), max(A.y,B.y), max(A.z,B.z))

function isapprox(A::T, B::T; rtol::Real=sqrt(eps(Float64)), atol::Real=0) where {T<:VectorLike}
    return (isapprox(A.x, B.x; rtol=rtol, atol=atol) &&
    isapprox(A.y, B.y; rtol=rtol, atol=atol) &&
    isapprox(A.z, B.z; rtol=rtol, atol=atol))
end

*(a::Number, v::T) where {T<:VectorLike} = T(a*v.x, a*v.y, a*v.z)
*(v::VectorLike, a::Number) = a*v
/(v::T, a::Number) where {T<:VectorLike} = T(v.x/a, v.y/a, v.z/a)

isnan(v::VectorLike) = isnan(v.x) || isnan(v.y) || isnan(v.z)

getindex(v::VectorLike, i::Integer) = i==1 ? v.x : i==2 ? v.y : i==3 ? v.z : throw(BoundsError())

iterate(V::VectorLike) = (V[1], 2)
iterate(V::VectorLike, state::Integer) = state <= 3 ? (getindex(V, state), state+1) : nothing

convert(::Type{Array{Float64,1}}, v::VectorLike) = [v.x, v.y, v.z]
convert(::Type{Vector3}, x::Array{T,1}) where {T<:Number} = Vector3(x...)
convert(::Type{Normal3}, x::Array{T,1}) where {T<:Number} = Normal3(x...)
convert(::Type{Point3}, x::Array{T,1}) where {T<:Number} = Point3(x...)
convert(::Type{Vector3}, v::Normal3) = Vector3(v.x, v.y, v.z)
convert(::Type{Vector3}, v::Point3) = Vector3(v.x, v.y, v.z)
convert(::Type{Normal3}, v::Vector3) = Normal3(v.x, v.y, v.z)
convert(::Type{Normal3}, v::Point3) = Normal3(v.x, v.y, v.z)
convert(::Type{Point3}, v::Vector3) = Point3(v.x, v.y, v.z)
convert(::Type{Point3}, v::Normal3) = Point3(v.x, v.y, v.z)

####################################
# TRANSFORMATIONS

struct Transformation
    M::SArray{Tuple{4,4},Float64,2,16}
    MInv::SArray{Tuple{4,4},Float64,2,16}
end

const EYE = @SMatrix [1.0 0.0 0.0 0.0 ; 0.0 1.0 0.0 0.0 ; 0.0 0.0 1.0 0.0 ; 0.0 0.0 0.0 1.0]
Transformation() = Transformation(EYE, EYE)
Transformation(M::Array{Float64,2}) = Transformation(SMatrix{4,4}(M), SMatrix{4,4}(inv(M)))
Transformation(M::Array{T,2}) where {T<:Real} = Transformation(convert(SArray{4,4}, M))
inv(T::Transformation) = Transformation(T.MInv, T.M)
*(T::Transformation, U::Transformation) = Transformation(T.M*U.M, U.MInv*T.MInv)
# TODO: check maths above

const IDENTITY_TRANSFORM = Transformation()

swaps_handedness(T::Transformation) = det(view(T.M, 1:3, 1:3)) < 0.0

translation(delta::Vector3) = translation(delta.x, delta.y, delta.z)
function translation(x::Real, y::Real, z::Real)
    M = zeros(Float64, 4, 4)
    M[1,4] = x
    M[2,4] = y
    M[3,4] = z
    M[1,1] = 1.0
    M[2,2] = 1.0
    M[3,3] = 1.0
    M[4,4] = 1.0
    return Transformation(M)
end
function translation(delta::Point3)
    M = Matrix{Float64}(1.0I, 4, 4)
    M[1,4] = delta.x
    M[2,4] = delta.y
    M[3,4] = delta.z
    return Transformation(M)
end

scaling(v::Vector3) = scaling(v.x, v.y, v.z)
function scaling(x::Real, y::Real, z::Real)
    M = zeros(4,4)
    M[1,1] = x
    M[2,2] = y
    M[3,3] = z
    M[4,4] = 1.0
    return Transformation(M)
end

function rotation(axis::Vector3, angle::Number)
    a = normalize(axis)
    s = sin(angle)
    c = cos(angle)
    M = zeros(Float64, 4, 4)
    @inbounds M[1,1] = a.x * a.x + (1.0 - a.x * a.x) * c
    @inbounds M[2,1] = a.x * a.y * (1.0 - c) + a.z * s
    @inbounds M[3,1] = a.x * a.z * (1.0 - c) - a.y * s

    @inbounds M[1,2] = a.y * a.x * (1.0 - c) - a.z * s
    @inbounds M[2,2] = a.y * a.y + (1.0 - a.y * a.y) * c
    @inbounds M[3,2] = a.y * a.z * (1.0 - c) + a.x * s

    @inbounds M[1,3] = a.z * a.x * (1.0 - c) + a.y * s
    @inbounds M[2,3] = a.z * a.y * (1.0 - c) - a.x * s
    @inbounds M[3,3] = a.z * a.z + (1.0 - a.z * a.z) * c

    @inbounds M[4,4] = 1.0

    return Transformation(M, M')
end

function (T::Transformation)(v::Vector3)

    @inbounds return Vector3(
        T.M[1,1]*v.x + T.M[1,2]*v.y + T.M[1,3]*v.z,
        T.M[2,1]*v.x + T.M[2,2]*v.y + T.M[2,3]*v.z,
        T.M[3,1]*v.x + T.M[3,2]*v.y + T.M[3,3]*v.z
    )
end

function (T::Transformation)(v::Point3)
    @inbounds invw = 1.0 / (T.M[4,1]*v.x + T.M[4,2]*v.y + T.M[4,3]*v.z + T.M[4,4])
    @inbounds p = Point3((T.M[1,1]*v.x + T.M[1,2]*v.y + T.M[1,3]*v.z + T.M[1,4]) * invw,
               (T.M[2,1]*v.x + T.M[2,2]*v.y + T.M[2,3]*v.z + T.M[2,4]) * invw,
               (T.M[3,1]*v.x + T.M[3,2]*v.y + T.M[3,3]*v.z + T.M[3,4]) * invw)
    return p
end

function (T::Transformation)(v::Normal3)
    @inbounds return Normal3(
        T.MInv[1,1]*v.x + T.MInv[2,1]*v.y + T.MInv[3,1]*v.z,
        T.MInv[1,2]*v.x + T.MInv[2,2]*v.y + T.MInv[3,2]*v.z,
        T.MInv[1,3]*v.x + T.MInv[2,3]*v.y + T.MInv[3,3]*v.z
    )
end

"""
    rotate_z_to(direction::Vector3)

Return a transformation that rotates the +Z axis to point in the given direction.

"""
function rotate_z_to(direction::Vector3)
    direction = normalize(direction)
    if direction ≈ Z_AXIS
        return Transformation()
    end
    if direction ≈ -Z_AXIS
        return rotation(X_AXIS, π)
    end

    c = direction.z
    s = -sqrt(1.0 - c^2)
    M = zeros(Float64, 4, 4)

    # TODO: this could probably be optimized a little. It's needed fairly often for some
    # radiation pressure camera models.
    a = normalize(Vector3(direction.y, -direction.x, 0.0))

    M[1,1] = a.x * a.x + (1.0 - a.x * a.x) * c
    M[1,2] = a.x * a.y * (1.0 - c) - a.z * s
    M[1,3] = a.x * a.z * (1.0 - c) + a.y * s

    M[2,1] = a.x * a.y * (1.0 - c) + a.z * s
    M[2,2] = a.y * a.y + (1.0 - a.y * a.y) * c
    M[2,3] = a.y * a.z * (1.0 - c) - a.x * s

    M[3,1] = a.x * a.z * (1.0 - c) - a.y * s
    M[3,2] = a.y * a.z * (1.0 - c) + a.x * s
    M[3,3] = a.z * a.z + (1.0 - a.z * a.z) * c

    M[4,4] = 1.0

    return Transformation(M, M')
end

end # module

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# -----------------------------------------------------------------------------
# Types

"""
The supertype for all sky coordinate systems.
"""
abstract type AbstractSkyCoords end

"""
    ICRSCoords(ra, dec)

[International Celestial Reference System](https://en.wikipedia.org/wiki/International_Celestial_Reference_System)

This is the current standard adopted by the International Astronomical Union notably due to its high level of accuracy compared to standard equatorial coordinate systems. What sets this apart from [`FK5Coords`](@ref) is that it is completely defined using extragalactic radio sources rather than a geocentric frame, which means the reference frame will not change due to Earth's motion.

# Coordinates
- `ra` - Right ascension in radians (0, 2π)
- `dec` - Declination in radians (-π/2, π/2)
"""
struct ICRSCoords{T<:AbstractFloat} <: AbstractSkyCoords
    ra::T
    dec::T
    ICRSCoords{T}(ra, dec) where {T<:AbstractFloat} = new(mod2pi(ra), dec)
end
ICRSCoords(ra::T, dec::T) where {T<:AbstractFloat} = ICRSCoords{T}(ra, dec)
ICRSCoords(ra::Real, dec::Real) = ICRSCoords(promote(float(ra), float(dec))...)
ICRSCoords(c::T) where {T<:AbstractSkyCoords} = convert(ICRSCoords, c)
ICRSCoords{F}(c::T) where {F,T<:AbstractSkyCoords} = convert(ICRSCoords{F}, c)

"""
    GalCoords(l, b)

[Galactic Coordinate System](https://en.wikipedia.org/wiki/Galactic_coordinate_system)

This coordinate system is defined based on the projection of the Milky Way galaxy onto our celestial sphere, with (0, 0) being approximately the center of our galaxy.

# Coordinates
- `l` - Galactic longitude in radians (-π, π)
- `b` - Galactic latitude in radians (-π/2, π/2)
"""
struct GalCoords{T<:AbstractFloat} <: AbstractSkyCoords
    l::T
    b::T
    GalCoords{T}(l, b) where {T<:AbstractFloat} = new(mod2pi(l), b)
end
GalCoords(l::T, b::T) where {T<:AbstractFloat} = GalCoords{T}(l, b)
GalCoords(l::Real, b::Real) = GalCoords(promote(float(l), float(b))...)
GalCoords(c::T) where {T<:AbstractSkyCoords} = convert(GalCoords, c)
GalCoords{F}(c::T) where {F,T<:AbstractSkyCoords} = convert(GalCoords{F}, c)

"""
    FK5Coords{equinox}(ra, dec)

[Equatorial Coordinate System](https://en.wikipedia.org/wiki/Equatorial_coordinate_system)

This coordinate system maps the celestial sphere based on a geocentric observer. Historically the oldest, this coordinate system has been shown to be inaccurate due to its definitions based on the Earth, which has long-scale precession causing the reference frame to change. Because of this, an equinox must be provided (typically 2000, commonly known as J2000) which defines the reference frame.

# Coordinates
- `ra` - Right ascension in radians (0, 2π)
- `dec` - Declination in radians (-π/2, π/2)
"""
struct FK5Coords{e,T<:AbstractFloat} <: AbstractSkyCoords
    ra::T
    dec::T
    FK5Coords{e,T}(ra, dec) where {T<:AbstractFloat,e} = new(mod2pi(ra), dec)
end
FK5Coords{e}(ra::T, dec::T) where {e,T<:AbstractFloat} = FK5Coords{e,T}(ra, dec)
FK5Coords{e}(ra::Real, dec::Real) where {e} =
    FK5Coords{e}(promote(float(ra), float(dec))...)
FK5Coords{e}(c::T) where {e,T<:AbstractSkyCoords} = convert(FK5Coords{e}, c)
FK5Coords{e,F}(c::T) where {e,F,T<:AbstractSkyCoords} = convert(FK5Coords{e,F}, c)

# Scalar coordinate conversions
Base.convert(::Type{T}, c::T) where {T<:AbstractSkyCoords} = c

function Base.convert(::Type{T}, c::S) where {T<:AbstractSkyCoords,S<:AbstractSkyCoords}
    r = rotmat(T, S) * coords2cart(c)
    lon, lat = cart2coords(r)
    T(lon, lat)
end

Base.:(==)(a::T, b::T) where {T<:AbstractSkyCoords} = lon(a) == lon(b) && lat(a) == lat(b)


text:
@inline function get_pr_branch_name(package_path)
    project_file_path = find_project_file(package_path)
    pkg = Pkg.Types.read_project(project_file_path)
    return RegistryTools.registration_branch(pkg)
end


text:
int main() {
    // bool operations with doubles
    double y = readDouble();
    printDouble(1.01325e5);
    printDouble(y);
    if (1.01325e5 == y && 1.01325e5 <= y && 1.01325e5 >= y) {
         printInt(1);
    }

    if (0.01325e5 > y && y < 0.01325e5 && y != 0.01325e5) {
    }
    else{
         printInt(2);
    }
    return 0;

}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
# bdf_utils
@inline function U!(k, U)
  @inbounds for r = 1:k
    U[1,r] = -r
    for j = 2:k
      U[j,r] = U[j-1,r] * ((j-1) - r)/j
    end
  end
  nothing
end

function R!(k, ρ, cache)
  @unpack R = cache
  @inbounds for r = 1:k
    R[1,r] = -r * ρ
    for j = 2:k
      R[j,r] = R[j-1,r] * ((j-1) - r * ρ)/j
    end
  end
  nothing
end

# This functions takes help of D2 array to create backward differences array D
# Ith row of D2 keeps Ith order backward differences (∇ⁱyₙ)
function backward_diff!(cache::OrdinaryDiffEqMutableCache, D, D2, k, flag=true)
  flag && copyto!(D[1], D2[1,1])
  for i = 2:k
    for j = 1:(k-i+1)
      @.. D2[i,j] = D2[i-1,j] - D2[i-1,j+1]
    end
    flag && copyto!(D[i], D2[i,1])
  end
end

function backward_diff!(cache::OrdinaryDiffEqConstantCache, D, D2, k, flag=true)
  flag && (D[1] = D2[1,1])
  for i = 2:k
    for j = 1:(k-i+1)
      D2[i,j] = D2[i-1,j] - D2[i-1,j+1]
    end
    flag && (D[i] = D2[i,1])
  end
end

# this function updates backward difference array D when stepsize gets change
# Formula -> D = D * (R * U)
# and it is taken from the paper -
# Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s paper
# E. Alberdi Celaya, J. J. Anza Aguirrezabala, and P. Chatzipantelidis
function reinterpolate_history!(cache::OrdinaryDiffEqMutableCache, D, R, k)
  @unpack tmp = cache
  fill!(tmp,zero(eltype(D[1])))
  for j = 1:k
    for k = 1:k
      @.. tmp += D[k] * R[k,j]
    end
    D[j] .= tmp
    fill!(tmp, zero(eltype(tmp)))
  end
end

function reinterpolate_history!(cache::OrdinaryDiffEqConstantCache, D, R, k)
  tmp = zero(D[1])
  for j = 1:k
    for k = 1:k
      tmp += D[k] * R[k,j]
    end
    D[j] = tmp
  end
end

global const γₖ = @SVector[sum(1//j for j in 1:k) for k in 1:6]

# this stepsize and order controller is taken from
# Implementation of an Adaptive BDF2 Formula and Comparison with the MATLAB Ode15s paper
# E. Alberdi Celaya, J. J. Anza Aguirrezabala, and P. Chatzipantelidis
function stepsize_and_order!(cache, est, estₖ₋₁, estₖ₊₁, h, k)
  zₛ = 1.2
  zᵤ = 0.1
  Fᵤ = 10

  expo = 1/(k+1)
  z = zₛ * ((est)^expo)
  F = inv(z)

  hₖ₋₁ = 0.0
  hₖ₊₁ = 0.0

  if z <= zₛ
    # step is successful
    # precalculations
    if z <= zᵤ
      hₖ = Fᵤ * h
    elseif zᵤ < z <= zₛ
      hₖ = F * h
    end

    if k > 1
      expo = 1/k
      zₖ₋₁ = 1.3 * ((estₖ₋₁)^expo)
      Fₖ₋₁ = inv(zₖ₋₁)
      if zₖ₋₁ <= 0.1
        hₖ₋₁ = 10 * h
      elseif 0.1 < zₖ₋₁ <= 1.3
        hₖ₋₁ = Fₖ₋₁ * h
      end
    end

    expo = 1/(k+2)
    zₖ₊₁ = 1.4 * ((estₖ₊₁)^expo)
    Fₖ₊₁ = inv(zₖ₊₁)

    if zₖ₊₁<= 0.1
      hₖ₊₁ = 10 * h
    elseif 0.1 < zₖ₊₁ <= 1.4
      hₖ₊₁ = Fₖ₊₁ * h
    end
    # adp order and step conditions
    if hₖ₋₁ > hₖ
      hₙ = hₖ₋₁
      kₙ = max(k-1,1)
    else
      hₙ = hₖ
      kₙ = k
    end
    if hₖ₊₁ > hₙ
      hₙ = hₖ₊₁
      kₙ = min(k+1,5)
    end
    if hₙ < h
      hₙ = h
      kₙ = k
    end
    cache.h = hₙ
    cache.order = kₙ
    return true
  else
    # step is not successful
    if cache.c >= 1  # postfail
      cache.h = h/2
      cache.order = k
      return
    end
    if 1.2 < z <= 10
      hₖ = F * h
    elseif z > 10
      hₖ = 0.1 * h
    end
    hₙ = hₖ
    kₙ = k
    if k > 1
      expo = 1/k
      zₖ₋₁ = 1.3 * ((estₖ₋₁)^expo)
      Fₖ₋₁ = inv(zₖ₋₁)
      if 1.3 < zₖ₋₁ <= 10
        hₖ₋₁ = Fₖ₋₁ * h
      elseif zₖ₋₁ > 10
        hₖ₋₁ = 0.1 * h
      end

      if hₖ₋₁ > hₖ
        hₙ = min(h,hₖ₋₁)
        kₙ = max(k-1,1)
      end
    end
    cache.h = hₙ
    cache.order = kₙ
    return false
  end
end


text:
module CellValuesReindex

using Gridap
using Gridap.Helpers
using Gridap.CellValuesGallery
using Gridap.CellMapApply: CellMapFromKernel
using Gridap.CellMapApply: IndexCellMapFromKernel
using Base: @pure

export reindex
export IndexLike

import Base: iterate
import Base: length
import Base: size
import Base: getindex
import Base: IndexStyle

const IndexLike = Union{Integer,CartesianIndex}

function reindex(values::IndexCellValue, indices::Vector{<:IndexLike})
  i = CellValueFromArray(indices)
  reindex(values,i)
end

function reindex(values::IndexCellValue, indices::CellValue{<:IndexLike})
  @notimplemented
end

function reindex(values::IndexCellValue, indices::IndexCellValue{<:IndexLike})
  @notimplemented
end

# Efficient implementations of reindex for concrete types

function reindex(cm::CellMapFromKernel,indices::CellValue{<:IndexLike})
  _reindex_cmfk(cm,indices)
end

function reindex(cm::IndexCellMapFromKernel,indices::CellValue{<:IndexLike})
  _reindex_cmfk(cm,indices)
end

function reindex(cm::IndexCellMapFromKernel,indices::IndexCellValue{<:IndexLike})
  _reindex_cmfk(cm,indices)
end

function _reindex_cmfk(cm,indices)
  k = cm.kernel
  rs = [ reindex(mi,indices) for mi in cm.cellvalues ]
  apply(k,rs...)
end

end # module


text:
export Gradient, RadialGradient

function Gradient(from::Point, to::Point, color1::Union{String, Colors.Colorant}, morecolors...)
    colors_unparsed = (color1, morecolors...)
    colors = [typeof(c) <: String ? parse(Colors.Colorant, c) : c for c in colors_unparsed]
    n = length(colors)
    stops = collect(range(0, 1, length=n))
    Gradient(from, to, stops, colors)
end

function RadialGradient(from::Circle, to::Circle, color1::Union{String, Colors.Colorant}, morecolors...)
    colors_unparsed = (color1, morecolors...)
    colors = [typeof(c) <: String ? parse(Colors.Colorant, c) : c for c in colors_unparsed]
    n = length(colors)
    stops = collect(range(0, 1, length=n))
    RadialGradient(from, to, stops, colors)
end

function RadialGradient(p::Point, r1::Real, r2::Real, color1::Union{String, Colors.Colorant}, morecolors...)
    RadialGradient(Circle(p, r1), Circle(p, r2), color1, morecolors...)
end

function RadialGradient(p::Point, r::Real, color1::Union{String, Colors.Colorant}, morecolors...)
    RadialGradient(Circle(p, 0), Circle(p, r), color1, morecolors...)
end


text:
# Mishchenko2019.jl - reproducing Mishchenko et al. 2019 PRL, with the Feynman
# variational polaron technique

println("Loading PolaronMobility...")
using PolaronMobility
println("Loading Plots for plotting...")
using Plots
gr(size=(500,375))

using Printf
import QuadGK.quadgk # one ring to integrate them all...

#function ImX(nurange,v,w,βred,α,ω,mb)

function savefigs(name)
    println("Saving PNG")
    savefig("$(name).png")
    println("Saving PDF")
    savefig("$(name).pdf")
# high res conversion of vector PDF -> suitable, anti-aliased, PNG for slides
    println("Convert PDF to high res PNG")
    run(`convert -density 300 -resize 1600x $(name).pdf $(name).pdf.png`)
end

# Specify problem
βred=5 # temperature, in thermodynamic Beta units
mb=0.12*9.1093837015E-31
ω=1
const eV = const q = const ElectronVolt = 1.602176487e-19;      # kg m2 / s2
const Boltzmann = const kB =  1.3806504e-23;                  # kg m2 / K s2
const hbar = const ħ = 1.05457162825e-34;          # kg m2 / s

function FHIP1962_Fig2()
    # FHIP1962 reproduction
    # Fig 2 FHIP1962
    α=5
    #in [6] #1:1:10
    v,w = feynmanvw(α)
    println("α: $(α) v: $(v) w: $(w)")

    nu=3.5:0.2:22
    println("Integrating ImX for nu=$nu range...")
    s=ImX(nu, v, w, βred,α,ω,mb)
    plot( s.nu,s.ImX,label="ImX",
             markersize=3,marker=:downtriangle, xlab="nu (units Omega)",ylab="ImX")
    savefigs("FHIP1962_Fig2")
end

function FHIP1962_Fig3()
    # Fig 3 FHIP1962
    α=7
    #in [6] #1:1:10
    v,w = feynmanvw(α)
    println("α: $(α) v: $(v) w: $(w)")

    nu=3.5:0.2:28
    println("Integrating ImX for nu=$nu range...")
    s=ImX(nu, v, w, βred,α,ω,mb)
    plot( s.nu,s.ImX,label="ImX",
             markersize=3,marker=:downtriangle, xlab="nu (units Omega)",ylab="ImX")
end

# Athernmal Feynman variational technique, to compare best to the DiagMC
# results, though nb. they have some kind of T-dep in the diagrams

function Mishchenko_FigureFour()
    # Specify problem
    #mb=1
    βred=8 # temperature, in thermodynamic Beta units
    ω=1

    nurange=0.0:0.05:20
    Trange=[8,4,2,1,0.5,0.25,0.125]

    #using Distributed
    #addprocs(6)
    #@distributed

    for α in [1,2,3,4,5,6,7,8,9,10,12,15,20] #1:1:10
       # athermal action
        v,w = feynmanvw(α)
        plot(xlab="nu (units Omega) alpha=$(α)",ylab="Mob", ylims=(0,1.0E12))

        for T in Trange
            βred=1/T
            @time s=ImX(nurange, v, w, βred,α,ω,mb)
            plot!( s.nu, s.μ, label="T=$(T)")
        end

        savefigs("Mishchenko-Fig4-AthermalAction_alpha_$(α)")

        plot(xlab="nu (units Omega) alpha=$(α)",ylab="Mob", ylims=(0,3.0E12))

        for T in Trange
            βred=1/T
            # Osaka finite temperature action
            v,w=feynmanvw(α, βred, verbose=true) # temperature dependent Action

            @time s=ImX(nurange, v, w, βred,α,ω,mb)
            plot!( s.nu, s.μ, label="T=$(T) v=$(round(v,sigdigits=2)) w=$(round(w,sigdigits=2))")
        end

        savefigs("Mishchenko-Fig4-OsakaFiniteTemperatureAction_alpha_$(α)")
    end
end
#wait()

function HellwarthMobility(v,w,βred, α, ω)
# Hellwarth1999 - directly do contour integration in Feynman1962, for
        # finite temperature DC mobility
        # Hellwarth1999 Eqn (2) and (1) - These are going back to the general
        # (pre low-T limit) formulas in Feynman1962.  to evaluate these, you
        # need to do the explicit contour integration to get the polaron
        # self-energy
        R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)

        b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)

        a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
        k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
        K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)

        #Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
        RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
        #μ=RHS^-1 * (q)/(ω*mb)
        #@printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
        RHS^-1
end

function Mischenko_MobilityComparison(α; Trange=0.1:0.01:25)
    println("Mishchenko Feynman athermal Hellwarth explicit integration mobility with α= $(α)")
    βred=8
    ω=1 #(50*kB)/ħ # in artifical units of T
    fout=open("MishchenkoMobility_$(α).dat","w")
    @printf(fout,"# T v w mob(SI)")

    for T in Trange
        β=ω/T #  ħ*ω/(kB*T)

        # Athermal variational
        v,w = feynmanvw(α)
        μ=HellwarthMobility(v,w, β, α, ω)
        @printf(fout,"%f    %f %f %f ",T,v,w,μ)

        # Thermal theory
        v,w = feynmanvw(α, β)
        μ=HellwarthMobility(v,w, β, α, ω)
        @printf(fout,"%f %f %f ",v,w,μ)
        @printf(fout,"\n")
    end
    close(fout)
end

for α in [2.4, 4, 6, 8, 10, 12, 14] # Mishchenko2019 Fig 2 and 3
    Mischenko_MobilityComparison(α)
end

println("That's me!")


text:
@testset "Maps and Reductions" begin
    @testset "sum" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                x = randn(T, sizes[1:N]...)
                test_frule(sum, x; fkwargs=(;dims=dims))
                test_rrule(sum, x; fkwargs=(;dims=dims))
            end
        end
    end  # sum

    @testset "sum abs2" begin
        sizes = (3, 4, 7)
        @testset "dims = $dims" for dims in (:, 1)
            @testset "Array{$N, $T}" for N in eachindex(sizes), T in (Float64, ComplexF64)
                x = randn(T, sizes[1:N]...)
                test_frule(sum, abs2, x; fkwargs=(;dims=dims))
                test_rrule(sum, abs2 ⊢ DoesNotExist(), x; fkwargs=(;dims=dims))
            end
        end
    end  # sum abs2
end


text:
using StaticArrays, OrdinaryDiffEq, BenchmarkTools
tspan = (0.0,10.0)
println("Testing speed of solving ODEs\n")
println("---Using StaticArrays: f(t, u) = du::SVector")
#--------------------#
@inline eom(t, u) = SVector{5}(10.0(u[2]-u[1]), u[1]*(28.0-u[3]) - u[2],
u[1]*u[2] - (8/3)*u[3], 2.0u[4]*(1-u[4]), 0.15*u[4]*(1-u[5]))

# println("ODEProblem creation:")
# display(@benchmark ODEProblem(eom,SVector{3}(rand(3)),$tspan))
prob = ODEProblem(eom,SVector{5}(10rand(5)),tspan)
println("Solution (for total time 10)")
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
display(@benchmark solve($prob,Tsit5(),reltol=1e-8,abstol=1e-8))
sleep(0.1)

println("\n---Using Base Arrays: f(t, u) = du::Vector")
#--------------------#
@inline eomB(t, u) = [10.0(u[2]-u[1]), u[1]*(28.0-u[3]) - u[2],
u[1]*u[2] - (8/3)*u[3], 2.0u[4]*(1-u[4]), 0.15*u[4]*(1-u[5])]

# println("ODEProblem creation:")
# display(@benchmark ODEProblem(eom,rand(3),$tspan))
prob = ODEProblem(eomB,10rand(5),tspan)
println("Solution (for total time 10)")
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
display(@benchmark solve($prob,Tsit5(),reltol=1e-8,abstol=1e-8))
sleep(0.1)

println("\n---Using Base & in-place: f!(t, u, du)::Vector")

@inline function eom!(t, u, du)
    du[1] = 10.0(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
    du[4] = 2.0u[4]*(1-u[4])
    du[5] = 0.15*u[4]*(1-u[5])
end

# println("ODEProblem creation:")
# display(@benchmark ODEProblem(eom!,rand(3),$tspan))
prob = ODEProblem(eom!,10rand(5),tspan)
println("Solution (for total time 10)")
sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
display(@benchmark solve($prob,Tsit5(),reltol=1e-8,abstol=1e-8))
sleep(0.1)

# println("\n---Using MVector & in-place: f!(t, u, du)::MVector")
# @inline function eomM!(t, u::MVector, du::MVector)
#     du[1] = 10.0(u[2]-u[1])
#     du[2] = u[1]*(28.0-u[3]) - u[2]
#     du[3] = u[1]*u[2] - (8/3)*u[3]
# end
#
# probM = ODEProblem(eomM!,MVector{3}(10rand(3)),tspan)
# println("Solution (for total time 1)")
# sol = solve(probM,Tsit5(),reltol=1e-8,abstol=1e-8)
# display(@benchmark solve($probM,Tsit5(),reltol=1e-8,abstol=1e-8))
# # #this doesnt work currently


text:
# Write tests
# Batch normalization may be needed once more food is collected and traded
module Nature

using ReinforcementLearning
using IntervalSets
using SparseArrays
using Distributions
using LinearAlgebra
using StableRNGs
using Random
using Flux
using Infiltrator
#using GLMakie
using DataStructures
using Serialization
using Logging
using TensorBoardLogger
import Base: size
import ReinforcementLearning.RLZoo: EnrichedAction

export NatureEnv, FoodGen, reset!, remove!, visualize,
step_through_env, MultiPPOManager, NatureHook, ($),
visualize_food, change_run, build_MultiPPOManager, PICK_FOOD_1, PLACE_FOOD_1

include("./structs.jl")
include("./utils.jl")
include("./act.jl")
include("./player.jl")
include("./food.jl")
include("./natureenv.jl")
include("./multippo.jl")
include("./build-model.jl")
include("./debug.jl")
#include("./visualize.jl")
include("./hooks.jl")

end


text:
export MaskedReLU

"""
$(TYPEDEF)

Represents a masked ReLU activation, with `mask` controlling how the ReLU is applied to
each output.

`p(x)` is shorthand for [`masked_relu(x, p.mask)`](@ref) when `p` is an instance of
`MaskedReLU`.

## Fields:
$(FIELDS)
"""
struct MaskedReLU{T<:Real} <: Layer
    mask::Array{T}
    tightening_algorithm::Union{TighteningAlgorithm,Nothing}
end

function MaskedReLU(mask::Array{T}) where {T<:Real}
    MaskedReLU{T}(mask, nothing)
end

function MaskedReLU(mask::Array{T}, ta::TighteningAlgorithm) where {T<:Real}
    MaskedReLU{T}(mask, ta)
end

function Base.show(io::IO, p::MaskedReLU)
    num_zeroed_units = count(p.mask .< 0)
    num_passthrough_units = count(p.mask .> 0)
    num_rectified_units = length(p.mask) - num_zeroed_units - num_passthrough_units
    print(
        io,
        "MaskedReLU with expected input size $(size(p.mask)). ($(num_zeroed_units) zeroed, $(num_passthrough_units) as-is, $(num_rectified_units) rectified).",
    )
end

(p::MaskedReLU)(x::Array{<:Real}) = masked_relu(x, p.mask)
(p::MaskedReLU)(x::Array{<:JuMPLinearType}) = (Memento.info(MIPVerify.LOGGER, "Applying $p ... ");
masked_relu(x, p.mask, nta = p.tightening_algorithm))


text:
# The get and set functions pertaining to embeddings

export getxy, hasxy
export set_vertex_color, set_line_color, get_vertex_color, get_line_color
export set_vertex_size, get_vertex_size

"""
`getxy(G)` returns (a copy of) the `xy`-embedding of `G`.

`getxy(G,v)` returns the `xy`-coordinates of the vertex `v`.
"""
function getxy(G::SimpleGraph)
    ensure_embed(G)
    return (deepcopy(G.cache[:xy]))
end

function getxy(G::SimpleGraph, v)
    ensure_embed(G)
    return G.cache[:xy][v]
end

"""
`hasxy(G::SimpleGraph)` returns `true` if an embedding has been
given to this graph.
"""
hasxy(G::SimpleGraph)::Bool = cache_check(G, :xy)

"""
`get_line_color(G)` returns the color for edges and vertex boundaries
"""
function get_line_color(G::SimpleGraph)
    ensure_embed(G)
    return G.cache[:line_color]
end

"""
`set_line_color(G::SimpleGraph, hue=:black)` sets the color of the graph's
edges and vertex boundaries.
"""
function set_line_color(G::SimpleGraph, hue = :black)
    ensure_embed(G)
    G.cache[:line_color] = hue
end

"""
`get_vertex_color(G,v)` returns the color assigned to vertex `v`.

`get_vertex_color(G)` returns a copy of the dictionary mapping vertices to colors.
"""
function get_vertex_color(G::SimpleGraph, v)
    ensure_embed(G)
    return G.cache[:vcolor][v]
end

function get_vertex_color(G::SimpleGraph)
    ensure_embed(G)
    return deepcopy(G.cache[:vcolor])
end

"""
`set_vertex_color(G,v,hue)` sets the color of vertex `v` to `hue`.

`set_vertex_color(G,hue)` sets the color of all vertices to `hue`.
If `hue` is omitted, we use `:white`.
"""
function set_vertex_color(G::SimpleGraph, v, hue)
    ensure_embed(G)
    G.cache[:vcolor][v] = hue
end

function set_vertex_color(G::SimpleGraph, hue)
    ensure_embed(G)
    d = G.cache[:vcolor]
    for v in G.V
        d[v] = hue
    end
end

set_vertex_color(G::SimpleGraph) = set_vertex_color(G, :white)

"""
`set_vertex_color(G::SimpleGraph, d::Dict, palette)` where `d` is a dictionary
mapping vertices to integers and `palette` is a list of colors.

Convert a mapping of vertices to integers into colors for the vertices.

Vertices are assigned colors as follows: vertex `v` gets color `palette[k]`
where `k=d[v]`. If `palette` is omitted, use the constant global variable
`colorize_hues`.
"""
function set_vertex_color(
    G::SimpleGraph,
    d::Dict{S,T},
    palette = colorize_hues,
) where {S,T<:Integer}
    colorize(G, d, palette)
end

"""
`set_vertex_size(G,s)` sets the radius of the graph's vertices.

`set_vertex_size(G)` restores the radius to the default value.
"""
function set_vertex_size(G::SimpleGraph, s::Int = DEFAULT_MARKER_SIZE)
    ensure_embed(G)
    G.cache[:vsize] = s
end

function get_vertex_size(G::SimpleGraph)
    ensure_embed(G)
    return G.cache[:vsize]
end


text:
__precompile__()

module Hospitalplanning
using JuMP
using Gurobi
using MathOptInterface
using Random
using Dates
using Query
using XLSX
using Random
using JuliaDB
using CSV
using DataFrames
using CSVFiles

include("Misc.jl")

include("GenerateSampleData.jl")
include("readwritedata.jl")
include("MIP/datastructures.jl")
include("MIP/SupportFunctions.jl")
include("MIP/masterproblem.jl")
include("MIP/subproblem.jl")
include("MIP/model_column.jl")

end # module


text:
using GLM, Statistics
import GLM: linkinv, linkfun
import Distributions: UnivariateDistribution

# root type: use this to define an interface shared by
# all instances of AbstractTrait
abstract type AbstractTraitModel end

__default_behavior(trait) = error("function not supported by $(typeof(trait))")

"Check the number of traits."
ntraits(trait::AbstractTraitModel) = __default_behavior(trait)

"Check the number of variance components."
nvc(trait::AbstractTraitModel) = __default_behavior(trait)

"Check the number of outcome categories."
noutcomecategories(trait::AbstractTraitModel) = __default_behavior(trait)

"Check the number of fixed covariates in linear model."
neffects(trait::AbstractTraitModel) = size(trait.X, 2)

"Check the number of people in the sample."
nsamplesize(trait::AbstractTraitModel) = size(trait.X, 1)

struct GLMTrait{distT, linkT, vecT1, vecT2, vecT3, matT, matT2} <: AbstractTraitModel
	X::matT             # all effects
	β::vecT1            # regression coefficients
	G::matT2        # genotypes CHANGE TO SNPARRAY
	γ::vecT1             # effect sizes for each snp
	η::vecT2            # linear predictor η = X*β
	μ::vecT3            # expected value of response μ = g^-1(η), g is the link function
	dist::Type{distT}   # univariate, exponential family of distributions
	link::linkT         # link function g(μ) = X*β = η
	function GLMTrait(X::matT, β::vecT1, G::matT2, γ::vecT1, η::vecT2, distribution::D, link::linkT; lb = -20, ub = 20) where {D, linkT, vecT1, vecT2, matT, matT2}
	# extract the base type without type parameters
	distT = Base.typename(typeof(distribution)).wrapper
	n = size(X, 1)
	μ = zeros(n)
	for i in eachindex(η)
	μ[i] = GLM.linkinv(link, η[i])
	end
	# make a new instance
	new{distT, linkT, vecT1, typeof(η), typeof(μ), matT, matT2}(X, β, G, γ, η, μ, distT, link)
	end
end

function GLMTrait(X::Matrix{T}, β::Vector{T}, G::SnpArray, γ::Vector{T},
	 distribution::D, link::linkT; lb = -20, ub = 20) where {D, linkT, T <: BlasReal}
	 distT = Base.typename(typeof(distribution)).wrapper
 	 n = size(X, 1)
	 η = Vector{T}(undef, n)
 	 clamp_eta!(η, X, β, distT; lb = lb , ub = ub)
	 genovec = zeros(Float64, size(G))
 	 Base.copyto!(genovec, @view(G[:, :]), model = ADDITIVE_MODEL, impute = true)
  return GLMTrait(X, β, genovec, γ, η, distribution, link)
end

function GLMTrait(X::Matrix{T}, β::Vector{T}, distribution::D, link::linkT; lb = -20, ub = 20) where {D, linkT, T <: BlasReal}
	distT = Base.typename(typeof(distribution)).wrapper
	n = size(X, 1)
	η = zeros(n)
	clamp_eta!(η, X, β, distT; lb = lb , ub = ub)
	genovec = zeros(Float64, n)
	copyto!(genovec, X[:, end])
    GLMTrait(X, β, genovec, [β[end]], η, distribution, link)
end

# clamp the trait values
  # specific to Poisson and NegativeBinomial
  function clamp_eta!(η::AbstractVecOrMat, X::Matrix{T}, β::AbstractVecOrMat,
	   distT::Union{Type{Poisson}, Type{NegativeBinomial}}; lb = -20, ub = 20) where T <: BlasReal
	  mul!(η, X, β)
	  η .= map(y -> y >= ub ? ub : y, η)
  end

  # specific to Bernoulli and Binomial
  function clamp_eta!(η::AbstractVecOrMat, X::Matrix{T}, β::AbstractVecOrMat,
	   distT::Union{Type{Bernoulli}, Type{Binomial}}; lb = -20, ub = 20) where T <: BlasReal
	  mul!(η, X, β)
	  clamp!(η, lb, ub)
	  η
  end

  # all others
  function clamp_eta!(η::AbstractVecOrMat, X::Matrix{T}, β::AbstractVecOrMat,
	   distT::Type{D}; lb = -20, ub = 20) where {D, T <: Real}
	  mul!(η, X, β)
	  η
  end

# better printing; customize how a type is summarized in a REPL
function Base.show(io::IO, x::GLMTrait)
    print(io, "Generalized Linear Model\n")
    print(io, "  * response distribution: $(x.dist)\n")
    print(io, "  * link function: $(typeof(x.link))\n")
    print(io, "  * sample size: $(nsamplesize(x))")
	print(io, "  * fixed effects: $(neffects(x))")
end

struct OrderedMultinomialTrait{matT, vecT1, vecT2, linkT} <: AbstractTraitModel
    X::matT             # all effects
    β::vecT1            # regression coefficients
    θ::vecT2            # must be increasing
	link::linkT	        # link function from GLM.jl
	# create a new instance
	function OrderedMultinomialTrait(X::matT, β::vecT1,
		 θ::vecT2, link::linkT)  where {matT, vecT1, vecT2, linkT}
    return new{matT, vecT1, vecT2, linkT}(X, β, θ, link)
  end
end

function Base.show(io::IO, x::OrderedMultinomialTrait)
    print(io, "Ordinal Multinomial Model\n")
    print(io, "  * number of fixed effects: $(neffects(x))\n")
	print(io, "  * number of ordinal multinomial outcome categories: $(noutcomecategories(x))\n")
    print(io, "  * link function: $(typeof(x.link))\n")
    print(io, "  * sample size: $(nsamplesize(x))")
end

noutcomecategories(trait::OrderedMultinomialTrait) = length(trait.θ) + 1

"""
VCMTrait
VCMTrait object is a model framework object that stores information about the simulation model
 of multiple traits, under the Variance Component Model Framework.
"""
struct VCMTrait{T <: Real} <: AbstractTraitModel
	X::Matrix{T}             # all effects
	β::Matrix{T}             # regression coefficients
	G::AbstractVecOrMat        # genotypes CHANGE TO SNPARRAY
	γ::Matrix{T}             # effect sizes for each snp
	vc::Vector{VarianceComponent}
	μ::Matrix{T}            # evaluated fixed effect
	Z::Matrix{T}
	function VCMTrait(X, β, G, γ, vc, μ)
		T = eltype(μ)
		Z  = Matrix{T}(undef, size(μ))
		new{T}(X, β, G, γ, vc, μ, Z)
	end
end

function VCMTrait(X::Matrix{T}, β::Matrix{T}, G::SnpArray, γ::Matrix{T},
	 vc::Vector{VarianceComponent}) where T <: BlasReal
	n, p, m, d = size(X, 1), size(X, 2), length(vc), size(β, 2)
	genovec = SnpLinAlg{Float64}(G, model=ADDITIVE_MODEL, center=true, scale=true);
	μ = zeros(n, d)
	μ_null = zeros(n, d)
	LinearAlgebra.mul!(μ_null, X, β)
	mul!(μ, genovec, γ)
	BLAS.axpby!(1.0, μ_null, 1.0, μ)
	# constructor
	VCMTrait(X, β, genovec, γ, vc, μ)
end

function VCMTrait(X::Matrix{T}, β::Matrix{T}, G::SnpArray,
	 γ::Matrix{T}, Σ::Vector{Matrix{T}}, V::Vector{Matrix{T}}) where T <: BlasReal
	vc = [VarianceComponent(Σ[i], V[i]) for i in 1:length(V)]
	return VCMTrait(X, β, G, γ, vc)
end

function VCMTrait(X::Matrix{T}, β::Matrix{T},
	 vc::Vector{VarianceComponent}) where T <: BlasReal
	n, d = size(X, 1), size(β, 2)
	μ = Matrix{T}(undef, n, d)
	mul!(μ, X, β)
	genovec = zeros(Float32, n, 1)
	copyto!(genovec, X[:, end])
	γ = Matrix{T}(undef, 1, d)
	copyto!(γ, β[end, :])
	return VCMTrait(X, β, genovec, γ, vc, μ)
end

function VCMTrait(X::Matrix{T}, β::Matrix{T},
	 Σ::Vector{Matrix{T}}, V::Vector{Matrix{T}}) where T <: BlasReal
	vc = [VarianceComponent(Σ[i], V[i]) for i in 1:length(V)]
	return VCMTrait(X, β, vc)
end

function VCMTrait(formulas::Vector{String}, df::DataFrame,
	 vc::Vector{VarianceComponent})
  d = length(formulas)
  n = size(df)[1]
  μ = zeros(n, d)
  found_covariates = Symbol[]
  for i in 1:d
    #calculate the mean vector
	μ[:, i], found_markers = mean_formula(formulas[i], df)
	union!(found_covariates, found_markers)
  end
  X = Matrix(df[:, found_covariates])
  p = length(found_covariates)
  β = Matrix{Float32}(undef, p, d)
  G = Matrix{Float32}(undef, n, 1)
  γ = Matrix{Float32}(undef, 1, p)
  return VCMTrait(X, β, G, γ, vc, μ)
end

function VCMTrait(formulas::Vector{String}, df::DataFrame,
	 Σ::Vector{Matrix{T}}, V::Vector{Matrix{T}}) where T <: BlasReal
	vc = [VarianceComponent(Σ[i], V[i]) for i in 1:length(V)]
	VCMTrait(formulas, df, vc)
end

#  Variance Component Model
function Base.show(io::IO, x::VCMTrait)
    print(io, "Variance Component Model\n")
    print(io, "  * number of traits: $(ntraits(x))\n")
	print(io, "  * number of variance components: $(nvc(x))\n")
    print(io, "  * sample size: $(nsamplesize(x))")
end

#make our new type implement the interface defined above
ntraits(trait::VCMTrait) = size(trait.vc[1].Σ, 1)
nvc(trait::VCMTrait) = length(trait.vc)

"""
GLMMTrait
GLMMTrait object is a model framework that stores information about the simulation of multiple traits, under the Generalized Linear Mixed Model Framework.
"""
struct GLMMTrait{distT, linkT, T <: BlasReal} <: AbstractTraitModel
   X::Matrix{T}             # all effects
   β::Matrix{T}            # regression coefficients
   μ::Matrix{T}            # mean of the glmm with random effects
   η::Matrix{T}
   Z::Matrix{T}             # place holder for aggregating random effects
   Y_vcm::Matrix{T}        # place holder for getting glmmm mean for simulation
   vc::Vector{VarianceComponent}
   dist::Type{distT}   # univariate, exponential family of distributions
   link::linkT         # link function g(μ) = X*β
   function GLMMTrait(X::Matrix{T}, β::Matrix{T}, vc::Vector{VarianceComponent},
   	 distribution::D, link::linkT; lb = -20, ub = 20) where {D, linkT, T<: BlasReal}
   	distT = Base.typename(typeof(distribution)).wrapper
   	Z = zeros(size(X, 1), size(β, 2))
   	Y_vcm = zeros(size(X, 1), size(β, 2))
    η = zeros(size(X, 1), size(β, 2))
   	clamp_eta!(η, X, β, distT; lb = lb , ub = ub)
   	μ = zero(η)
   	for i in eachindex(η)
   		μ[i] = GLM.linkinv(link, η[i])
   	end
     return new{distT, linkT, T}(X, β, μ, η, Z, Y_vcm, vc, distT, link)
   end
end

# better printing; customize how a type is summarized in a REPL
function Base.show(io::IO, x::GLMMTrait)
    print(io, "Generalized Linear Mixed Model\n")
    print(io, "  * response distribution: $(x.dist)\n")
    print(io, "  * link function: $(typeof(x.link))\n")
	print(io, "  * number of variance components: $(nvc(x))\n")
    print(io, "  * sample size: $(nsamplesize(x))")
end
# make our new type implement the interface defined above
ntraits(trait::GLMMTrait) = size(trait.vc[1].Σ, 1)
nvc(trait::GLMMTrait) = length(trait.vc)