text:
using Genie, Genie.Renderer.Html, Stipple, StipplePlotly, Colors, ColorSchemes

Genie.config.log_requests = false

function rgb(pix::RGB{Float64})
  R = round(Int, 255 * clamp(pix.r, 0.0, 1.0))
  G = round(Int, 255 * clamp(pix.g, 0.0, 1.0))
  B = round(Int, 255 * clamp(pix.b, 0.0, 1.0))
  return "rgb($R,$G,$B)"
end

# See: https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/
function ColorScale(scheme::Symbol, N = 101)
  x = permutedims(0.0:(1.0/(N - 1)):1.0)
  cs = get(colorschemes[scheme], x, :clamp)
  cs_rgb = rgb.(cs)
  return vcat(x, cs_rgb)
end

function to_plotly_standard(x)
  if isnothing(x)
    m = "null"
  elseif ismissing(x)
    m = "null"
  elseif isa(x, AbstractString)
    m = x
  elseif isinf(x)
    m = "null"
  elseif isnan(x)
    m = "null"
  else
    m = x
  end
  return m
end

xx = ["L0", "L63", "L128", "L191", "L255"]
yy = ["L255", "L128", "L0"]
zz = [ 100.0 0.0  missing;
       31.30235455584154  94.9495352852204 0.7150792978758869;
       63.172349721957175 37 -72.02122607433563;
       51.19876097916769  18.676411380287483 314.15;
       87.03992589163836  88.74217734803698 62.122844059857215];

cs = :algae # color scheme

function textcolor_mask(z,zmin,zmax,cs; halfluminance=(0.5^2.2))
  m = zeros(Float64, size(z))
  for i=1:length(m)
    if !ismissing(z[i])
      c = get(colorschemes[cs], Float64(z[i]), (zmin,zmax))
      if xyY(c).Y < halfluminance
        m[i] = 1.0
      end
    end
  end
  m
end

function annotate_map(x, y, z, zmin, zmax, cs)
  m = textcolor_mask(z,zmin,zmax,cs; halfluminance=(0.5^2.2))
  anv = Vector{PlotAnnotation}(undef, 0)
  anv = PlotAnnotation[]
  for i=1:size(z,1)
    for j=1:size(z,2)
      c = rgb(RGB(m[i,j], m[i,j], m[i,j]))
      an = PlotAnnotation(visible=!ismissing(z[i,j]), x=x[i], y=y[j], xref="x", yref="y", text= "$(round(z[i,j]; sigdigits=3))", ax=0, ay=0, showarrow=false, font=Font(color=c, size=16) )
      push!(anv, an)
    end
  end
  anv
end

anv = annotate_map(xx, yy, zz, 0.0, 100.0, cs);

pd(name, x, y, z; zmin = 0.0, zmax = 100.0, cs = cs) = PlotData(
  z =  to_plotly_standard.(round.(z; sigdigits=3)),
  zmin = zmin,
  zmax = zmax,
  x = x,
  y = y,
  plot = StipplePlotly.Charts.PLOT_TYPE_HEATMAP,
  name = name,
  colorscale = ColorScale(cs),
  colorbar = ColorBar("Z-data", 18, "right"),
  hoverongaps = false,
  hoverinfo = "x+y+z"
)

pl(title; annotations::Union{Nothing, Vector{PlotAnnotation}} = nothing) = PlotLayout(
  plot_bgcolor = "#FFFFFF",
  title = PlotLayoutTitle(text=title, font=Font(24)),
  margin_b = 25,
  margin_t = 80,
  margin_l = 60,
  margin_r = 40,
  xaxis = [PlotLayoutAxis(xy = "x", index = 1,
    title = "From",
    font = Font(18),
    ticks = "outside top",
    side = "top",
    position = 1.0,
    showline = true,
    showgrid = false,
    zeroline = false,
    mirror = "all",
    ticklabelposition = "outside top")],
  yaxis = [PlotLayoutAxis(xy = "y", index = 1,
    showline = true,
    zeroline = false,
    mirror = "all",
    showgrid = false,
    title = "To",
    font = Font(18),
    ticks = "outside",
    scaleanchor = "x",
    scaleratio = 1,
    constrain = "domain",
    constraintoward = "top")],
  annotations = annotations
)

Base.@kwdef mutable struct Model <: ReactiveModel
  data::R{Vector{PlotData}} = [pd("Random 1", xx, yy, zz)], READONLY
  layout::R{PlotLayout} = pl(""; annotations=anv), READONLY
  config::R{PlotConfig} = PlotConfig(), READONLY
end

model = Stipple.init(Model())

function ui()
  page(
    vm(model), class="container", [
      plot(:data, layout = :layout, config = :config)
    ]
  ) |> html
end

route("/", ui)

up()


text:
using FileIO,Optim
include("model.jl")
img = load("raw/509.tiff")
x0 = Float32[290 265 385; 310 360 315] / 512
ygold = convert(Array{Float32}, img.data)
#ygold[ygold.<0.05] = 0.0
x = similar(ygold,(2,3))
x6 = reshape(x, 6)
minit(x)
ypred = similar(ygold)
dx = similar(x)

f(x)=Float32(mlossxy(reshape(x,(2,3)),ygold))

function g(x,dx)
    x33=reshape(x,(2,3))
    dx33=reshape(dx,(2,3))
    mforw(x33,ypred)
    mback(x33,ypred,ygold,dx33)
    dx33
end
# optimize(f, g, x, LBFGS())

ext: jl

max_stars_repo_name: JuliaTagBot/NeuronTracking

max_issues_repo_name: JuliaTagBot/NeuronTracking

max_forks_repo_name: JuliaTagBot/NeuronTracking


text:
using Test
include("Utilities.jl")

@testset "Utilities" begin
    G1 = MatrixNetwork(load_matrix_network("cores_example")) # Simple networks
    G2 = MatrixNetwork([1,2,2,3,3], [2,4,5,4,5])
    uG2 = make_undirected(G2)

    M1 = MatrixNetwork([1,1,3,3,3,3,3,4,4,6,6,6], [3,1,1,2,2,4,6,6,6,1,2,7])
    M2 = MatrixNetwork([1,1,2,2,2,3,3,4,4,5,5,6,6,6,7,7,7],
                       [2,6,1,1,4,2,5,1,3,3,4,1,6,7,4,4,5])

    @testset "Degree Distribution" begin
        for i = 2:10 # Testing common subgraphs: path, cycle, clique
            @test deg_distr(make_undirected(path(i))) == [1; [2 for j = 1:i-1]; 1]
            @test deg_distr(make_undirected(cycle(i+1))) == [2 for j = 1:i+1]
            @test deg_distr(clique(i)) == [i-1 for j = 1:i]
        end

        @test deg_distr(G1) == [1,4,2,6,5,3,3,4,3,5,4,1,1,3,3,0,2,2,2,3,1]
        @test deg_distr(uG2) == [1,3,2,2,2]
    end

    @testset "Nodes by Degree" begin # Getting sorted list of nodes by their degree
        @test nodes_by_deg(G1) == [16, 1, 12, 13, 21, 3, 17, 18, 19, 6, 7, 9, 14, 15, 20, 2, 8, 11, 5, 10, 4]
        @test nodes_by_deg(G1; descending=true) == [4, 5, 10, 2, 8, 11, 6, 7, 9, 14, 15, 20, 3, 17, 18, 19, 1, 12, 13, 21, 16]

        @test nodes_by_deg(uG2) == [1, 3, 4, 5, 2]
        @test nodes_by_deg(uG2; descending=true) == [2, 3, 4, 5, 1]
    end

    @testset "Get Edges" begin
        function edge_compare(A::MatrixNetwork)
            ei, ej = directed_edges(A)
            return get_edges(A) == [[ei[x],ej[x]] for x = 1:length(ei)]
        end
        @test edge_compare(G1)
        @test edge_compare(G2)
        @test edge_compare(uG2)
        @test edge_compare(M1)
        @test edge_compare(M2)
    end

    @testset "Simplify" begin
        simplify_test_simple(A::MatrixNetwork) = get_edges(simplify(A)) == get_edges(A) # Checks that simplifying a simple graph keeps it simple
        @test simplify_test_simple(G1)
        @test simplify_test_simple(G2)

        function simplify_test(A::MatrixNetwork)
            M = Array(sparse(simplify(A)))
            return Diagonal(M) == zeros(Int64, (A.n, A.n)) && minimum(M) == 0 && maximum(M) == 1
        end

        @test simplify_test(M1)
        @test simplify_test(M2)
    end

    @testset "Directed -> Undirected" begin
        function undirected_test(A::MatrixNetwork)
            U = sparse(make_undirected(A))
            return U == U' == sparse(A) + sparse(A)'
        end

        @test undirected_test(G1)
        @test undirected_test(G2)
        @test undirected_test(uG2)

        @test undirected_test(M1)
        @test undirected_test(M2)
    end

    @testset "Clustering Coefficients" begin
        @test clustering(G1) == (0.3555555555555555, 0.44)
        @test_throws ErrorException clustering(G2)

        @test clustering(make_undirected(M1)) == (0.5898904756952629, 0.585695488070208)
        @test clustering(make_undirected(M2)) == (0.4950303252150675, 0.48426879640604437)
    end

end


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

using Test

# interpreted but inferred/optimized top-level expressions with vars
let code = """
           while true
               try
                   this_is_undefined_29213
                   ed = 0
                   break
               finally
                   break
               end
           end
           print(42)
           """
    @test read(`$(Base.julia_cmd()) --startup-file=no --compile=min -e $code`, String) == "42"
end

let code = "Threads.atomic_add!(Threads.Atomic{Int}(40), 2)"
    @test read(`$(Base.julia_cmd()) --startup-file=no --compile=min -E $code`, String) == "40\n"
end

let p = Pipe(),
    c = pipeline(`$(Base.julia_cmd()) --startup-file=no --compile=min -E 'error()'`, stderr=p)
    proc = run(c, wait=false)
    readline(p)
    @test readline(p) == "Stacktrace:"
    wait(proc)
    close(p)
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: AsForMe/julia

max_forks_repo_name: AsForMe/julia


text:
module SfPlanetary

using Unitful, UnitfulAstro, UnitfulAngles
using ..SfGeometry

export AbstractBody, Body, Orbit, Rotation, Star
export satellites

import ..SfUnits: Angle, to_angle

@derived_dimension ThermalFlux Unitful.𝐌*Unitful.𝐓^-3

abstract type AbstractBody end

const satellites_of = Dict{AbstractBody, Vector{AbstractBody}}()

"""
    satellites(body::AbstractBody)

Return a list of sattelites associated with `body`, which may be empty.
"""
satellites(body::AbstractBody) = haskey(satellites_of, body) ? satellites_of[body] : Vector{AbstractBody}()

"""
    add_satellite(primary::AbstractBody, satellite::AbstractBody)

Associate `satellite` with `primary`. Should not usually be called by user code.
"""
function add_satellite(primary::AbstractBody, satellite::AbstractBody)
	if !haskey(satellites_of, primary)
		satellites_of[primary] = [ satellite ]
	else
		push!(satellites_of[primary], satellite)
	end
end

struct Orbit
	parent::AbstractBody
	semi_major_axis::Unitful.Length
	eccentricity::Real
	inclination::Angle # relative to equator of parent body
	ascending_node::Union{Missing, Angle}
	periapsis::Union{Missing, Angle}
end

Orbit(parent::AbstractBody,
		semi_major_axis::Unitful.Length, eccentricity::Real,
		inclination::Real, ascending_node::Union{Missing, Real} = missing, periapsis::Union{Missing, Real} = missing) =
	Orbit(parent, semi_major_axis, eccentricity, to_angle(inclination), to_angle(ascending_node), to_angle(periapsis))

# needed to avoid some awkward ambiguous method resolution issues :(
Orbit(parent::AbstractBody, semi_major_axis::Unitful.Length, eccentricity::Real, inclination::Angle) =
	Orbit(parent, semi_major_axis, eccentricity, to_angle(inclination), missing, missing)

"""
    Rotation(moment_of_inertia::Union{Nothing, Real}, rotation_period::Unitful.Time, axial_tilt::Angle)

Describe the rotation of a body about its axis, with an optional `moment_of_inertia` and an `axial_tilt` relative to either its orbit or the local ecliptic depending on context.
"""
struct Rotation
	moment_of_inertia::Union{Missing, Real}
	rotation_period::Unitful.Time
	axial_tilt::Angle
end

struct Star <: AbstractBody
	name::String
	mass::Unitful.Mass
	shape::Shape
	surface_temperature::Unitful.Temperature
	absolute_magnitude::Real
	spectral_class::String

	orbit::Union{Nothing, Orbit}
	rotation::Union{Nothing, Rotation}

	function Star(name::String, mass::Unitful.Mass, shape::Shape, temp::Unitful.Temperature, amag::Real, class::String,
			orbit::Union{Nothing, Orbit} = nothing,
			rotation::Union{Nothing, Rotation} = nothing)
		b = new(name, mass, shape, temp, amag, class, orbit, rotation)

		if orbit != nothing
			add_satellite(orbit.parent, b)
		end

		return b
	end
end

"""
    Star(name::String, mass::Unitful.Mass, radius::Unitful.Length, temp::Unitful.Temperature, amag::Real, class::String)

Construct a spherical, non-rotating, solitary star with the given characteristics.
"""
Star(name::String, mass::Unitful.Mass, radius::Unitful.Length, temp::Unitful.Temperature, amag::Real, class::String) =
	Star(name, mass, Sphere(radius), temp, amag, class)

struct Body <: AbstractBody
	name::String
	mass::Unitful.Mass
	shape::Shape
	bond_albedo::Union{Missing, Real}
	geometric_albedo::Union{Missing, Real}
	orbit::Union{Nothing, Orbit}
	rotation::Union{Nothing, Rotation}

"""
    Body(name::String, mass::Unitful.Mass, shape::Shape, bond_albedo::Union{Missing, Real}, geometric_albedo::Union{Missing, Real}, orbit::Union{Nothing, Orbit}, rotation::Union{Nothing, Rotation})

Construct a new body with the given characteristics, with an optional `orbit` defining its relationship with a parent body and optional `rotation` defining the nature of its day.
"""
	function Body(name::String, mass::Unitful.Mass, shape::Shape, bond_albedo::Union{Missing, Real}, geometric_albedo::Union{Missing, Real}, orbit::Union{Nothing, Orbit}, rotation::Union{Nothing, Rotation})
		b = new(name, mass, shape, bond_albedo, geometric_albedo, orbit, rotation)

		if orbit != nothing
			add_satellite(orbit.parent, b)
		end

		return b
	end
end

"""
    Body(name::String, mass::Unitful.Mass, radius::Unitful.Length, bond_albedo::Union{Nothing, Real} = nothing)

Construct a spherical, non-rotating, solitary astronomical body with the given characteristics.
"""
Body(name::String, mass::Unitful.Mass, radius::Unitful.Length, bond_albedo::Union{Missing, Real} = missing, geometric_albedo::Union{Missing, Real} = missing) =
	Body(name, mass, Sphere(radius), bond_albedo, geometric_albedo, nothing, nothing)

Rotation(moment_of_inertia::Union{Missing, Real}, rotation_period::Unitful.Time, axial_tilt::Real) =
	Rotation(moment_of_inertia, rotation_period, to_angle(axial_tilt))

#####################################################################################

import ..SfGravity: gravity, planetary_mass, planetary_radius, orbital_period, orbital_radius, orbital_velocity, escape_velocity,
	hill_sphere, gravitational_binding_energy, roche_limit, barycentric_distance
import ..SfRelativity: relativistic_kinetic_energy
import ..SfPhysics: kinetic_energy
import ..SfMatter: density, mass
import ..SfGeometry: spherical_cap_solid_angle, volume, radius, equatorial_radius, polar_radius, area, cross_sectional_area, shape
import ..SfAstronomy: absolute_magnitude, apparent_magnitude, diffuse_sphere_q

import PhysicalConstants.CODATA2018: σ, G, k_B # σ = Stefan-Boltzmann constant, k_B Boltzmann constant

export mass, gravity, orbital_period, orbital_radius, orbital_velocity, escape_velocity, hill_sphere,
	relativistic_kinetic_energy, kinetic_energy, stellar_luminosity, stellar_irradiance, planetary_equilibrium_temperature,
	gravitational_binding_energy, roche_limit, volume, density, radius, equatorial_radius,
	area, cross_sectional_area, absolute_magnitude, apparent_magnitude, star, stellar_distance

"""
    shape(body::AbstractBody)

The shape of `body`, generally some kind of [`Ellipsoid`](@ref).
"""
shape(body::AbstractBody) = body.shape

"""
	gravity(body::AbstractBody)

Calculate the surface gravity of `body`.

# Example

```julia-repl
julia> gravity(SfSolarSystem.moon)
1.625143043265411976549419802968796018897730001404175282685973623864077610575064 m s^-2
```"""
gravity(body::AbstractBody) = gravity(body.mass, equatorial_radius(body))
planetary_mass(body::AbstractBody) = body.mass
planetary_radius(body::AbstractBody) = equatorial_radius(body)

# kepler's third
orbital_period(orbit::Orbit) = orbital_period(orbit.parent.mass, orbit.semi_major_axis)
orbital_period(body::AbstractBody) = orbital_period(body.orbit)
orbital_radius(orbit::Orbit) = orbit.semi_major_axis
orbital_radius(body::AbstractBody) = orbital_radius(body.orbit)

orbital_velocity(orbit::Orbit) = orbital_velocity(orbit.semi_major_axis, orbital_period(orbit), orbit.eccentricity)
orbital_velocity(body::AbstractBody) = orbital_velocity(body.orbit)

"""
    escape_velocity(body::AbstractBody)

Calculate the escape velocity from the surface of `body`.

If `body` is spheroidal, the escaping object is taken to be on its equator, otherwise the average radius of `body` is used.
"""
escape_velocity(body::AbstractBody) = escape_velocity(body.mass, radius(body))

"""
    escape_velocity(orbit::Orbit)

Calculate the escape velocity at the given `orbit`, based on the mass of the primary and the semimajor axis.
"""
escape_velocity(orbit::Orbit) = escape_velocity(orbit.parent.mass, orbit.semi_major_axis)

"""
    hill_sphere(body::AbstractBody)

Calculate the radius of the Hill sphere of `body`, using its current orbit and parent.
"""
hill_sphere(body::AbstractBody) = hill_sphere(body.orbit.parent.mass, body.mass, body.orbit.semi_major_axis, body.orbit.eccentricity)

"""
    hill_sphere(orbit::Orbit, mass::Unitful.Mass)

Calculate the radius of the Hill sphere of a body of `mass` in `orbit`, using the parent body described by that orbit.
"""
hill_sphere(orbit::Orbit, mass::Unitful.Mass) = hill_sphere(orbit.parent.mass, mass, orbit.semi_major_axis, orbit.eccentricity)

"""
    hill_sphere(parent::AbstractBody, satellite::AbstractBody, sma::Unitful.Length, ecc::Real = 1)

Calculate the radius of the Hill sphere of `satellite` as if it were orbiting `parent` with a semi-major axis of `sma` and eccentricity `ecc`.
"""
hill_sphere(parent::AbstractBody, satellite::AbstractBody, sma::Unitful.Length, ecc::Real = 1) = hill_sphere(mass(parent), mass(satellite), sma, ecc)

"""
    hill_sphere(parent::AbstractBody, satellite::AbstractBody, orbit::Orbit)

Calculate the radius of the Hill sphere of `satellite` as if it were orbiting `parent` in orbit `orbit`.
"""
hill_sphere(parent::AbstractBody, satellite::AbstractBody, orbit::Orbit) = hill_sphere(mass(parent), mass(satellite), orbit.semi_major_axis, orbit.eccentricity)

"""
    gravitational_binding_energy(body::AbstractBody)

Calculate the gravitational binding energy of `body`.
"""
gravitational_binding_energy(body::AbstractBody) = gravitational_binding_energy(body.mass, radius(body))

"""
    kinetic_energy(mass::Unitful.Mass, orbit::Orbit)

Compute the average kinetic energy of a body of mass `mass` in orbit `orbit`.
"""
kinetic_energy(mass::Unitful.Mass, orbit::Orbit) = kinetic_energy(mass, orbital_velocity(orbit))

"""
    kinetic_energy(body::AbstractBody)

Compute the average kinetic energy of `body` in its current orbit.
"""
kinetic_energy(body::AbstractBody) = kinetic_energy(body.mass, body.orbit)

"""
    roche_limit(primary::AbstractBody, satellite::AbstractBody)

Calculate the Roche limit of `satellite` as it approaches `primary`.
"""
roche_limit(primary::AbstractBody, satellite::AbstractBody) = roche_limit(equatorial_radius(primary), density(primary), density(satellite))

"""
    roche_limit(body::AbstractBody)

Calculate the Roche limit of `body` as it approache the object it orbits.

If `body` does not orbit anything, an error will be raised.
"""
roche_limit(body::AbstractBody) = roche_limit(equatorial_radius(body.orbit.parent), density(body.orbit.parent), density(body))

"""
    barycentric_distance(body::AbstractBody)

Find the offset from the barycentre of `body` and its parent to the centre of the parent.

If `body` does not orbit anything, an error will be raised.
"""
barycentric_distance(body::AbstractBody) = barycentric_distance(mass(body), mass(body.orbit.parent), orbital_radius(body))

"""
    volume(body::AbstractBody)

Calculate the volume of `body`.
"""
volume(body::AbstractBody) = volume(body.shape) |> u"km^3"

"""
    area(body::AbstractBody)

Calculate the surface area of `body`.

Some shapes, such as triaxial ellipsoids, may not be supported yet.
"""
area(body::AbstractBody) = area(body.shape) |> u"km^2"

"""
    cross_sectional_area(body::AbstractBody)

Calculate the approximate cross sectional area of `body`, where possible.

Results are well defined for spherical and spheroidal bodies, but may throw an error for other shapes such as triaxial ellipsoids.
"""
cross_sectional_area(body::AbstractBody) = cross_sectional_area(body.shape) |> u"km^2"

"""
    radius(body::AbstractBody)

Mean radius of `body`.
"""
radius(body::AbstractBody) = radius(body.shape)

"""
    equatorial_radius(body::AbstractBody)

Returns the equatorial radius of `body` when it is spheroidal, or an error for bodies without a well defined equatorial radius.
"""
equatorial_radius(body::AbstractBody) = equatorial_radius(body.shape)

"""
    polar_radius(body::AbstractBody)

Returns the polar radius of `body` when it is spheroidal, or an error for bodies without a well defined polar radius.
"""
polar_radius(body::AbstractBody) = polar_radius(body.shape)

"""
    radius(orbit::Orbit)

Returns the semi-major axis of `orbit`.
"""
radius(orbit::Orbit) = orbit.semi_major_axis

"""
    mass(body::AbstractBody)

Return the mass of `body`.
"""
mass(body::AbstractBody) = body.mass

star(s::Star) = s

"""
    star(body::AbstractBody)

Recuse upwards through the hierarchy of orbits to find the [`Star`](@ref) that `body` ultimately orbits. If there is no parent star, `nothing` is returned.
"""
star(body::AbstractBody) = body.orbit == nothing ? nothing : star(body.orbit.parent)

stellar_distance(::Star) = 0u"AU"

"""
	stellar_distance(body::AbstractBody)

Returns the approximate average distance of `body` to its parent star, if any. If there is no parent star, an error is raised. If `body` is a star, 0 is returned.

This returns the average radius of the top-level orbit, and so is only a loose approximation for eccentric.
"""
function stellar_distance(body::AbstractBody)
	if body.orbit == nothing
		error("Body does not ultimately orbit a star")
	elseif typeof(body.orbit.parent) == Star
		return radius(body.orbit)
	else
		return stellar_distance(body.orbit.parent)
	end
end

"""
	stellar_luminosity(r_star::Unitful.Length, t_surface::Unitful.Temperature, r_star_p = r_star)

Approximate the luminosity of a star with radius `r_star` and surface temperature `t_surface`.

Optionally, a polar radius may be given, otherwise the star is assumed to have a circular cross-section.
"""
stellar_luminosity(r_star::Unitful.Length, t_surface::Unitful.Temperature, r_star_p = r_star) = 4π * r_star * r_star_p * σ * t_surface^4 |>u"W"

"""
    stellar_luminosity(star::Star)

Approximate the luminosity of the given `star`.
"""
stellar_luminosity(star::Star) = stellar_luminosity(equatorial_radius(star), star.surface_temperature, polar_radius(star))

"""
	stellar_irradiance = function(l_stellar::Unitful.Power, r_orbit::Unitful.Length, r_body::Unitful.Length)

Compute the proportion of a star's luminosity that falls upon a body with a circular silhouette of the radius `r_body` at a distance of `r_orbit`.
"""
function stellar_irradiance(l_stellar::Unitful.Power, r_orbit::Unitful.Length, r_body::Unitful.Length)
	Ω = spherical_cap_solid_angle(r_orbit, r_body) # steradians
	return l_stellar * Ω / 4π
end

"""
	stellar_irradiance(body::AbstractBody)

Approximate the irradiance delivered to `body` by the star it ultimately orbits.

This function will recursively ascend its tree of parent bodies until a star is found. If the ultimate parent is not
a star, an error will be raised. The distance from `star` to the topmost parent is used to compute the irradiance,
as a proxy for average orbital distance.

This method will deadlock if used on a binary or multiple system at present.
"""
function stellar_irradiance(body::AbstractBody)
	s = star(body)

	if s == nothing
		error("Body does not ultimately orbit a star")
	end

	return stellar_irradiance(stellar_luminosity(s), orbital_radius(body), radius(body))
end

"""
    planetary_equilibrium_temperature(irradiance::ThermalFlux, bond_albedo::Real)

Compute the equilibrium temperature for a body with `bond_albedo` subject to stellar `irradiance`.
"""
planetary_equilibrium_temperature(irradiance::ThermalFlux, bond_albedo::Real) = ((irradiance * (1 - bond_albedo)) / 4σ)^.25 |> u"K"

"""
    planetary_equilibrium_temperature(l_stellar::Unitful.Power, r_orbit::Unitful.Length, r_body::Unitful.Length, bond_albedo::Real)

Approximate the planetary equilibrium temperature of a body of radius `r_body` and `bond_albedo`, orbiting a star of luminosity `l_stellar` at a distance of `r_orbit`.
"""
planetary_equilibrium_temperature(l_stellar::Unitful.Power, r_orbit::Unitful.Length, r_body::Unitful.Length, bond_albedo::Real) =
	planetary_equilibrium_temperature(stellar_irradiance(s_l, r_orbit, r_body) / (π * r_body^2), bond_albedo) |> u"K"

"""
    planetary_equilibrium_temperature(body::Body)

Approximate the planetary equilibrium temperature of `body`, using the irradiance of the star it ultimately orbits.

If the body does not ultimately orbit a star, an error will be raised. Binary and multiple star systems will result in deadlocks.
"""
planetary_equilibrium_temperature(body::Body) =
	planetary_equilibrium_temperature(stellar_irradiance(body) / cross_sectional_area(body), body.bond_albedo) |> u"K"

"""
    absolute_magnitude(body::Body)

Absolute magnitude of `body`, assuming it is a diffusely reflecting sphere with a known geometric albedo.
"""
absolute_magnitude(body::Body) = absolute_magnitude(radius(body), body.geometric_albedo)

"""
    apparent_magnitude(observer::AbstractBody, observed::Body, separation::Unitful.Length, phase::Angle)

Apparent magnitude of `observed` as seen from `observer` at a distance of `separation` with `phase` angle.

`observed` is assumed to be a diffusely reflecting sphere. Values will be incorrect for non-spherical bodies.
"""
apparent_magnitude(observer::AbstractBody, observed::Body, separation::Unitful.Length, phase::Angle) =
	apparent_magnitude(absolute_magnitude(observed), stellar_distance(observed), separation, stellar_distance(observer), diffuse_sphere_q(phase))

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"""
    DiscreteEvents

A Julia package for generating and simulating discrete events. It runs on Julia
`VERSION ≥ v"1.0"`. Multithreading requires `VERSION ≥ v"1.3"`.

Its Github repo is [https://github.com/pbayer/DiscreteEvents.jl].

The current stable, registered version is installed with
```julia
pkg> add DiscreteEvents
```

The development version is installed with:
```julia
pkg> add("https://github.com/pbayer/DiscreteEvents.jl")
```
"""
module DiscreteEvents

"Gives the package version."
const version = v"0.3.4"

using Unitful, Random, DataStructures, Logging, .Threads, Distributions
import Unitful: FreeUnits, Time

include("types.jl")
include("components.jl")
include("fclosure.jl")
include("schedule.jl")
include("events.jl")
include("clock.jl")
include("process.jl")
include("threads.jl")
include("timer.jl")
include("printing.jl")
include("resources.jl")
include("utils.jl")
include("macros.jl")

export  Clock, PClock, RTClock, setUnit!, 𝐶,
        Action, Timing, at, after, every, before, until,
        tau, sample_time!, fun, event!, periodic!, register!,
        incr!, run!, stop!, resume!, sync!, resetClock!,
        Prc, process!, interrupt!, delay!, wait!, now!,
        fork!, collapse!, pclock, diagnose,
        createRTClock, stopRTClock,
        PrcException,
        Resource,
        onthread, pseed!,
        @process, @event, @periodic, @delay, @wait, @run!

Random.seed!(123)
rng = MersenneTwister(2020)
𝐶.state == Undefined() ? init!(𝐶) : nothing

end # module


text:
"""
    LIF

A leaky-integrate-fire neuron described by the following differential equation

``\\frac{\\mathrm{d}v}{\\mathrm{d}t} = \\frac{R}{\\tau} I - \\lambda``

# Fields
- `voltage::VT`: membrane potential
- `current::VT`: injected (unprocessed) current
- `lastt::IT`: the last time this neuron processed a spike
- `τm::VT`: membrane time constant
- `vreset::VT`: reset voltage potential
- `R::VT`: resistive constant (typically = 1)
"""
mutable struct LIF{VT<:Real, IT<:Integer} <: AbstractCell
    # required fields
    voltage::VT
    current::VT

    # model specific fields
    lastt::IT
    τm::VT
    vreset::VT
    R::VT
end

Base.show(io::IO, ::MIME"text/plain", neuron::LIF) =
    print(io, """LIF:
                     voltage: $(neuron.voltage)
                     current: $(neuron.current)
                     τm:      $(neuron.τm)
                     vreset:  $(neuron.vreset)
                     R:       $(neuron.R)""")
Base.show(io::IO, neuron::LIF) =
    print(io, "LIF(voltage = $(neuron.voltage), current = $(neuron.current))")

"""
    LIF(τm, vreset, vth, R = 1.0)

Create a LIF neuron with zero initial voltage and empty current queue.
"""
LIF(τm::Real, vreset::Real, R::Real = 1.0) = LIF{Float32, Int}(vreset, 0, 0, τm, vreset, R)

isactive(neuron::LIF, t::Integer; dt::Real = 1.0) = (neuron.current > 0)
getvoltage(neuron::LIF) = neuron.voltage

"""
    excite!(neuron::LIF, current)
    excite!(neurons::AbstractArray{<:LIF}, current)

Excite a `neuron` with external `current`.
"""
excite!(neuron::LIF, current) = (neuron.current += current)
excite!(neurons::T, current) where T<:AbstractArray{<:LIF} = (neurons.current .+= current)

"""
    spike!(neuron::LIF, t::Integer; dt::Real = 1.0)
    spike!(neurons::AbstractArray{<:LIF}, spikes; dt::Real = 1.0)

Record a output spike from the threshold function with the `neuron` body.
Sets `neuron.lastt`.
"""
spike!(neuron::LIF, t::Integer; dt::Real = 1.0) = (t > 0) && (neuron.voltage = neuron.vreset)
function spike!(neurons::T, spikes; dt::Real = 1.0) where T<:AbstractArray{<:LIF}
    map!((x, y, s) -> (s > 0) ? y : x, neurons.voltage, neurons.voltage, neurons.vreset, spikes)
end

"""
    evaluate!(neuron::LIF, t::Integer; dt::Real = 1.0)
    (neuron::LIF)(t::Integer; dt::Real = 1.0)
    evaluate!(neurons::AbstractArray{<:LIF}, t::Integer; dt::Real = 1.0)

Evaluate the neuron model at time `t`. Return the resulting membrane potential.
"""
function evaluate!(neuron::LIF, t::Integer; dt::Real = 1.0)
    neuron.voltage = SpikingNNFunctions.Neuron.lif((t - neuron.lastt) * dt, neuron.current, neuron.voltage; R = neuron.R, tau = neuron.τm)
    neuron.lastt = t
    neuron.current = 0

    return neuron.voltage
end
(neuron::LIF)(t::Integer; dt::Real = 1.0) = evaluate!(neuron, t; dt = dt)
function evaluate!(neurons::T, t::Integer; dt::Real = 1.0) where T<:AbstractArray{<:LIF}
    SpikingNNFunctions.Neuron.lif!((t .- neurons.lastt) .* dt, neurons.current, neurons.voltage; R = neurons.R, tau = neurons.τm)
    neurons.lastt .= t
    neurons.current .= 0

    return neurons.voltage
end

"""
    reset!(neuron::LIF)
    reset!(neurons::AbstractArray{<:LIF})

Reset `neuron` by setting the membrane potential to `neuron.vreset`.
"""
function reset!(neuron::LIF)
    neuron.lastt = 0
    neuron.voltage = neuron.vreset
    neuron.current = 0
end
function reset!(neurons::T) where T<:AbstractArray{<:LIF}
    neurons.lastt .= 0
    neurons.voltage .= neurons.vreset
    neurons.current .= 0
end


text:
function mungenames!(gwas::DataFrame)
    rename!(uppercase, gwas)
    colnames = Dict(
        "SNPID" => "SNP",
        "MARKER" => "SNP",
        "MARKERNAME" => "SNP",
        "RSID" => "SNP",
        "RS_ID" => "SNP",
        "RS_NUMBER" => "SNP",
        "RS_NUMBERS" => "SNP",
        "ID" => "SNP",
        "#CHR" => "CHR",
        "#CHROM" => "CHR",
        "CHROMOSOME" => "CHR",
        "CHROM" => "CHR",
        "HG19CHR" => "CHR",
        "CHR.ICOGS" => "CHR",
        "POSITION" => "BP",
        "POS" => "BP",
        "POSITION(HG19)" => "BP",
        "BP_HG19" => "BP",
        "BASE_PAIR_LOCATION" => "BP",
        "POS_B37" => "BP",
        "POSGRCH37" => "BP",
        "POSITION.ICOGS" => "BP",
        "ALLELE1" => "A1",
        "ALLELE_1" => "A1",
        "EFFECT_ALLELE" => "A1",
        "TESTED_ALLELE" => "A1",
        "TESTEDALLELE" => "A1",
        "ALLELE1" => "A1",
        "REF" => "A1",
        "REFERENCE_ALLELE" => "A1",
        "EA" => "A1",
        "INC_ALLELE" => "A1",
        "EFFECT.META" => "A1",
        "ALLELE2" => "A2",
        "ALLELE_2" => "A2",
        "OTHER_ALLELE" => "A2",
        "OTHERALLELE" => "A2",
        "NON_EFFECT_ALLELE" => "A2",
        "NONEFFECT_ALLELE" => "A2",
        "ALT" => "A2",
        "NEA" => "A2",
        "BASELINE.META" => "A2",
        "DEC_ALLELE" => "A2",
        "ZSCORE" => "Z",
        "GC_ZSCORE" => "Z",
        "Z-SCORE" => "Z",
        "T-STATISTIC" => "Z",
        "PVAL" => "P",
        "PVALUE" => "P",
        "P-VALUE" => "P",
        "P-VAL" => "P",
        "P_DGC" => "P",
        "P_VALUE" => "P",
        "GC.PVALUE" => "P",
        "GC_PVALUE" => "P",
        "P.VALUE" => "P",
        "P.2GC" => "P",
        "ALL_INV_VAR_META_P" => "P",
        "P.META" => "P",
        "OR(A1)" => "OR",
        "OR(MINALLELE)" => "OR",
        "SE_DGC" => "SE",
        "STDERR" => "SE",
        "STANDARD_ERROR" => "SE",
        "ALL_INV_VAR_META_SEBETA" => "SE",
        "STDERRLOGOR" => "SE",
        "SE.2GC" => "SE",
        "SDE.META" => "SE",
        "EFFECT" => "BETA",
        "EFFECTS" => "BETA",
        "STDBETA" => "BETA",
        "B" => "BETA",
        "LOGOR" => "BETA",
        "LOG_ODDS" => "BETA",
        "ALL_INV_VAR_META_BETA" => "BETA",
        "BETA.META" => "BETA",
        "EAF" => "FRQ",
        "AF" => "FRQ",
        "FREQ1" => "FRQ",
        "FREQ" => "FRQ",
        "EAF_A1" => "FRQ",
        "EAF_UKB" => "FRQ",
        "MAF" => "FRQ",
        "EFFECT_ALLELE_FREQ" => "FRQ",
        "EFFECT_ALLELE_FREQUENCY" => "FRQ",
        "FREQ_TESTED_ALLELE_IN_HRS" => "FRQ",
        "EAF_EUR_UKB" => "FRQ",
        "FREQ_HAPMAP" => "FRQ",
        "EAF_HRC" => "FRQ",
        "NUMBER_OF_STUDIES" => "NSTUDY",
        "NSTUDIES" => "NSTUDY",
        "N_STUDIES" => "NSTUDY",
        "N_STUDY" => "NSTUDY",
        "NCASES" => "N_CAS",
        "NCASE" => "N_CAS",
        "N_CASE" => "N_CAS",
        "NCAS" => "N_CAS",
        "NCA" => "N_CAS",
        "N_CASES" => "N_CAS",
        "CASES_N" => "N_CAS",
        "NCONTROLS" => "N_CON",
        "NCONTROL" => "N_CON",
        "N_CONTROLS" => "N_CON",
        "CONTROLS_N" => "N_CON",
        "NCON" => "N_CON",
        "NCO" => "N_CON",
        "N_CONTROL" => "N_CON",
        "ALL_META_N" => "N",
        "N_ANALYZED" => "N",
        "IMPINFO" => "INFO",
        "INFO_UKB" => "INFO",
        "MEDIAN_INFO" => "INFO",
        "MININFO" => "INFO"
    )
    if in("SNP", [get(colnames, name, 0) for name in names(gwas)]) && "SNP" in names(gwas)
        rename!(gwas, "SNP" => "EXTRAID")
    end
    for name in names(gwas)
        haskey(colnames, name) ? rename!(gwas, name => colnames[name]) : nothing
    end
end

function mungesnpid!(gwas::DataFrame)
    if !("SNP" in names(gwas))
        gwas.SNP = string.(gwas.CHR, ":", gwas.BP)
    end
    if any(ismissing.(gwas.SNP))
        ind = ismissing.(gwas.SNP)
        gwas.SNP = string.(gwas.SNP)
        gwas.SNP[ind] = string.(gwas.CHR[ind], ":", gwas.BP[ind])
    end
end

function mungetypes!(gwas::DataFrame)
    eltype(gwas.CHR) <: AbstractString ? nothing : gwas.CHR = string.(gwas.CHR)
    if eltype(gwas.BP) != Int && eltype(gwas.BP) <: AbstractString
        gwas.BP = parse.(Int, gwas.BP)
    elseif eltype(gwas.BP) != Int && eltype(gwas.BP) <: Real
        gwas.BP = convert.(Int, gwas.BP)
    end
end

function mungealleles!(gwas::DataFrame)
    gwas.A1 = uppercase.(gwas.A1)
    gwas.A2 = uppercase.(gwas.A2)
end

function mungezscore!(gwas::DataFrame)
    if "BETA" in names(gwas) && "SE" in names(gwas)
        gwas.Z = gwas.BETA ./ gwas.SE
        filter!(x -> !isnan(x.Z), gwas)
    elseif "OR" in names(gwas) && "SE" in names(gwas)
        gwas.Z = log.(gwas.OR) ./ gwas.SE
        filter!(x -> !isnan(x.Z), gwas)
    elseif "BETA" in names(gwas)
        gwas.Z .= 0.0
        for i in 1:nrow(gwas)
            if gwas.BETA[i] >= 0
                gwas.Z[i] = sqrt(cquantile(Chisq(1), gwas.P[i]))
            else
                gwas.Z[i] = -sqrt(cquantile(Chisq(1), gwas.P[i]))
            end
        end
    elseif "OR" in names(gwas)
        gwas.BETA = log.(gwas.OR)
        gwas.Z .= 0.0
        for i in 1:nrow(gwas)
            if gwas.BETA[i] >= 0
                gwas.Z[i] = sqrt(cquantile(Chisq(1), gwas.P[i]))
            else
                gwas.Z[i] = -sqrt(cquantile(Chisq(1), gwas.P[i]))
            end
        end
    elseif "OVERALL" in names(gwas)
        gwas.Z .= 0.0
        for i in 1:nrow(gwas)
            if gwas.OVERALL[i] == "+"
                gwas.Z[i] = sqrt(cquantile(Chisq(1), gwas.P[i]))
            else
                gwas.Z[i] = -sqrt(cquantile(Chisq(1), gwas.P[i]))
            end
        end
    end
end

function mungepvalue!(gwas::DataFrame)
    filter!(x -> 0 <= x.P <= 1, gwas)
    gwas.P = clamp.(gwas.P, floatmin(0.0), 1)
end

"""
    mungesumstats!(gwas::DataFrame)
    mungesumstats!(gwas::Vector{DataFrame})

Munge `gwas` by harmonizing the names of columns, their types, and P values, among others.
"""
function mungesumstats!(gwas::Vector{DataFrame})
    for i in eachindex(gwas)
        mungenames!(gwas[i])
        ind = findall(in(["SNP", "CHR", "BP", "A1", "A2", "Z", "P", "BETA", "OR", "SE", "OVERALL"]), names(gwas[i]))
        select!(gwas[i], ind)
        mungesnpid!(gwas[i])
        dropmissing!(gwas[i])
        mungetypes!(gwas[i])
        if "A1" in names(gwas[i]) && "A2" in names(gwas[i])
            mungealleles!(gwas[i])
        end
        mungezscore!(gwas[i])
        mungepvalue!(gwas[i])
        if "A1" in names(gwas[i]) && "A2" in names(gwas[i]) && "Z" in names(gwas[i])
            select!(gwas[i], :SNP, :CHR, :BP, :A1, :A2, :Z, :P)
        elseif "A1" in names(gwas[i]) && "A2" in names(gwas[i])
            select!(gwas[i], :SNP, :CHR, :BP, :A1, :A2, :P)
        elseif "Z" in names(gwas[i])
            select!(gwas[i], :SNP, :CHR, :BP, :Z, :P)
        else
            select!(gwas[i], :SNP, :CHR, :BP, :P)
        end
    end
end

mungesumstats!(gwas::DataFrame) = mungesumstats!([gwas])

findlocus(CHR::AbstractVector, BP::AbstractVector, chr::AbstractString, range1::Real, range2::Real) =
    findall((CHR .== chr) .& (BP .>= range1) .& (BP .<= range2))

findlocus(gwas::DataFrame, chr::AbstractString, range1::Real, range2::Real) =
    findlocus(gwas.CHR, gwas.BP, chr, range1, range2)

findlocus(ref::SnpData, chr::AbstractString, range1::Real, range2::Real) =
    findlocus(ref.snp_info.chromosome, ref.snp_info.position, chr, range1, range2)

"""
    findgwasloci(gwas::DataFrame; p::Real)
    findgwasloci(gwas::Vector{DataFrame}; p::Real)

Find genome-wide significant loci for `gwas` that are separated from each
other by at least 1 Mb.

Alternatively, find genome-wide significant loci across multiple `gwas` that
are all separated by at least 1 Mb. `p` determines the genome-wide significance threshold,
which is 5e-8 by default.
"""
function findgwasloci(gwas::DataFrame; p::Real = 5e-8)
    loci = DataFrame(CHR = String[], BP = Int[], P = Float64[])
    df = gwas[findall(x -> x < p, gwas.P), [:CHR, :BP, :P]]
    while nrow(df) > 0
        ind = argmin(df.P)
        push!(loci, df[ind, :])
        ind = findlocus(df.CHR, df.BP, df.CHR[ind], df.BP[ind] - 1e6, df.BP[ind] + 1e6)
        isnothing(ind) ? nothing : deleteat!(df, ind)
    end
    return loci
end

function findgwasloci(gwas::Vector{DataFrame}; kwargs...)
    loci = Vector{DataFrame}(undef, length(gwas))
    for i in eachindex(gwas)
        loci[i] = findgwasloci(gwas[i]; kwargs...)
    end
    for i in 1:length(gwas)
        for j in (i + 1):length(gwas)
            for k in 1:nrow(loci[i])
                ind = findlocus(loci[j].CHR, loci[j].BP, loci[i].CHR[k], loci[i].BP[k] - 1e6, loci[i].BP[k] + 1e6)
                isnothing(ind) ? nothing : deleteat!(loci[j], ind)
            end
        end
    end
    return vcat(loci...)
end

function findsnps(CHR₁::AbstractVector, BP₁::AbstractVector, CHR₂::AbstractVector, BP₂::AbstractVector)
    ind = Vector{Union{Missing, Int}}(undef, length(CHR₁))
    for i in 1:length(ind)
        j = findfirst((CHR₂ .== CHR₁[i]) .& (BP₂ .== BP₁[i]))
        isnothing(j) ? ind[i] = missing : ind[i] = j
    end
    return ind
end

function findsnps(
    CHR₁::AbstractVector,
    BP₁::AbstractVector,
    A1₁::AbstractVector,
    A2₁::AbstractVector,
    CHR₂::AbstractVector,
    BP₂::AbstractVector,
    A1₂::AbstractVector,
    A2₂::AbstractVector;
    matchalleles::Bool = true
)

    if matchalleles
        ind = Matrix{Union{Missing, Int}}(undef, length(CHR₁), 2)
        for i in 1:size(ind, 1)
            j = findfirst((CHR₂ .== CHR₁[i]) .& (BP₂ .== BP₁[i]))
            if isnothing(j)
                ind[i, :] .= missing
            else
                if A1₂[j] == A1₁[i] && A2₂[j] == A2₁[i]
                    ind[i, 1], ind[i, 2] = j, missing
                elseif A2₂[j] == A1₁[i] && A1₂[j] == A2₁[i]
                    ind[i, 1], ind[i, 2] = missing, j
                else
                    ind[i, :] .= missing
                end
            end
        end
        return ind
    else
        return findsnps(CHR₁, BP₁, CHR₂, BP₂)
    end
end

function findsnps(gwas::DataFrame, ref::SnpData; matchalleles::Bool = true)
    if matchalleles
        return findsnps(gwas.CHR, gwas.BP, gwas.A1, gwas.A2,
            ref.snp_info.chromosome, ref.snp_info.position, ref.snp_info.allele1, ref.snp_info.allele2)
    else
        return findsnps(gwas.CHR, gwas.BP, ref.snp_info.chromosome, ref.snp_info.position)
    end
end

function findsnps(ref::SnpData, gwas::DataFrame; matchalleles::Bool = true)
    if matchalleles
        return findsnps(ref.snp_info.chromosome, ref.snp_info.position, ref.snp_info.allele1, ref.snp_info.allele2,
            gwas.CHR, gwas.BP, gwas.A1, gwas.A2)
    else
        return findsnps(ref.snp_info.chromosome, ref.snp_info.position, gwas.CHR, gwas.BP)
    end
end

function findsnps(gwas₁::DataFrame, gwas₂::DataFrame; matchalleles::Bool = true)
    if matchalleles
        return findsnps(gwas₁.CHR, gwas₁.BP, gwas₁.A1, gwas₁.A2, gwas₂.CHR, gwas₂.BP, gwas₂.A1, gwas₂.A2)
    else
        return findsnps(gwas₁.CHR, gwas₁.BP, gwas₂.CHR, gwas₂.BP)
    end
end

"""
    findclosestgene(chr::AbstractString, bp::Real, gencode::DataFrame; start::Bool, proteincoding::Bool)
    findclosestgene(df::DataFrame, gencode::DataFrame; start::Bool, proteincoding::Bool)

Find the closest gene(s) to a genomic coordinate or a list of genomic coordinates using `gencode`.

Optionally, the closest gene can be defined from the gene start site using `start`,
and only protein coding genes can be considered using `proteincoding`.
The default `start` and `proteincoding` are `false`.
"""
function findclosestgene(
    chr::AbstractString,
    bp::Real,
    gencode::DataFrame;
    start::Bool = false,
    proteincoding::Bool = false
)

    if proteincoding
        df = filter(x -> (x.seqnames == chr) && (x.feature == "gene") && (x.gene_type == "protein_coding"), gencode)
    else
        df = filter(x -> (x.seqnames == chr) && (x.feature == "gene"), gencode)
    end
    if start
        df.dist₁ .= 1
        for i in 1:nrow(df)
            df.strand[i] == "+" ? df.dist₁[i] = abs(df.start[i] - bp) : df.dist₁[i] = abs(df.end[i] - bp)
        end
        ind₁ = argmin(df.dist₁)
        return df.gene_name[ind₁], df.dist₁[ind₁]
    else
        df.dist₁ = abs.(df.start .- bp)
        df.dist₂ = abs.(df.end .- bp)
        ind₁ = argmin(df.dist₁)
        ind₂ = argmin(df.dist₂)
        if df.dist₁[ind₁] < df.dist₂[ind₂]
            return df.gene_name[ind₁], df.dist₁[ind₁]
        elseif df.dist₁[ind₁] > df.dist₂[ind₂]
            return df.gene_name[ind₂], df.dist₂[ind₂]
        elseif ind₁ == ind₂
            return df.gene_name[ind₁], df.dist₁[ind₁]
        elseif df.strand[ind₁] == "+" && df.strand[ind₂] == "+"
            return df.gene_name[ind₁], df.dist₁[ind₁]
        elseif df.strand[ind₁] == "-" && df.strand[ind₂] == "-"
            return df.gene_name[ind₂], df.dist₁[ind₂]
        elseif df.gene_type[ind₁] == "protein_coding" && df.gene_type[ind₂] != "protein_coding"
            return df.gene_name[ind₁], df.dist₁[ind₁]
        elseif df.gene_type[ind₁] != "protein_coding" && df.gene_type[ind₂] == "protein_coding"
            return df.gene_name[ind₂], df.dist₁[ind₂]
        else
            return df.gene_name[ind₁], df.dist₁[ind₁]
        end
    end
end

function findclosestgene(CHR::AbstractVector, BP::AbstractVector, gencode::DataFrame; kwargs...)
    gencodeₛ = select(gencode, :seqnames, :start, :end, :strand, :gene_type, :feature, :gene_name)
    filter!(x -> x.feature == "gene", gencodeₛ)
    storage = DataFrame(CHR = String[], BP = Int[], gene = String[], distance = Int[])
    for i in eachindex(CHR)
        gene, dist = findclosestgene(CHR[i], BP[i], gencodeₛ; kwargs...)
        push!(storage, [CHR[i], BP[i], gene, dist])
    end
    return storage
end

findclosestgene(df::DataFrame, gencode::DataFrame; kwargs...) = findclosestgene(df.CHR, df.BP, gencode; kwargs...)

function findclosestgenes(
    chr::AbstractString,
    bp::Real,
    gencode::DataFrame;
    start::Bool = false,
    proteincoding::Bool = false,
    n::Real = 5
)

    if proteincoding
        df = filter(x -> (x.seqnames == chr) && (x.feature == "gene") && (x.gene_type == "protein_coding"), gencode)
    else
        df = filter(x -> (x.seqnames == chr) && (x.feature == "gene"), gencode)
    end
    if start
        df.dist₁ .= 1
        for i in 1:nrow(df)
            df.strand[i] == "+" ? df.dist₁[i] = abs(df.start[i] - bp) : df.dist₁[i] = abs(df.end[i] - bp)
        end
        sort!(df, :dist₁)
        select!(df, :gene_name, :dist₁)
        rename!(df, :dist₁ => :distance)
        return first(df, n)
    else
        df.dist₁ = abs.(df.start .- bp)
        df.dist₂ = abs.(df.end .- bp)
        df.dist = min.(df.dist₁, df.dist₂)
        sort!(df, :dist)
        select!(df, :gene_name, :dist)
        rename!(df, :dist => :distance)
        return first(df, n)
    end
end

function getsnpinfo(snp::AbstractString, SNP::AbstractVector, CHR::AbstractVector, BP::AbstractVector)
    ind = findfirst(isequal(snp), SNP)
    isnothing(ind) ? nothing : (CHR[ind], BP[ind])
end

getsnpinfo(snp::AbstractString, df::DataFrame) = getsnpinfo(snp, df.SNP, df.CHR, df.BP)

getsnpinfo(snp::AbstractString, ref::SnpData) =
    getsnpinfo(snp, ref.snp_info.snpid, ref.snp_info.chromosome, ref.snp_info.position)

function getsnpinfo(chr::AbstractString, bp::Real, SNP::AbstractVector, CHR::AbstractVector, BP::AbstractVector)
    ind = findfirst((CHR .== chr) .& (BP .== bp))
    isnothing(ind) ? nothing : SNP[ind]
end

getsnpinfo(chr::AbstractString, bp::Real, df::DataFrame) = getsnpinfo(chr, bp, df.SNP, df.CHR, df.BP)

getsnpinfo(chr::AbstractString, bp::Real, ref::SnpData) =
    getsnpinfo(chr, bp, ref.snp_info.snpid, ref.snp_info.chromosome, ref.snp_info.position)

function flipalleles(gwas::Vector{DataFrame}; outer = false)
    if outer
        df = outerjoin(gwas..., on = [:CHR, :BP], makeunique = true)
        for i in 1:(length(gwas) - 1)
            for j in 1:nrow(df)
                if ismissing(df.A1[j]) || ismissing(df.A2[j]) || ismissing(df[j, "A1_$i"]) || ismissing(df[j, "A2_$i"])
                    continue
                elseif df.A1[j] == df[j, "A1_$i"] && df.A2[j] == df[j, "A2_$i"]
                    continue
                elseif df.A1[j] == df[j, "A2_$i"] && df.A2[j] == df[j, "A1_$i"]
                    df[j, "Z_$i"] = -df[j, "Z_$i"]
                else
                    df[j, "Z_$i"] = NaN
                end
            end
        end
    else
        df = innerjoin(gwas..., on = [:CHR, :BP], makeunique = true)
        for i in 1:(length(gwas) - 1)
            for j in 1:nrow(df)
                if df.A1[j] == df[j, "A1_$i"] && df.A2[j] == df[j, "A2_$i"]
                    continue
                elseif df.A1[j] == df[j, "A2_$i"] && df.A2[j] == df[j, "A1_$i"]
                    df[j, "Z_$i"] = -df[j, "Z_$i"]
                else
                    df[j, "Z_$i"] = NaN
                end
            end
        end
    end
    return df
end

flipalleles(gwas₁::DataFrame, gwas₂::DataFrame; kwargs...) = flipalleles([gwas₁, gwas₂]; kwargs...)

GRCh38_highld = [
    (chr = "1", range1 = 47761740, range2 = 51761740),
    (chr = "1", range1 = 125169943, range2 = 125170022),
    (chr = "1", range1 = 144106678, range2 = 144106709),
    (chr = "1", range1 = 181955019, range2 = 181955047),
    (chr = "2", range1 = 85919365, range2 = 100517106),
    (chr = "2", range1 = 135275091, range2 = 135275210),
    (chr = "2", range1 = 182427027, range2 = 189427029),
    (chr = "3", range1 = 47483505, range2 = 49987563),
    (chr = "3", range1 = 83368158, range2 = 86868160),
    (chr = "5", range1 = 44464140, range2 = 51168409),
    (chr = "5", range1 = 129636407, range2 = 132636409),
    (chr = "6", range1 = 25391792, range2 = 33424245),
    (chr = "6", range1 = 57788603, range2 = 58453888),
    (chr = "6", range1 = 61109122, range2 = 61424451),
    (chr = "6", range1 = 139637169, range2 = 142137170),
    (chr = "7", range1 = 54964812, range2 = 66897578),
    (chr = "8", range1 = 8105067, range2 = 12105082),
    (chr = "8", range1 = 43025699, range2 = 48924888),
    (chr = "8", range1 = 110918594, range2 = 113918595),
    (chr = "9", range1 = 64198500, range2 = 64200392),
    (chr = "10", range1 = 36671065, range2 = 43184546),
    (chr = "11", range1 = 88127183, range2 = 91127184),
    (chr = "12", range1 = 32955798, range2 = 41319931),
    (chr = "20", range1 = 33948532, range2 = 36438183),
]

GRCh37_highld = [
    (chr = "1", range1 = 48000000, range2 = 52000000),
    (chr = "2", range1 = 86000000, range2 = 100500000),
    (chr = "2", range1 = 134500000, range2 = 138000000),
    (chr = "2", range1 = 183000000, range2 = 190000000),
    (chr = "3", range1 = 47500000, range2 = 50000000),
    (chr = "3", range1 = 83500000, range2 = 87000000),
    (chr = "3", range1 = 89000000, range2 = 97500000),
    (chr = "5", range1 = 44500000, range2 = 50500000),
    (chr = "5", range1 = 98000000, range2 = 100500000),
    (chr = "5", range1 = 129000000, range2 = 132000000),
    (chr = "5", range1 = 135500000, range2 = 138500000),
    (chr = "6", range1 = 24000000, range2 = 36000000),
    (chr = "6", range1 = 57000000, range2 = 64000000),
    (chr = "6", range1 = 140000000, range2 = 142500000),
    (chr = "7", range1 = 55000000, range2 = 66000000),
    (chr = "8", range1 = 7000000, range2 = 13000000),
    (chr = "8", range1 = 43000000, range2 = 50000000),
    (chr = "8", range1 = 112000000, range2 = 115000000),
    (chr = "10", range1 = 37000000, range2 = 43000000),
    (chr = "11", range1 = 46000000, range2 = 57000000),
    (chr = "11", range1 = 87500000, range2 = 90500000),
    (chr = "12", range1 = 33000000, range2 = 40000000),
    (chr = "12", range1 = 109500000, range2 = 112000000),
    (chr = "20", range1 = 32000000, range2 = 34500000),
    (chr = "17", range1 = 42900000, range2 = 45100000)
]
# https://genome.sph.umich.edu/wiki/Regions_of_high_linkage_disequilibrium_(LD)

GRCh38_totlength = Dict(
    "1" => 248956422,
    "2" => 242193529,
    "3" => 198295559,
    "4" => 190214555,
    "5" => 181538259,
    "6" => 170805979,
    "7" => 159345973,
    "8" => 145138636,
    "9" => 138394717,
    "10" => 133797422,
    "11" => 135086622,
    "12" => 133275309,
    "13" => 114364328,
    "14" => 107043718,
    "15" => 101991189,
    "16" => 90338345,
    "17" => 83257441,
    "18" => 80373285,
    "19" => 58617616,
    "20" => 64444167,
    "21" => 46709983,
    "22" => 50818468,
    "X" => 156040895,
    "Y" => 57227415,
)

GRCh37_totlength = Dict(
    "1" => 249250621,
    "2" => 243199373,
    "3" => 198022430,
    "4" => 191154276,
    "5" => 180915260,
    "6" => 171115067,
    "7" => 159138663,
    "8" => 146364022,
    "9" => 141213431,
    "10" => 135534747,
    "11" => 135006516,
    "12" => 133851895,
    "13" => 115169878,
    "14" => 107349540,
    "15" => 102531392,
    "16" => 90354753,
    "17" => 81195210,
    "18" => 78077248,
    "19" => 59128983,
    "20" => 63025520,
    "21" => 48129895,
    "22" => 51304566,
    "X" => 155270560,
    "Y" => 59373566,
)
# https://www.ncbi.nlm.nih.gov/grc/human/data?asm=GRCh37

gwas = Dict(
    "scz" => (url = "https://figshare.com/ndownloader/files/28169757",
        PMID = "",
        title = "Schizophrenia (PGC3)",
        file = "PGC3_SCZ_wave3_public.v2.tsv.gz"),
    "bd" => (url = "https://figshare.com/ndownloader/files/26603681",
        PMID = "34002096",
        title = "Bipolar (Mullins et al. 2021)",
        file = "pgc-bip2021-all.vcf.tsv.gz"),
    "asd" => (url = "https://figshare.com/ndownloader/files/28169292",
        PMID = "30804558",
        title = "Autism (Grove et al. 2019)",
        file = "iPSYCH-PGC_ASD_Nov2017.gz"),
    "sczvsbd" => (url = "https://figshare.com/ndownloader/files/28169358",
        PMID = "29906448",
        title = "SCZ vs BD (Ruderfer et al. 2018)",
        file = "SCZvsBD.sumstats.gz"),
    "adhd" => (url = "https://figshare.com/ndownloader/files/28169253",
        PMID = "30478444",
        title = "ADHD (Demontis et al. 2019)",
        file = "daner_adhd_meta_filtered_NA_iPSYCH23_PGC11_sigPCs_woSEX_2ell6sd_EUR_Neff_70.meta.gz"),
    "anorexia" => (url = "https://figshare.com/ndownloader/files/28169271",
        PMID = "31308545",
        title = "Anorexia (Watson et al. 2019)",
        file = "pgcAN2.2019-07.vcf.tsv.gz"),
    "amd" => (url = "http://csg.sph.umich.edu/abecasis/public/amd2015/Fritsche_2015_AdvancedAMD.txt.gz",
        PMID = "26691988",
        title = "AMD (Fritsche et al. 2016)",
        file = "Fritsche_2015_AdvancedAMD.txt.gz"),
    "menarche" => (url = "https://www.reprogen.org/Menarche_1KG_NatGen2017_WebsiteUpload.zip",
        PMID = "28436984",
        title = "Age at menarche (Day et al. 2017)",
        file = "Menarche_1KG_NatGen2017_WebsiteUpload.txt"),
    "menopause" => (url = "https://www.reprogen.org/reprogen_ANM_201K_170621.txt.gz",
        PMID = "34349265",
        title = "Age at menoapuse (Ruth et al. 2021)",
        file = "reprogen_ANM_201K_170621.txt.gz"),
    "alz" => (url = "https://ctg.cncr.nl/documents/p1651/AD_sumstats_Jansenetal_2019sept.txt.gz",
        PMID = "30617256",
        title = "Alzheimer disease (Jansen et al. 2019)",
        file = "AD_sumstats_Jansenetal_2019sept.txt.gz"),
    "cross" => (url = "https://figshare.com/ndownloader/files/28169382",
        PMID = "31835028",
        title = "Cross disorder (PGC 2019)",
        file = "pgc_cdg2_meta_no23andMe_oct2019_v2.txt.daner.txt.gz"),
    "mdd" => (url = "https://datashare.ed.ac.uk/bitstream/handle/10283/3203/PGC_UKB_depression_genome-wide.txt",
        PMID = "30718901",
        title = "Major depression (Howard et al. 2019)",
        file = "PGC_UKB_depression_genome-wide.txt"),
    "lupus" => (url = "http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST003001-GCST004000/GCST003156/bentham_2015_26502338_sle_efo0002690_1_gwas.sumstats.tsv.gz",
        PMID = "26502338",
        title = "Lupus (Bentham et al. 2015)",
        file = "bentham_2015_26502338_sle_efo0002690_1_gwas.sumstats.tsv.gz"),
    "height" => (url = "https://portals.broadinstitute.org/collaboration/giant/images/6/63/Meta-analysis_Wood_et_al%2BUKBiobank_2018.txt.gz",
        PMID = "30124842",
        title = "Height (Yengo et al. 2018)",
        file = "Meta-analysis_Wood_et_al%2BUKBiobank_2018.txt.gz"),
    "weight" => (url = "https://portals.broadinstitute.org/collaboration/giant/images/c/c8/Meta-analysis_Locke_et_al%2BUKBiobank_2018_UPDATED.txt.gz",
        PMID = "30124842",
        title = "Weight (Yengo et al. 2018)",
        file = "Meta-analysis_Locke_et_al%2BUKBiobank_2018_UPDATED.txt.gz"),
    "left-hand" => (url = "https://evansgroup.di.uq.edu.au/GWAS_RESULTS/HANDEDNESS/LeftHandedness_MetaAnalysis_UKBB_IHC.txt.gz",
        PMID = "32989287",
        title = "Left-handedness (Cuellar-Partida et al. 2020)",
        file = "LeftHandedness_MetaAnalysis_UKBB_IHC.txt.gz"),
    "tsh" => (url = "http://csg.sph.umich.edu/willer/public/TSH2020/2020_Zhou_et_al_TSH_meta_HUNT-MGI-ThyroidOmics.txt",
        PMID = "32769997",
        title = "TSH level (Zhou et al. 2020)",
        file = "2020_Zhou_et_al_TSH_meta_HUNT-MGI-ThyroidOmics.txt"),
    "ra" => (url = "http://plaza.umin.ac.jp/~yokada/datasource/files/GWASMetaResults/RA_GWASmeta_European_v2.txt.gz",
        PMID = "24390342",
        title = "Rheumatoid arthritis (Okada et al. 2014)",
        file = "RA_GWASmeta_European_v2.txt.gz"),
    "als" => (url = "https://surfdrive.surf.nl/files/index.php/s/E5RetKw10hC3jXy/download?path=%2F&files=ALS_GWAS_RVB_SMR_sumstats.tar.gz",
        PMID = "34873335",
        title = "ALS (van Rheenen et al. 2021)",
        file = "ALS_sumstats_EUR_only.txt.gz"),
    "parkinson" => (url = "https://drive.google.com/file/d/1FZ9UL99LAqyWnyNBxxlx6qOUlfAnublN",
        PMID = "31701892",
        title = "Parkinson's disease (Nalls et al. 2019)",
        file = "nallsEtAl2019_excluding23andMe_allVariants.tab"),
    "cp" => (url = "https://thessgac.com/papers/",
        PMID = "30038396",
        title = "Cognitive performance (Lee et al. 2018)",
        file = "GWAS_CP_all.txt"),
    "ea" => (url = "https://thessgac.com/papers/",
        PMID = "35361970",
        title = "Educational attainment (Okbay et al. 2022)",
        file = "EA4_additive_excl_23andMe.txt.gz"),
    "birth" => (url = "https://thessgac.com/papers/",
        PMID = "27798627",
        title = "Age at first birth (Barban et al. 2016)",
        file = "AgeFirstBirth_Pooled.txt"),
    "children" => (url = "https://thessgac.com/papers/",
        PMID = "27798627",
        title = "Number of children (Barban et al. 2016)",
        file = "NumberChildrenEverBorn_Pooled.txt"),
    "risk" => (url = "https://thessgac.com/papers/",
        PMID = "30643258",
        title = "Risk tolerance (Linner et al. 2019)",
        file = "RISK_GWAS_MA_UKB+replication.txt"),
    "neuroticism" => (url = "https://ctg.cncr.nl/documents/p1651/sumstats_neuroticism_ctg_format.txt.gz",
        PMID = "29942085",
        title = "Neuroticism (Nagel et al. 2018)",
        file = "sumstats_neuroticism_ctg_format.txt.gz"),
    "cad" => (url = "http://www.cardiogramplusc4d.org/data-downloads/",
        PMID = "26343387",
        title = "CAD (Nikpay et al. 2015)",
        file = "cad.add.160614.website.txt"),
    "afib" => (url = "https://personal.broadinstitute.org/ryank/AF_HRC_GWAS_ALLv11.zip",
        PMID = "29892015",
        title = "Atrial fibrillation (Roselli et al. 2018)",
        file = "AF_HRC_GWAS_ALLv11.txt"),
    "stroke" => (url = "https://www.kp4cd.org/datasets/stroke",
        PMID = "29531354",
        title = "Stroke (Malik et al. 2018)",
        file = "metastroke.all.chr.bp"),
    "ibd" => (url = "https://www.ibdgenetics.org/downloads.html",
        PMID = "26192919",
        title = "IBD (Liu et al. 2015)",
        file = "EUR.IBD.gwas_info03_filtered.assoc"),
    "ckd" => (url = "http://ckdgen.imbi.uni-freiburg.de/files/Wuttke2019/CKD_overall_EA_JW_20180223_nstud23.dbgap.txt.gz",
        PMID = "31152163",
        title = "CKD (Wuttke et al. 2019)",
        file = "CKD_overall_EA_JW_20180223_nstud23.dbgap.txt.gz"),
    "t2d" => (url = "https://diagram-consortium.org/downloads.html",
        PMID = "30297969",
        title = "Type II diabetes (Mahajan et al. 2018)",
        file = "Mahajan.NatGenet2018b.T2D.European.txt"),
    "smoking" => (url = "https://conservancy.umn.edu/bitstream/handle/11299/201564/CigarettesPerDay.txt.gz",
        PMID = "30643251",
        title = "Cigarettes per day (Liu et al. 2019)",
        file = "CigarettesPerDay.txt.gz"),
    "alcohol" => (url = "https://conservancy.umn.edu/bitstream/handle/11299/201564/DrinksPerWeek.txt.gz",
        PMID = "30643251",
        title = "Drinks per week (Liu et al. 2019)",
        file = "DrinksPerWeek.txt.gz"),
    "intelligence" => (url = "https://ctg.cncr.nl/documents/p1651/SavageJansen_IntMeta_sumstats.zip",
        PMID = "29942086",
        title = "Intelligence (Savage et al. 2018)",
        file = "SavageJansen_2018_intelligence_metaanalysis.txt"),
    "insomnia" => (url = "https://ctg.cncr.nl/documents/p1651/Insomnia_sumstats_Jansenetal.txt.gz",
        PMID = "30804565",
        title = "Insomnia (Jansen et al. 2019)",
        file = "Insomnia_sumstats_Jansenetal.txt.gz"),
    "icv" => (url = "https://ctg.cncr.nl/documents/p1651/meta_analysis_BV_Jansenetal_2020.sumstats.txt.gz",
        PMID = "33154357",
        title = "Intracranial volume (Jansen et al. 2020)",
        file = "meta_analysis_BV_Jansenetal_2020.sumstats.txt.gz"),
    "covid" => (url = "https://storage.googleapis.com/covid19-hg-public/20210415/results/20210607/COVID19_HGI_A2_ALL_leave_23andme_20210607.b37.txt.gz",
        PMID = "34237774",
        title = "COVID-19 (2021)",
        file = "COVID19_HGI_A2_ALL_leave_23andme_20210607.b37.txt.gz"),
    "acne" => (url = "http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST90091001-GCST90092000/GCST90092000/GCST90092000_buildGRCh37.tsv",
        PMID = "35132056",
        title = "Acne (Mitchell et al. 2022)",
        file = "GCST90092000_buildGRCh37.tsv.gz"),
    "breast" => (url = "https://bcac.ccge.medschl.cam.ac.uk/bcacdata/oncoarray/oncoarray-and-combined-summary-result/gwas-summary-associations-breast-cancer-risk-2020/",
        PMID = "32424353",
        title = "Breast cancer (Zhang et al. 2020)",
        file = "icogs_onco_gwas_meta_overall_breast_cancer_summary_level_statistics.txt"),
    "prostate" => (url = "https://bcac.ccge.medschl.cam.ac.uk/bcacdata/oncoarray/oncoarray-and-combined-summary-result/gwas-summary-associations-breast-cancer-risk-2020/",
        PMID = "23065704",
        title = "Prostate cancer (Schumacher et al. 2018)",
        file = "meta_v3_onco_euro_overall_ChrAll_1_release.txt")
)

function downloadgwas(path::AbstractString; pheno::Union{AbstractVector, AbstractString, Nothing} = nothing)
    noteasy = ["als", "menarche", "parkinson", "intelligence", "t2d", "ibd", "stroke", "afib",
        "cad", "cp", "ea", "birth", "children", "risk", "acne"]
    if isnothing(pheno)
        for key in keys(gwas)
            key in noteasy ? continue : nothing
            @info "Downloading summary statistics for $(key)."
            url = gwas[key].url
            file = gwas[key].file
            isdir(path) || mkdir(path)
            isfile("$(path)/$(file)") || download(url, "$(path)/$(file)")
        end
    elseif isa(pheno, AbstractVector)
        for key in pheno
            if !(key in keys(gwas)) || (key in noteasy)
                @error "Cannot download summary statistics for $(key)."
            else
                @info "Downloading summary statistics for $(key)."
                url = gwas[key].url
                file = gwas[key].file
                isdir(path) || mkdir(path)
                isfile("$(path)/$(file)") || download(url, "$(path)/$(file)")
            end
        end
    else
        if !(pheno in keys(gwas)) || (pheno in noteasy)
            @error "Cannot download summary statistics for $(pheno)."
        else
            @info "Downloading summary statistics for $(pheno)."
            url = gwas[pheno].url
            file = gwas[pheno].file
            isdir(path) || mkdir(path)
            isfile("$(path)/$(file)") || download(url, "$(path)/$(file)")
        end
    end
end


text:
module GaussianMixtureTest

using Distributions, StatsBase, StatsFuns, Yeppp, LinearAlgebra

#import Yeppp: add!, exp!, log!
import StatsBase: RealArray, RealVector, RealArray, IntegerArray, IntegerVector, IntegerMatrix, IntUnitRange

export gmm, gmmrepeat, asymptoticdistribution, kstest

include("estimate.jl")
include("arithmetic.jl")
include("test.jl")
end # module


text:
module AdversarialAttacks

using Zygote
using Adapt
using LinearAlgebra: norm

export fgsm, fgsm!, pgd, pgd!

include("utils.jl")
include("attacks.jl")

end


text:
import Printf: @sprintf

"""
Change of the derivative.
"""
struct Event
    """Position (to sort the queues)"""
    x::Float64

    """Delta slope"""
    slope::Float64
end

"""
    intercept(e)

Intercept (also called *offset*); implicitly encoded in the [`Event`](@ref).
"""
@inline intercept(e::Event)::Float64 = -e.x * e.slope

Base.isless(e1::Event, e2::Event) = isless(e1.x, e2.x)

Base.zero(::Type{Event}) = Event(NaN, NaN, NaN)

Base.isnan(e::Event) = isnan(e.x)

Base.show(io::IO, e::Event) =
    print(io,
          if isnan(e) "Event(NaN)" else
              @sprintf("% .2f -> % .2ft + % .3f",
                       e.x, e.slope, intercept(e))
          end)

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.

export addlazycallback, addcutcallback, addheuristiccallback, addinfocallback

abstract JuMPCallback
type LazyCallback <: JuMPCallback
    f::Function
    fractional::Bool
end
type CutCallback <: JuMPCallback
    f::Function
end
type HeuristicCallback <: JuMPCallback
    f::Function
end
type InfoCallback <: JuMPCallback
    f::Function
    when::Symbol
end

type CallbackAbort <: Exception end
export CallbackAbort

Base.copy{T<:JuMPCallback}(c::T) = T(copy(c))
Base.copy(c::LazyCallback) = LazyCallback(copy(c.f), c.fractional)

function addlazycallback(m::Model, f::Function; fractional::Bool=false)
    m.internalModelLoaded = false
    push!(m.callbacks, LazyCallback(f,fractional))
end
function addcutcallback(m::Model, f::Function)
    m.internalModelLoaded = false
    push!(m.callbacks, CutCallback(f))
end
function addheuristiccallback(m::Model, f::Function)
    m.internalModelLoaded = false
    push!(m.callbacks, HeuristicCallback(f))
end

function _unspecifiedstate()
    warn("""Info Callbacks without a when clause will currently default
         to fire only in the :Intermediate state to preserve its behavior
         with previous versions of JuMP.

         This behavior will be deprecated in subsequent versions of JuMP, so
         please rewrite all invocations of addinfocallbacks(m, f) to
         addinfocallbacks(m, f, when = :Intermediate) instead.
         """)
    :Intermediate
end

function addinfocallback(m::Model, f::Function; when::Symbol = _unspecifiedstate())
    m.internalModelLoaded = false
    push!(m.callbacks, InfoCallback(f, when))
end

function lazycallback(d::MathProgBase.MathProgCallbackData, m::Model, cbs::Vector{LazyCallback})
    state = MathProgBase.cbgetstate(d)
    @assert state == :MIPSol || state == :MIPNode
    anyfrac = mapreduce(|, cbs) do cb
        cb.fractional
    end
    if state == :MIPSol
        MathProgBase.cbgetmipsolution(d,m.colVal)
    elseif anyfrac
        MathProgBase.cbgetlpsolution(d,m.colVal)
    else
        return :Continue
    end
    try
        for cb in cbs
            if state == :MIPSol || cb.fractional
                cb.f(d)
            end
        end
    catch y
        if isa(y, CallbackAbort)
            return :Exit
        else
            rethrow(y)
        end
    end
    :Continue
end

function attach_callbacks(m::Model, cbs::Vector{LazyCallback})
    cb = d -> lazycallback(d,m,cbs)
    if applicable(MathProgBase.setlazycallback!, m.internalModel, cb)
        MathProgBase.setlazycallback!(m.internalModel, cb)
    else
        error("Solver does not support lazy callbacks")
    end
end

function cutcallback(d::MathProgBase.MathProgCallbackData, m::Model, cbs::Vector{CutCallback})
    state = MathProgBase.cbgetstate(d)
    @assert state == :MIPNode
    MathProgBase.cbgetlpsolution(d,m.colVal)
    try
        for cb in cbs
            cb.f(d)
        end
    catch y
        if isa(y, CallbackAbort)
            return :Exit
        else
            rethrow(y)
        end
    end
    :Continue
end

function attach_callbacks(m::Model, cbs::Vector{CutCallback})
    cb = d -> cutcallback(d,m,cbs)
    if applicable(MathProgBase.setcutcallback!, m.internalModel, cb)
        MathProgBase.setcutcallback!(m.internalModel, cb)
    else
        error("Solver does not support cut callbacks")
    end
end

function heurcallback(d::MathProgBase.MathProgCallbackData, m::Model, cbs::Vector{HeuristicCallback})
    state = MathProgBase.cbgetstate(d)
    @assert state == :MIPNode
    MathProgBase.cbgetlpsolution(d,m.colVal)
    try
        for cb in cbs
            cb.f(d)
        end
    catch y
        if isa(y, CallbackAbort)
            return :Exit
        else
            rethrow(y)
        end
    end
    :Continue
end

function attach_callbacks(m::Model, cbs::Vector{HeuristicCallback})
    cb = d -> heurcallback(d,m,cbs)
    if applicable(MathProgBase.setheuristiccallback!, m.internalModel, cb)
        MathProgBase.setheuristiccallback!(m.internalModel, cb)
    else
        error("Solver does not support heuristic callbacks")
    end
end

function infocallback(d::MathProgBase.MathProgCallbackData, m::Model, cbs::Vector{InfoCallback})
    state = MathProgBase.cbgetstate(d)
    if state == :MIPSol
        MathProgBase.cbgetmipsolution(d,m.colVal)
    elseif state == :MIPNode
        MathProgBase.cbgetlpsolution(d,m.colVal)
    end
    try
        for cb in cbs
            if cb.when == state
                cb.f(d)
            end
        end
    catch y
        if isa(y, CallbackAbort)
            return :Exit
        else
            rethrow(y)
        end
    end
    :Continue
end

function attach_callbacks(m::Model, cbs::Vector{InfoCallback})
    cb = d -> infocallback(d,m,cbs)
    if applicable(MathProgBase.setinfocallback!, m.internalModel, cb)
        MathProgBase.setinfocallback!(m.internalModel, cb)
    else
        error("Solver does not support info callbacks")
    end
end

function registercallbacks(m::Model)
    isempty(m.callbacks) && return # might as well avoid allocating the indexedVector

    cbtypes = unique(map(typeof, m.callbacks))
    for typ in cbtypes
        cbs::Vector{typ} = filter(m.callbacks) do cb
            isa(cb, typ)
        end
        attach_callbacks(m, cbs)
    end

    # prepare storage for callbacks
    m.indexedVector = IndexedVector(Float64, m.numCols)
end

# TODO: Should this be somewhere else?
const sensemap = Dict(:(<=) => '<', :(==) => '=', :(>=) => '>')

## Lazy constraints
export @lazyconstraint

macro lazyconstraint(args...)
    cbdata = esc(args[1])
    x = args[2]
    extra = vcat(args[3:end]...)
    # separate out keyword arguments
    kwargs = filter(ex->isexpr(ex,:kw), extra) # filtering expressions corresponding to kw args specs
    extra = filter(ex->!isexpr(ex,:kw), extra) # others

    localcut_val = false # by default, the lazy constraint is global
    for ex in kwargs
        kwarg = ex.args[1]
        if kwarg == :localcut
            localcut_val = esc(ex.args[2])   # excepted if otherwise specified ...
        else error("in @lazyconstraint($(join(args,','))), invalid keyword argument: $(kwarg)")
        end
    end

    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    if isexpr(x, :call) && length(x.args) == 3 # simple comparison
        lhs = :($(x.args[2]) - $(x.args[3])) # move everything to the lhs
        newaff, parsecode = parseExprToplevel(lhs, :aff)
        sense, vectorized = _canonicalize_sense(x.args[1])
        code = quote
            aff = zero(AffExpr)
            $parsecode
            constr = constructconstraint!($newaff, $(quot(sense)))
        end
        if vectorized
            return quote
                $code
                for con in constr
                    addlazyconstraint($cbdata, con, localcut=$localcut_val)
                end
            end
        else
            return quote
                $code
                addlazyconstraint($cbdata, constr, localcut=$localcut_val)
            end
        end
    else
        error("Syntax error in @lazyconstraint, expected one-sided comparison.")
    end
end

function addlazyconstraint(cbdata::MathProgBase.MathProgCallbackData, constr::LinearConstraint; localcut::Bool=false)
    if length(constr.terms.vars) == 0
        MathProgBase.cbaddlazy!(cbdata, Cint[], Float64[], sensemap[sense(constr)], rhs(constr))
        return
    end
    assert_isfinite(constr.terms)
    m::Model = constr.terms.vars[1].m
    indices, coeffs = merge_duplicates(Cint, constr.terms, m.indexedVector, m)
    if localcut
        MathProgBase.cbaddlazylocal!(cbdata, indices, coeffs, sensemap[sense(constr)], rhs(constr))
    else
        MathProgBase.cbaddlazy!(cbdata, indices, coeffs, sensemap[sense(constr)], rhs(constr))
    end
end

addlazyconstraint(cbdata::MathProgBase.MathProgCallbackData,constr::QuadConstraint; localcut::Bool=false) = error("Quadratic lazy constraints are not supported.")

## User cuts
export addusercut
export @usercut

macro usercut(args...)
    cbdata = esc(args[1])
    x = args[2]
    extra = vcat(args[3:end]...)
    # separate out keyword arguments
    kwargs = filter(ex->isexpr(ex,:kw), extra) # filtering expressions corresponding to kw args specs
    extra = filter(ex->!isexpr(ex,:kw), extra) # others

    localcut_val = false # by default, the user cut is global
    for ex in kwargs
        kwarg = ex.args[1]
        if kwarg == :localcut
            localcut_val = esc(ex.args[2])   # excepted if otherwise specified ...
        else error("in @usercut($(join(args,','))), invalid keyword argument: $(kwarg)")
        end
    end

    #cbdata = esc(cbdata)
    if VERSION < v"0.5.0-dev+3231"
        x = comparison_to_call(x)
    end
    if isexpr(x, :call) && length(x.args) == 3 # simple comparison
        lhs = :($(x.args[2]) - $(x.args[3])) # move everything to the lhs
        newaff, parsecode = parseExprToplevel(lhs, :aff)
        sense, vectorized = _canonicalize_sense(x.args[1])
        vectorized && error("Cannot add vectorized constraint in cut callback")
        quote
            aff = zero(AffExpr)
            $parsecode
            constr = constructconstraint!($newaff, $(quot(sense)))
            addusercut($cbdata, constr, localcut=$localcut_val)
        end
    else
        error("Syntax error in @usercut, expected one-sided comparison.")
    end
end

function addusercut(cbdata::MathProgBase.MathProgCallbackData, constr::LinearConstraint; localcut::Bool=false)
    if length(constr.terms.vars) == 0
        MathProgBase.cbaddcut!(cbdata, Cint[], Float64[], sensemap[sense(constr)], rhs(constr))
        return
    end
    assert_isfinite(constr.terms)
    m::Model = constr.terms.vars[1].m
    indices, coeffs = merge_duplicates(Cint, constr.terms, m.indexedVector, m)
    if localcut
        MathProgBase.cbaddcutlocal!(cbdata, indices, coeffs, sensemap[sense(constr)], rhs(constr))
    else
        MathProgBase.cbaddcut!(cbdata, indices, coeffs, sensemap[sense(constr)], rhs(constr))
    end
end

## User heuristic
export addsolution, setsolutionvalue

addsolution(cbdata::MathProgBase.MathProgCallbackData) = MathProgBase.cbaddsolution!(cbdata)
function setsolutionvalue(cbdata::MathProgBase.MathProgCallbackData, v::Variable, x)
    MathProgBase.cbsetsolutionvalue!(cbdata, convert(Cint, v.col), x)
end


text:
using IterTools
using SparseArrays
using StatsBase
using LinearAlgebra

include("visualize.jl")

mutable struct PredatorPreyHexWorld
    G
    m::Int64
    n::Int64
    alpha::Float64
    num_vertices::Int64
    L::Int64

    prey_penalty::Float64
    actions
    S
    A

    p
    r
    R
    T
    U
    discount_factor

    """Create a new hex world
    """
    function PredatorPreyHexWorld(num_predators::Int64, num_prey::Int64, alpha::Float64; prey_penalty::Float64=1.0)
        default_graph =
            [Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [2, 3, 4, 5, 6, 7]),
            Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [9, 10, 3, 1, 7, 8]),
            Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [10, 11, 12, 4, 1, 2]),
            Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [3, 12, 13, 14, 5, 1]),
            Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [1, 4, 14, 15, 16, 6]),
            Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [7, 1, 5, 16, 17, 18]),
            Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> [8, 2, 1, 6, 18, 19]),
            Dict(["NE", "NW", "W", "SW", "SE"] .=> [9, 2, 7, 19, 20]),
            Dict(["NW", "W", "SW"] .=> [10, 2, 8]),
            Dict(["NW", "W", "SW", "SE"] .=> [11, 3, 2, 9]),
            Dict(["W", "SW", "SE"] .=> [12, 3, 10]),
            Dict(["E", "W", "SW", "SE"] .=> [11, 13, 4, 3]),
            Dict(["E", "SW", "SE"] .=> [12, 14, 4]),
            Dict(["E", "NE", "SW", "SE"] .=> [4, 13, 15, 5]),
            Dict(["E", "NE", "SE"] .=> [5, 14, 16]),
            Dict(["E", "NE", "NW", "SE"] .=> [6, 5, 15, 17]),
            Dict(["E", "NE", "NW"] .=> [18, 8, 16]),
            Dict(["E", "NE", "NW", "W"] .=> [19, 7, 6, 17]),
            Dict(["E", "NE", "NW", "W"] .=> [20, 8, 7, 18]),
            Dict(["NW", "W"] .=> [8, 19])]

        return new(
            default_graph,
            num_predators,
            num_prey,
            alpha,
            length(default_graph),
            num_predators + num_prey,
            prey_penalty,
            ["E", "NE", "NW", "W", "SW", "SE"],
            vec(collect(product([1:length(default_graph) for _ in 1:(num_predators + num_prey)]...))),
            vec(collect(product([["E", "NE", "NW", "W", "SW", "SE"] for _ in 1:(num_predators + num_prey)]...))),
            nothing,
            nothing,
            nothing,
            nothing,
            nothing,
            0.9
        )
    end
end

function to_number(action)
    d = Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> 1:6)
    m = Array{Int64}(undef, length(action))
    for (i, direction) in enumerate(action)
        m[i] = d[direction]
    end
    return m
end

function get_index(a)
    s = 0
    l = 0
    if isa(a, String)
        return Dict(["E", "NE", "NW", "W", "SW", "SE"] .=> 1:6)[a]
    elseif isa(a[1], String)
        a = to_number(a)
        l = 6
    else
        l = 20
    end
    for i = 1:length(a)
        s += (a[i] - 1) * l ^ (i - 1)
    end
    s += 1
    return s
end

"""Return available actions at loc
"""
function avail_actions(world::PredatorPreyHexWorld, loc)
    return keys(world.G[loc])
end

"""Return adjacent vertices of loc
"""
function adj(world::PredatorPreyHexWorld, loc)
    return values(world.G[loc])
end

function actions(world::PredatorPreyHexWorld, state)
    return vec(collect(product([avail_actions(world, state[i]) for i=1:length(state)]...)))
end

function result(world::PredatorPreyHexWorld, state, action)
    return [world.G[state[i]][a] for (i, a) in enumerate(action)]
end

""" Return a vector of predators currently in `location`
"""
function predators_at(world::PredatorPreyHexWorld, state, location)
    r = Vector{Int64}()
    for i=1:world.m
        if state[i] == location && !(i in r)
            push!(r, i)
        end
    end
    return r
end

""" Return a vector of prey currently in `location`
"""
function prey_at(world::PredatorPreyHexWorld, state, location)
    r = Vector{Int64}()
    for i=(world.m + 1):world.L
        if state[i] == location && !(i in r)
            push!(r, i)
        end
    end
    return r
end

"""Reward for `i` if it executes `action` under `state`
"""
function reward(world::PredatorPreyHexWorld, i, state, action)
    if !(action in actions(world, state))
        return 0
    end

    nstate = result(world, state, action)
    loc = nstate[i]
    predators = predators_at(world, nstate, loc)
    preys = prey_at(world, nstate, loc)

    # i is predator
    if i <= world.m
        return -1 + length(preys) / length(predators)
    # i is prey
    else
        return length(predators) == 0 ? 0 : -world.prey_penalty
    end
end

""" Obtain a random location in the hex world where no agent is currently
occupying
"""
function free_location(world::PredatorPreyHexWorld, state)
    occupied = sort(unique(state[1:world.L]))
    list = []
    for i=1:world.num_vertices
        if !(i in occupied)
            push!(list, i)
        end
    end
    return list[rand(1:length(list))]
end

"""Return the probability distribution over the actions
"""
function predator_policy(world::PredatorPreyHexWorld, i, state)
    available_actions = avail_actions(world, state[i])
    v = world.G[state[i]]
    k = 0
    n = 0

    d = Dict(world.actions .=> [0.0 for _=1:length(world.actions)])

    for (a, ns) in v
        count = length(prey_at(world, state, ns))
        d[a] = count
        k += count
        n += count
        if count == 0
            n += 1
        end
    end

    for a in world.actions
        if !(a in available_actions)
            continue
        end

        # No prey around
        if k == 0
            # Each available cell gets equal probability
            d[a] = 1 / n
        else
            if d[a] != 0
                d[a] = d[a] * world.alpha / k
            else
                d[a] = (1 - world.alpha) / (n - k)
            end
        end
    end
    return d
end

function prey_policy(world::PredatorPreyHexWorld, i, state)
    d = Dict(world.actions .=> [0.0 for _=1:length(world.actions)])
    available_actions = avail_actions(world, state[i])
    v = world.G[state[i]]

    k = 0
    n = 0
    for (a, ns) in v
        count = length(predators_at(world, state, ns))
        d[a] +=  count
        if count == 0
            k += 1
        else
            n += 1 / count
        end
    end

    for a in world.actions
        if !(a in available_actions)
            continue
        end

        # All cells are safe
        if n == 0
            d[a] = 1 / length(available_actions)
        else
            if d[a] == 0
                d[a] = world.alpha / k
            else
                d[a] = (1 - world.alpha) / (d[a] * n)
            end
        end
    end
    return d
end

"""Probability distribution of the actions that agent `i` will choose at
`state`

Return a `Dict{String, Float64}` `d` where:
* `d["E"]`: the probability that agent `i` will select action `E` at this state
* `d["NW"]`: the probability that agent `i` will select action `NW` at this
state
* and so on
"""
function common_policy(world::PredatorPreyHexWorld, i, state)
    if i <= world.m
        return predator_policy(world, i, state)
    else
        return prey_policy(world, i, state)
    end
end

function is_applicable(world::PredatorPreyHexWorld, next_state, state, action)
    for i=1:length(state)
        e = avail_actions(world, state[i])
        v = adj(world, state[i])

        if !(next_state[i] in v) || !(e in action)
            return false
        end
    end
    return true
end

function compute_p(world::PredatorPreyHexWorld)
    m = Array{Float64}(undef, length(world.A), length(world.S), length(world.actions), world.L)
    p(a, s, act, i) =
        act in a ? prod([(j == i ? 1 : common_policy(world, j, s)[a[j]]) for j=1:world.L]) : 0
    for i=1:world.L
        for act in world.actions
            for a in world.A
                for s in world.S
                    m[get_index(a), get_index(s), get_index(act), i] = p(a, s, act, i)
                end
            end
        end
    end
    return m
end

function compute_r(world::PredatorPreyHexWorld)
    m = Array{Float64}(undef, length(world.A), length(world.S), world.L)
    for i in 1:world.L
        m[:, :, i] = [reward(world, i, s, a) for a in world.A, s in world.S]
    end
    return m
end

function compute_R(world::PredatorPreyHexWorld)
    r = world.r
    p = world.p

    res = Array{Float64}(undef, length(world.S), length(world.actions), world.L)
    for i=1:world.L
        for act in world.actions
            for s in world.S
                sum = 0
                for a in world.A
                    sum += p[get_index(a), get_index(s), get_index(act), i] * r[get_index(a), get_index(s), i]
                end
                res[get_index(s), get_index(act), i] = sum
            end
        end
    end
    return res
end

function compute_T(world::PredatorPreyHexWorld)
    p = world.p
    t = Array{Float64}(undef, length(world.S), length(world.S), length(world.actions), world.L)
    for i=1:world.L
        for act in world.actions
            for ns in world.S
                for s in world.S
                    sum = 0
                    guarantee = [a for a in actions(world, s) if a[i] == act]
                    for a in guarantee
                        sum += p[get_index(a), get_index(s), get_index(act), i]
                    end
                    t[get_index(ns), get_index(s), get_index(act), i] = sum
                end
            end
        end
    end
    return t
end

function compute_U(world::PredatorPreyHexWorld)
    U = Array{Float64}(undef, length(world.S), length(world.actions), world.L)

    for i=1:world.L
        for act in world.actions
            U[:, get_index(act), i] = (I - world.discount_factor * world.T[:, :, get_index(act), i]) \ world.R[:, get_index(act), i]
        end
    end
    return U
end

function compute!(world::PredatorPreyHexWorld)
    world.p = compute_p(world)
    println("Done p, $(size(world.p))")
    world.r = compute_r(world)
    println("Done r, $(size(world.r))")
    world.R = compute_R(world)
    println("Done R, $(size(world.R))")
    world.T = compute_T(world)
    println("Done T, $(size(world.T))")
    world.U = compute_U(world)
    println("Done U, $(size(world.U))")
end

function best_response(world::PredatorPreyHexWorld, state, i)
    max = -Inf
    select = []
    for a in avail_actions(world, state[i])
        u =  world.U[get_index(state), get_index(a), i]
        if u > max
            max = u
            select = [a]
        elseif u == max
            push!(select, a)
        end
    end
    return sample(select)
end

function next_action(world::PredatorPreyHexWorld, state)
    return [best_response(world, state, i) for i in 1:world.L]
end

function apply!(world::PredatorPreyHexWorld, state, action)
    new_state = result(world, state, action)

    prey_locations = state[(world.m + 1):world.L]
    for loc in prey_locations
        prey = prey_at(world, new_state, loc)

        if length(prey) != 0
            for p in prey
                new_state[p + world.m] = free_location(world, new_state)
            end
        end
    end

    return new_state
end

function run(world::PredatorPreyHexWorld, init_state, folder, times::Int64)
    compute!(world)
    state = init_state
    for i=1:times
        visualize_state("assets/$folder/step$(lpad(i, 3, "0")).png", state, world.m, world.n, i)
        action = next_action(world, state)
        next_state = apply!(world, state, action)
        state = next_state
    end
end

function test()
    num_tests = 2
    times = 30
    folders = ["test1", "test2"]
    alpha = [0.75, 1.0]
    prey_penalty = 100.0

    for i=1:num_tests
        println("Running test $i...")
        world = PredatorPreyHexWorld(1, 1, alpha[i]; prey_penalty=prey_penalty)
        init_state = sample(1:world.num_vertices, world.L, replace=false)
        run(world, init_state, folders[i], times)
    end
end

test()


text:
using DEBmicroTrait
using CSV, DataFrames

df_metabolites  = CSV.read("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/exudation_properties.csv", DataFrame, missingstring="N/A")

N_metabolites      = 83
dGcox              = zeros(N_metabolites)
dGcat              = zeros(N_metabolites)
dGAn               = zeros(N_metabolites)
λ_base             = zeros(N_metabolites)
N_C                = zeros(N_metabolites)
eta                = zeros(N_metabolites)
chemFormBiom       = [1, 1.64, 0.2, 0.38, 0, 0, 0]

for i in 1:N_metabolites
    elementstring = convert(String,df_metabolites.Formula[i])
    N_C[i] = DEBmicroTrait.extract_composition(elementstring)[1]
    out = DEBmicroTrait.get_lambda(elementstring, chemFormBiom)
    dGcox[i] = out[2][3]
    dGcat[i] = out[2][5]
    dGAn[i]  = out[2][8]
    λ_base[i]     = out[1][1]
    eta      = @. dGAn/(λ_base*dGcat)
end

df_out = DataFrame()
df_out.monomer = df_metabolites.Name
df_out.lambda = λ_base
df_out.eta = eta
df_out.ontology = df_metabolites.Ontology

CSV.write("/Users/glmarschmann/.julia/dev/DEBmicroTrait/files/exudates_eta.csv", df_out)

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


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

module CoreDocs

import ..esc, ..push!, ..getindex, ..unsafe_load, ..Csize_t, ..@nospecialize

function doc!(source::LineNumberNode, mod::Module, str, ex)
    @nospecialize str ex
    push!(DOCS, Core.svec(mod, ex, str, source.file, source.line))
    nothing
end
const DOCS = Array{Core.SimpleVector,1}()

isexpr(@nospecialize(x), h::Symbol) = isa(x, Expr) && x.head === h

lazy_iterpolate(s::AbstractString) = Expr(:call, Core.svec, s)
lazy_iterpolate(@nospecialize x) = isexpr(x, :string) ? Expr(:call, Core.svec, x.args...) : x

function docm(source::LineNumberNode, mod::Module, str, x)
    out = Expr(:call, doc!, QuoteNode(source), mod, lazy_iterpolate(str), QuoteNode(x))
    if isexpr(x, :module)
        out = Expr(:toplevel, out, x)
    elseif isexpr(x, :call)
    else
        out = Expr(:block, x, out)
    end
    return esc(out)
end
docm(source::LineNumberNode, mod::Module, x) =
    isexpr(x, :->) ? docm(source, mod, x.args[1], x.args[2].args[2]) : error("invalid '@doc'.")

end

max_stars_repo_name: helgee/julia

max_issues_repo_name: helgee/julia

max_forks_repo_name: helgee/julia


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

JLLWrappers.@generate_main_file_header("GAP_pkg_json")
JLLWrappers.@generate_main_file("GAP_pkg_json", UUID("b348cf08-4407-5d37-b307-1564396f8b07"))
end  # module GAP_pkg_json_jll


text:
using Documenter, LuxurySparse

# TODO: use Literate to process examples
# using Literate
# preprocess tutorial scripts

# make documents
makedocs(
    modules = [LuxurySparse],
    clean = false,
    format = :html,
    sitename = "LuxurySparse.jl",
    linkcheck = !("skiplinks" in ARGS),
    analytics = "UA-89508993-1",
    pages = [
        "Home" => "index.md",
        "Tutorial" => "tutorial.md",
        "Manual" => "luxurysparse.md",
    ],
    html_prettyurls = !("local" in ARGS),
    html_canonical = "https://quantumbfs.github.io/LuxurySparse.jl/latest/",
)

deploydocs(
    repo = "github.com/QuantumBFS/LuxurySparse.jl.git",
    target = "build",
    julia = "1.0",
    deps = nothing,
    make = nothing,
)


text:
# ### Categorical Histogram Plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

# https://jp.mathworks.com/matlabcentral/fileexchange/63452-matlab-plot-gallery-categorical-histogram-plot?focused=7852078&tab=example

using StatsBase, Random

Random.seed!(2018)
s = ["Poor", "Fair", "Good", "Excellent"]
data = sample(s, Weights([1,1,2,2]), 100)
datamap = countmap(data)
bar((x -> datamap[x]).(s), xticks=(1:4, s))