text:
# Format forced into left/right children, and additional markov probability matrix functionality added.
mutable struct PhyloNode
    label::String
    left::Nullable{PhyloNode}
    right::Nullable{PhyloNode}
    length::Real
    prob::Array{Float64,2}
end

const EMPTY_NODE = PhyloNode("", Nullable{PhyloNode}(), Nullable{PhyloNode}(), 0.0, Array{Float64}(0,0))

struct PhyloTree
   root::PhyloNode
   order::Vector{String}
   index::Dict{String,Int}
end

const EMPTY_TREE = PhyloTree(EMPTY_NODE, String[], Dict{String,Int}())

### 'parsenewick' function based on: https://gist.github.com/porterjamesj/7672080
### License is CC-BY-NC
function parsenewick(newick::String)
    newick = rstrip(newick,';')
    newick = replace(newick,":","+")
    newick_expr = parse(newick)
    root  = parsenewick(newick_expr)
    order = collect(root)
    index = index_dict(order)
    PhyloTree(root, order, index)
end

parsenewick(newick::Symbol) = PhyloNode(string(newick), Nullable{PhyloNode}(), Nullable{PhyloNode}(), -1, Array{Float64}(0,0))

function parsenewick(newick::Expr)
    if newick.head == :tuple
        children = [parsenewick(child) for child in newick.args]
        left  = Nullable(children[1])
        right = Nullable(children[2])
        name = ""
        length = -1
    elseif newick.head == :call
        if newick.args[1] == :+
            # + indicates length
            length = newick.args[3]
            if typeof(newick.args[2]) == Expr
                if newick.args[2].head == :tuple
                    children = [parsenewick(child) for child in newick.args[2].args]
                    left  = Nullable(children[1])
                    right = Nullable(children[2])
                    name = ""
                elseif newick.args[2].head == :call && newick.args[2].args[1] == :*
                    # * indicates naming
                    name = string(newick.args[2].args[3])
                    children = [parsenewick(child) for child in newick.args[2].args[2].args]
                    left  = Nullable(children[1])
                    right = Nullable(children[2])
                end
            elseif typeof(newick.args[2]) == Symbol || typeof(newick.args[2]) == Int
                # tip node
                name = string(newick.args[2])
                left  = Nullable{PhyloNode}()
                right = Nullable{PhyloNode}()
            end
        elseif newick.args[1] == :*
            # bare * indicates a node with name but no length
            name = string(newick.args[3])
            children = [parsenewick(child) for child in newick.args[2].args]
            left  = Nullable(children[1])
            right = Nullable(children[2])
            length = -1
        end
    end
    PhyloNode(name,left,right,length, Array{Float64}(0,0))
end

# Tree collection
function Base.collect( root::PhyloNode )
   res = Vector{String}()
   function dfs_names!( vec::Vector{String}, node::PhyloNode )
      if !isnull(node.left) && !isnull(node.right)
         dfs_names!( vec, node.left.value )
         dfs_names!( vec, node.right.value )
      else
         push!( vec, node.label )
      end
      vec
   end
   dfs_names!( res, root )
end

function index_dict( labels::Vector{String} )
   ret = Dict{String,Int}()
   for i in 1:length(labels)
      ret[labels[i]] = i
   end
   ret
end

# Matrix setting and branch length modification functions
# Set P = exp(Q*b) for each node
set_prob_mat!(tree::PhyloTree, Q::Array{Float64,2}) = set_prob_mat!(tree.root, Q)

function set_prob_mat!(node::PhyloNode, Q::Array{Float64,2})
   if !isnull(node.left) && !isnull(node.right) # internal node
      left  = set_prob_mat!(node.left.value,  Q)
      right = set_prob_mat!(node.right.value, Q)
   end
   if node.length > 0
      node.prob = expm(Q*node.length)
   end
   node
end

extend_branches!(tree::PhyloTree, multiplier::Float64) = extend_branches!(tree.root, multiplier)

function extend_branches!(node::PhyloNode, multiplier::Float64)
   if !isnull(node.left) && !isnull(node.right) # internal node
      left  = extend_branches!(node.left.value, multiplier)
      right = extend_branches!(node.right.value, multiplier)
   end
   if node.length > 0
      node.length *= multiplier
   end
   node
end

# My implementation of Felsenstein's pruning algorithm
# to calculate the L(D|T) = L(k) = Sigma_x(pi(x)*L(x))
# from the tree topology, branch lengths, instantaneous matrix (Q),
# and the data at site i
function likelihood( tree::PhyloTree, single::Vector{DNA};
                     background=GTR_SINGLE_PI, gapdenom::Float64=1.0 )
   function _likelihood( node::PhyloNode, single::Vector{DNA}, branchlength::Bool )
      if !isnull(node.left) && !isnull(node.right)
          if branchlength
             const left_res  = node.prob * _likelihood( node.left.value,  single, true )
             const right_res = node.prob * _likelihood( node.right.value, single, true )
          else
             const left_res  = _likelihood( node.left.value,  single, true )
             const right_res = _likelihood( node.right.value, single, true )
          end
          return left_res .* right_res
      else
         const ind = tree.index[node.label]
         const v = isgap( single[ind] ) ? (return ones(4) / gapdenom) : zeros(4)
         v[ trailing_zeros(reinterpret(UInt8, single[ind])) + 1 ] = 1.0
         return node.prob * v
      end
   end
   sum(_likelihood( tree.root, single, false ) .* background)
end

function likelihood( tree::PhyloTree, first::Vector{DNA}, last::Vector{DNA};
                     background=GTR_PAIRED_PI )
   index( l::DNA, r::DNA ) = index( reinterpret(UInt8, l), reinterpret(UInt8, r) )
   function index( l::UInt8, r::UInt8 )
      const lidx = trailing_zeros(l) << 2
      const ridx = trailing_zeros(r)
      (lidx | ridx) + 1
   end
   function _likelihood( node::PhyloNode, first::Vector{DNA}, last::Vector{DNA}, branchlength::Bool )
      if !isnull(node.left) && !isnull(node.right)
          if branchlength
             const left_res  = node.prob * _likelihood( node.left.value,  first, last, true )
             const right_res = node.prob * _likelihood( node.right.value, first, last, true )
          else
             const left_res  = _likelihood( node.left.value,  first, last, true )
             const right_res = _likelihood( node.right.value, first, last, true )
          end
          return left_res .* right_res
      else
         const ind = tree.index[node.label]
         const v = isgap(first[ind]) || isgap(last[ind])  ? (return ones(16)) : zeros(16)
         v[ index(first[ind], last[ind]) ] = 1.0
         return node.prob * v
      end
   end
   sum(_likelihood( tree.root, first, last, false ) .* background)
end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    LearnedModel(m, θ)

An object that stores a learning `model`
along with its learned parameters `θ`.
"""
struct LearnedModel
  model
  θ
end

"""
    learn(task, sdata, model)

Learn the `task` with `sdata` using a learning `model`.
"""
function learn(task::AbstractLearningTask, sdata, model)
  if issupervised(task)
    X = view(sdata, collect(features(task)))
    y = sdata[label(task)]
    θ, _, __ = MI.fit(model, 0, X, y)
  else
    X = view(sdata, collect(features(task)))
    θ, _, __ = MI.fit(model, 0, X)
  end

  LearnedModel(model, θ)
end

"""
    perform(task, sdata, lmodel)

Perform the `task` with `sdata` using a *learned* `lmodel`.
"""
function perform(task::AbstractLearningTask, sdata, lmodel)
  # unpack model and learned parameters
  model, θ = lmodel.model, lmodel.θ

  # apply model to the data
  X = view(sdata, collect(features(task)))
  ŷ = MI.predict(model, θ, X)

  # post-process result
  var = outputvars(task)[1]
  val = if issupervised(task)
    isprobabilistic(model) ? mode.(ŷ) : ŷ
  else
    ŷ
  end

  georef((; var=>val), domain(sdata))
end


text:
@base type Options
    N::Int #Number of years to calculate (from 2010 onwards)
    tstep::Int #Years per Period
    α::Float64 #Elasticity of marginal utility of consumption
    ρ::Float64 #Initial rate of social time preference per year ρ
    γₑ::Float64 #Capital elasticity in production function (Capital share)
    pop₀::Int #Initial world population (millions)
    popadj::Float64 #Growth rate to calibrate to 2050 pop projection
    popasym::Int #Asymptotic population (millions)
    δk::Float64 #Depreciation rate on capital (per year)
    q₀::Float64 #Initial world gross output (trill 2005 USD)
    k₀::Float64 #Initial capital value (trill 2005 USD)
    a₀::Float64 #Initial level of total factor productivity
    ga₀::Float64 #Initial growth rate for TFP per 5 years
    δₐ::Float64 #Decline rate of TFP per 5 years
    gσ₁::Float64 #Initial growth of sigma (continuous per year)
    δσ::Float64 #Decline rate of decarbonization per period
    eland₀::Float64 #Carbon emissions from land 2010 (GtCO2 per year)
    deland::Float64 #Decline rate of land emissions (per period)
    mat₀::Float64 #Initial Concentration in atmosphere 2010 (GtC)
    mu₀::Float64 #Initial Concentration in upper strata 2010 (GtC)
    ml₀::Float64 #Initial Concentration in lower strata 2010 (GtC)
    mateq::Float64 #Equilibrium concentration atmosphere  (GtC)
    mueq::Float64 #Equilibrium concentration in upper strata (GtC)
    mleq::Float64 #Equilibrium concentration in lower strata (GtC)
    ϕ₁₂::Float64 #Carbon cycle transition matrix coefficient
    ϕ₂₃::Float64 #Carbon cycle transition matrix coefficient
    t2xco2::Float64 #Equilibrium temp impact (oC per doubling CO2)
    fₑₓ0::Float64 #2010 forcings of non-CO2 GHG (Wm-2)
    fₑₓ1::Float64 #2100 forcings of non-CO2 GHG (Wm-2)
    tocean₀::Float64 #Initial lower stratum temp change (C from 1900)
    tatm₀::Float64 #Initial atmospheric temp change (C from 1900)
    ξ₁::Float64 #Climate equation coefficient for upper level
    ξ₃::Float64 #Transfer coefficient upper to lower stratum
    ξ₄::Float64 #Transfer coefficient for lower level
    η::Float64 #Forcings of equilibrium CO2 doubling (Wm-2)
    ψ₁::Float64 #Damage intercept
    ψ₂::Float64 #Damage quadratic term
    ψ₃::Float64 #Damage exponent
    θ₂::Float64 #Exponent of control cost function
    pback::Float64 #Cost of backstop 2005$ per tCO2 2010
    gback::Float64 #Initial cost decline backstop cost per period
    limμ::Float64 #Upper limit on control rate after 2150
    fosslim::Float64 #Maximum cumulative extraction fossil fuels (GtC)
    scale1::Float64 #Multiplicative scaling coefficient
    scale2::Float64 #Additive scaling coefficient
end

@base type Parameters
    ϕ₁₁::Float64 # Carbon cycle transition matrix coefficient
    ϕ₂₁::Float64 # Carbon cycle transition matrix coefficient
    ϕ₂₂::Float64 # Carbon cycle transition matrix coefficient
    ϕ₃₂::Float64 # Carbon cycle transition matrix coefficient
    ϕ₃₃::Float64 # Carbon cycle transition matrix coefficient
    σ₀::Float64 # Carbon intensity 2010 (kgCO2 per output 2005 USD 2010)
    λ::Float64 # Climate model parameter
    gₐ::Array{Float64,1} # Growth rate of productivity from 0 to N
    Etree::Array{Float64,1} # Emissions from deforestation
    L::Array{Float64,1} # Level of population and labor
    A::Array{Float64,1} # Level of total factor productivity
    gσ::Array{Float64,1} # Change in sigma (cumulative improvement of energy efficiency)
    σ::Array{Float64,1} # CO2-equivalent-emissions output ratio
    θ₁::Array{Float64,1} # Adjusted cost for backstop
    fₑₓ::Array{Float64,1} # Exogenous forcing for other greenhouse gases
end

@base type Variables
    μ::Array{JuMP.Variable,1} # Emission control rate GHGs
    FORC::Array{JuMP.Variable,1} # Increase in radiative forcing (watts per m2 from 1900)
    Tₐₜ::Array{JuMP.Variable,1} # Increase temperature of atmosphere (degrees C from 1900)
    Tₗₒ::Array{JuMP.Variable,1} # Increase temperatureof lower oceans (degrees C from 1900)
    Mₐₜ::Array{JuMP.Variable,1} # Carbon concentration increase in atmosphere (GtC from 1750)
    Mᵤₚ::Array{JuMP.Variable,1} # Carbon concentration increase in shallow oceans (GtC from 1750)
    Mₗₒ::Array{JuMP.Variable,1} # Carbon concentration increase in lower oceans (GtC from 1750)
    E::Array{JuMP.Variable,1} # Total CO2 emissions (GtCO2 per year)
    C::Array{JuMP.Variable,1} # Consumption (trillions 2005 US dollars per year)
    K::Array{JuMP.Variable,1} # Capital stock (trillions 2005 US dollars)
    CPC::Array{JuMP.Variable,1} # Per capita consumption (thousands 2005 USD per year)
    I::Array{JuMP.Variable,1} # Investment (trillions 2005 USD per year)
    S::Array{JuMP.Variable,1} # Gross savings rate as fraction of gross world product
    RI::Array{JuMP.Variable,1} # Real interest rate (per annum)
    Y::Array{JuMP.Variable,1} # Gross world product net of abatement and damages (trillions 2005 USD per year)
    YGROSS::Array{JuMP.Variable,1} # Gross world product GROSS of abatement and damages (trillions 2005 USD per year)
    YNET::Array{JuMP.Variable,1} # Output net of damages equation (trillions 2005 USD per year)
    DAMAGES::Array{JuMP.Variable,1} # Damages (trillions 2005 USD per year)
    MCABATE::Array{JuMP.Variable,1} # Marginal cost of abatement (2005$ per ton CO2)
    CCA::Array{JuMP.Variable,1} # Cumulative industrial carbon emissions (GTC)
    PERIODU::Array{JuMP.Variable,1} # One period utility function
    UTILITY::JuMP.Variable # Welfare function
end

#TODO: Consider extending this to all equations rather than the ones useful for results
@base type Equations
    eeq::Array{JuMP.ConstraintRef,1} # Emissions Equation
end

#NOTE: Monetary units herein are based on differing years depinding on version.
# For example: v2013R uses 2005 as its base year and v2016R uses 2010.
# Rule of thumb for these values is that the starting value `Results.years[1]`
# is the baseline.
@base type Results
    years::Array{Int64,1} #Simulation units to true year values
    Mₐₜ::Array{Float64,1} # Carbon concentration increase in atmosphere (GtC from 1750)
    Mₐₜppm::Array{Float64,1} # Carbon concentration increase in atmosphere (ppm from 1750)
    Mᵤₚ::Array{Float64,1} # Carbon concentration increase in shallow oceans (GtC from 1750)
    Mₗₒ::Array{Float64,1} # Carbon concentration increase in lower oceans (GtC from 1750)
    CCA::Array{Float64,1} # Cumulative industrial carbon emissions (GTC)
    CCAratio::Array{Float64,1} # Cumulative industrial emissions ratio to maximum possible
    Tₐₜ::Array{Float64,1} # Increase temperature of atmosphere (degrees C from 1900)
    FORC::Array{Float64,1} # Increase in radiative forcing (watts per m2 from 1900)
    Tₗₒ::Array{Float64,1} # Increase temperatureof lower oceans (degrees C from 1900)
    YGROSS::Array{Float64,1} # Gross world product GROSS of abatement and damages (trillions USD per year)
    DAMAGES::Array{Float64,1} # Damages (trillions USD per year)
    YNET::Array{Float64,1} # Output net of damages equation (trillions USD per year)
    Y::Array{Float64,1} # Gross world product net of abatement and damages (trillions USD per year)
    E::Array{Float64,1} # Total CO2 emissions (GtCO2 per year)
    I::Array{Float64,1} # Investment (trillions USD per year)
    K::Array{Float64,1} # Capital stock (trillions US dollars)
    MPK::Array{Float64,1} #First Period gross MPK
    C::Array{Float64,1} # Consumption (trillions US dollars per year)
    CPC::Array{Float64,1} # Per capita consumption (thousands USD per year)
    PERIODU::Array{Float64,1} # One period utility function
    UTILITY::Float64 # Welfare function
    S::Array{Float64,1} # Gross savings rate as fraction of gross world product
    μ::Array{Float64,1} # Emission control rate GHGs
    RI::Array{Float64,1} # Real interest rate (per annum)
    scc::Array{Float64,1} # Social Cost of Carbon
end


text:
module TestIsLess
using Test
using CategoricalArrays

pool = CategoricalPool([1, 2, 3])

v1 = CategoricalValue(1, pool)
v2 = CategoricalValue(2, pool)
v3 = CategoricalValue(3, pool)

@testset "values from unordered CategoricalPool" begin
    @test_throws ArgumentError v1 < v1
    @test_throws ArgumentError v1 < v2
    @test_throws ArgumentError v1 < v3
    @test_throws ArgumentError v2 < v1
    @test_throws ArgumentError v2 < v2
    @test_throws ArgumentError v2 < v3
    @test_throws ArgumentError v3 < v1
    @test_throws ArgumentError v3 < v2
    @test_throws ArgumentError v3 < v3

    @test_throws ArgumentError v1 <= v1
    @test_throws ArgumentError v1 <= v2
    @test_throws ArgumentError v1 <= v3
    @test_throws ArgumentError v2 <= v1
    @test_throws ArgumentError v2 <= v2
    @test_throws ArgumentError v2 <= v3
    @test_throws ArgumentError v3 <= v1
    @test_throws ArgumentError v3 <= v2
    @test_throws ArgumentError v3 <= v3

    @test_throws ArgumentError v1 > v1
    @test_throws ArgumentError v1 > v2
    @test_throws ArgumentError v1 > v3
    @test_throws ArgumentError v2 > v1
    @test_throws ArgumentError v2 > v2
    @test_throws ArgumentError v2 > v3
    @test_throws ArgumentError v3 > v1
    @test_throws ArgumentError v3 > v2
    @test_throws ArgumentError v3 > v3

    @test_throws ArgumentError v1 >= v1
    @test_throws ArgumentError v1 >= v2
    @test_throws ArgumentError v1 >= v3
    @test_throws ArgumentError v2 >= v1
    @test_throws ArgumentError v2 >= v2
    @test_throws ArgumentError v2 >= v3
    @test_throws ArgumentError v3 >= v1
    @test_throws ArgumentError v3 >= v2
    @test_throws ArgumentError v3 >= v3

    @test isless(v1, v1) === false
    @test isless(v1, v2) === true
    @test isless(v1, v3) === true
    @test isless(v2, v1) === false
    @test isless(v2, v2) === false
    @test isless(v2, v3) === true
    @test isless(v3, v1) === false
    @test isless(v3, v2) === false
    @test isless(v3, v3) === false

    @testset "comparison with values of different types" begin
        @test isless(v1, 1) === false
        @test isless(v1, 2) === true
        @test_throws KeyError isless(v1, 10)
        @test_throws KeyError isless(v1, "a")
        @test isless(1, v1) === false
        @test_throws KeyError isless("a", v1)
        @test_throws ArgumentError v1 < 1
        @test_throws ArgumentError v1 < 10
        @test_throws ArgumentError v1 < 2
        @test_throws ArgumentError v1 < "a"
        @test_throws ArgumentError v1 <= 1
        @test_throws ArgumentError v1 <= 2
        @test_throws ArgumentError v1 <= "a"
        @test_throws ArgumentError v1 > 1
        @test_throws ArgumentError v1 > 2
        @test_throws ArgumentError v1 > "a"
        @test_throws ArgumentError v1 >= 1
        @test_throws ArgumentError v1 >= 2
        @test_throws ArgumentError v1 >= "a"
    end

    @testset "comparison with missing" begin
        @test isless(v1, missing)
        @test !isless(missing, v1)
        @test ismissing(v1 < missing)
        @test ismissing(v1 <= missing)
        @test ismissing(v1 > missing)
        @test ismissing(v1 >= missing)
        @test ismissing(missing < v1)
        @test ismissing(missing <= v1)
        @test ismissing(missing > v1)
        @test ismissing(missing >= v1)
    end
end

@testset "values from ordered CategoricalPool" begin
    @test ordered!(pool, true) === pool
    @test isordered(pool) === true

    @test (v1 < v1) === false
    @test (v1 < v2) === true
    @test (v1 < v3) === true
    @test (v2 < v1) === false
    @test (v2 < v2) === false
    @test (v2 < v3) === true
    @test (v3 < v1) === false
    @test (v3 < v2) === false
    @test (v3 < v3) === false

    @test (v1 <= v1) === true
    @test (v1 <= v2) === true
    @test (v1 <= v3) === true
    @test (v2 <= v1) === false
    @test (v2 <= v2) === true
    @test (v2 <= v3) === true
    @test (v3 <= v1) === false
    @test (v3 <= v2) === false
    @test (v3 <= v3) === true

    @test (v1 > v1) === false
    @test (v1 > v2) === false
    @test (v1 > v3) === false
    @test (v2 > v1) === true
    @test (v2 > v2) === false
    @test (v2 > v3) === false
    @test (v3 > v1) === true
    @test (v3 > v2) === true
    @test (v3 > v3) === false

    @test (v1 >= v1) === true
    @test (v1 >= v2) === false
    @test (v1 >= v3) === false
    @test (v2 >= v1) === true
    @test (v2 >= v2) === true
    @test (v2 >= v3) === false
    @test (v3 >= v1) === true
    @test (v3 >= v2) === true
    @test (v3 >= v3) === true

    @test isless(v1, v1) === false
    @test isless(v1, v2) === true
    @test isless(v1, v3) === true
    @test isless(v2, v1) === false
    @test isless(v2, v2) === false
    @test isless(v2, v3) === true
    @test isless(v3, v1) === false
    @test isless(v3, v2) === false
    @test isless(v3, v3) === false

    @testset "comparison with values of different types" begin
        @test isless(v1, 1) === false
        @test isless(v1, 2) === true
        @test_throws KeyError isless(v1, 10)
        @test_throws KeyError isless(v1, "a")
        @test isless(1, v1) === false
        @test_throws KeyError isless("a", v1)
        @test (v1 < 1) === false
        @test (v1 < 2) === true
        @test_throws KeyError v1 < 10
        @test_throws KeyError v1 < "a"
        @test (v1 <= 1) === true
        @test (v1 <= 2) === true
        @test_throws KeyError v1 <= "a"
        @test (v1 > 1) === false
        @test (v1 > 2) === false
        @test_throws KeyError v1 > "a"
        @test (v1 >= 1) === true
        @test (v1 >= 2) === false
        @test_throws KeyError v1 >= "a"
    end

    @testset "comparison with missing" begin
        @test isless(v1, missing)
        @test !isless(missing, v1)
        @test ismissing(v1 < missing)
        @test ismissing(v1 <= missing)
        @test ismissing(v1 > missing)
        @test ismissing(v1 >= missing)
        @test ismissing(missing < v1)
        @test ismissing(missing <= v1)
        @test ismissing(missing > v1)
        @test ismissing(missing >= v1)
    end
end

@testset "comparisons with reordered levels" begin
    @test levels!(pool, [2, 3, 1]) === pool
    @test levels(pool) == [2, 3, 1]

    @test (v1 < v1) === false
    @test (v1 < v2) === false
    @test (v1 < v3) === false
    @test (v2 < v1) === true
    @test (v2 < v2) === false
    @test (v2 < v3) === true
    @test (v3 < v1) === true
    @test (v3 < v2) === false
    @test (v3 < v3) === false

    @test (v1 <= v1) === true
    @test (v1 <= v2) === false
    @test (v1 <= v3) === false
    @test (v2 <= v1) === true
    @test (v2 <= v2) === true
    @test (v2 <= v3) === true
    @test (v3 <= v1) === true
    @test (v3 <= v2) === false
    @test (v3 <= v3) === true

    @test (v1 > v1) === false
    @test (v1 > v2) === true
    @test (v1 > v3) === true
    @test (v2 > v1) === false
    @test (v2 > v2) === false
    @test (v2 > v3) === false
    @test (v3 > v1) === false
    @test (v3 > v2) === true
    @test (v3 > v3) === false

    @test (v1 >= v1) === true
    @test (v1 >= v2) === true
    @test (v1 >= v3) === true
    @test (v2 >= v1) === false
    @test (v2 >= v2) === true
    @test (v2 >= v3) === false
    @test (v3 >= v1) === false
    @test (v3 >= v2) === true
    @test (v3 >= v3) === true

    @test isless(v1, v1) === false
    @test isless(v1, v2) === false
    @test isless(v1, v3) === false
    @test isless(v2, v1) === true
    @test isless(v2, v2) === false
    @test isless(v2, v3) === true
    @test isless(v3, v1) === true
    @test isless(v3, v2) === false
    @test isless(v3, v3) === false

    @testset "comparison with values of different types" begin
        @test isless(v1, 1) === false
        @test isless(v1, 2) === false
        @test isless(v2, 1) === true
        @test_throws KeyError isless(v1, 10)
        @test_throws KeyError isless(v1, "a")
        @test isless(1, v1) === false
        @test isless(2, v1) === true
        @test_throws KeyError isless("a", v1)
        @test (v1 < 1) === false
        @test (v1 < 2) === false
        @test (v2 < 1) === true
        @test_throws KeyError v1 < 10
        @test_throws KeyError v1 < "a"
        @test (v1 <= 1) === true
        @test (v1 <= 2) === false
        @test (v2 <= 1) === true
        @test_throws KeyError v1 <= "a"
        @test (v1 > 1) === false
        @test (v1 > 2) === true
        @test (v2 > 1) === false
        @test_throws KeyError v1 > "a"
        @test (v1 >= 1) === true
        @test (v1 >= 2) === true
        @test (v2 >= 1) === false
        @test_throws KeyError v1 >= "a"
    end

    @test ordered!(pool, false) === pool
    @test isordered(pool) === false

    @test_throws ArgumentError v1 < v1
    @test_throws ArgumentError v1 < v2
    @test_throws ArgumentError v1 < v3
    @test_throws ArgumentError v2 < v1
    @test_throws ArgumentError v2 < v2
    @test_throws ArgumentError v2 < v3
    @test_throws ArgumentError v3 < v1
    @test_throws ArgumentError v3 < v2
    @test_throws ArgumentError v3 < v3

    @test_throws ArgumentError v1 <= v1
    @test_throws ArgumentError v1 <= v2
    @test_throws ArgumentError v1 <= v3
    @test_throws ArgumentError v2 <= v1
    @test_throws ArgumentError v2 <= v2
    @test_throws ArgumentError v2 <= v3
    @test_throws ArgumentError v3 <= v1
    @test_throws ArgumentError v3 <= v2
    @test_throws ArgumentError v3 <= v3

    @test_throws ArgumentError v1 > v1
    @test_throws ArgumentError v1 > v2
    @test_throws ArgumentError v1 > v3
    @test_throws ArgumentError v2 > v1
    @test_throws ArgumentError v2 > v2
    @test_throws ArgumentError v2 > v3
    @test_throws ArgumentError v3 > v1
    @test_throws ArgumentError v3 > v2
    @test_throws ArgumentError v3 > v3

    @test_throws ArgumentError v1 >= v1
    @test_throws ArgumentError v1 >= v2
    @test_throws ArgumentError v1 >= v3
    @test_throws ArgumentError v2 >= v1
    @test_throws ArgumentError v2 >= v2
    @test_throws ArgumentError v2 >= v3
    @test_throws ArgumentError v3 >= v1
    @test_throws ArgumentError v3 >= v2
    @test_throws ArgumentError v3 >= v3

    @test isless(v1, v1) === false
    @test isless(v1, v2) === false
    @test isless(v1, v3) === false
    @test isless(v2, v1) === true
    @test isless(v2, v2) === false
    @test isless(v2, v3) === true
    @test isless(v3, v1) === true
    @test isless(v3, v2) === false
    @test isless(v3, v3) === false
end

@testset "ordering comparisons between pools fail" begin
    pool2 = CategoricalPool([1, 2, 3])
    ordered!(pool2, true)

    v = CategoricalValue(1, pool2)

    @test_throws ArgumentError v < v1
    @test_throws ArgumentError v <= v1
    @test_throws ArgumentError v > v1
    @test_throws ArgumentError v >= v1
    @test_throws ArgumentError isless(v, v1)
end

end


text:
abstract type AbstractLoadForm end

abstract type AbstractControllableLoadForm <: AbstractLoadForm end

struct InterruptibleLoad <: AbstractControllableLoadForm end

include("electric_loads/load_variables.jl")
include("electric_loads/shedding_constraints.jl")
include("electric_loads/controlableload_cost.jl")


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

plot(x=[0.236033, 0.346517, 0.312707, 0.00790928, 0.488613, 0.210968, 0.951916, 0.999905, 0.251662, 0.986666],
     y=[0.555751, 0.437108, 0.424718, 0.773223, 0.28119, 0.209472, 0.251379, 0.0203749, 0.287702, 0.859512],
     Scale.x_continuous(format=:scientific),
     Scale.y_continuous(format=:plain))


text:
PJane = 25000
PJoe = 40000
P = PJane + PJoe
i = 0.1
m = 1
N = 6
iJane = 0.09

F = P * (1 + i / m)^(N * m)

iJoe = ((F - PJane * (1 + i)^N) / PJoe)^(1/N) - 1


text:
"""
Created on Tue May 19 2020

@author: Yoann Pradat

Functions used for assessing convergence of NMF algorithms.
"""

"""
    _update_metrics!(metrics, cur_iter, V, W, H)

Helper function for tracking cost evolution in NMF algorithms.
"""
function _update_metrics!(metrics::NMFMetrics, V::Matrix{T}, W::Matrix{T}, H::Matrix{T}, stopping_crit::Symbol) where T <: Real
    #### compute connectivity matrices only if required
    if stopping_crit == :conn
        prv_conn = metrics.cur_conn
        cur_conn  = _connectivity_matrix(H)
        push!(metrics.unchange_conn, prv_conn == cur_conn)
        metrics.cur_conn = cur_conn
    end

    #### compute indices and cost
    prv_cost = metrics.cur_cost

    #### compute new cost and new connectivity matrix
    cur_cost  = metrics.func_cost(V, W, H)
    cur_divg  = metrics.func_divg(V, W, H)
    rel_cost  = abs((cur_cost - prv_cost)/prv_cost)

    #### update metrics
    push!(metrics.absolute_cost, cur_cost)
    push!(metrics.absolute_divg, cur_divg)
    push!(metrics.relative_cost, rel_cost)

    metrics.cur_iter += 1
    metrics.cur_cost = cur_cost
    metrics.cur_divg = cur_divg
end

"""
    _build_conn(labels)

Helper function for `_connectivity_matrix`.
"""
function _build_conn(labels::Array{Int64,1})
    indices_1 = repeat(transpose(labels), length(labels), 1)
    indices_2 = repeat(labels, 1, length(labels))

    convert(Matrix{Float64}, indices_1 .== indices_2)
end

"""
    _connectivity_matrix(H)

Build connectivity matrix from H in which rows are latent factors and columns individuals.
"""
function _connectivity_matrix(H::Matrix{T}; clustering::Symbol=:max) where T <: Real
    if clustering == :max
        labels = argmax(H, dims=1)
        labels = (c -> c[1]).(labels)
        labels = vec(labels)
    elseif clustering == :kmeans
        throw(ArgumentError("clustering method not implemented"))
    end
    _build_conn(labels)
end

"""
    _check_stopping_crit(metrics, global_params)

Helper function for checking convergence in NMF algorithms.
"""
function _check_stopping_crit(metrics::NMFMetrics, global_params::AbstractNMFParams)
    stopping_crit_met = false

    if global_params.stopping_crit == :conn
        #### stop when the connectivity matrix remains constant over successive iterations.
        if metrics.cur_iter > global_params.stopping_iter
            range_iter = (metrics.cur_iter-global_params.stopping_iter):(metrics.cur_iter-1)
            stopping_crit_met = all(metrics.unchange_conn[range_iter])
        end

    elseif global_params.stopping_crit == :rel_cost
        #### stop when relative change in cost is below tolerance over successive iterations
        if metrics.cur_iter > global_params.stopping_iter
            range_iter = (metrics.cur_iter-global_params.stopping_iter):(metrics.cur_iter-1)
            stopping_crit_met = maximum(metrics.relative_cost[range_iter]) < global_params.stopping_tol
        end

    elseif global_params.stopping_crit == :cost
        #### stop when absolute change in cost is below tolerance over successive iterations
        if metrics.cur_iter > global_params.stopping_iter
            range_iter = (metrics.cur_iter-global_params.stopping_iter):(metrics.cur_iter-1)
            stopping_crit_met = maximum(metrics.absolute_cost[range_iter]) < global_params.stopping_tol
        end
    end

    stopping_crit_met || metrics.cur_iter == global_params.max_iter
end

"""
    _check_sanity_iter(W, H)

Helper function for checking iteration sanity in NMF algorithms.
"""

function _check_sanity_iter(W::Matrix{T}, H::Matrix{T}) where T <: Real
    sane = true
    cs_W = sum(W, dims=1)
    cs_H = sum(H, dims=1)
    rs_W = sum(W, dims=2)
    rs_H = sum(H, dims=2)

    if any(cs_W .== 0)
        @warn "Sanity check: One or more columns of W are null vectors, leading to division by zero in the update of H and W. Try reducing regularization if any, removing the scaling of W or modifying the initiazation method. The problem might as well be ill-posed."
        sane = false
    end

    if any(rs_W .== 0)
        @warn "Sanity check: One or more rows of W are null vectors, leading to division by zero in the update of H and W. Try reducing regularization if any, removing the scaling of W or modifying the initiazation method. The problem might as well be ill-posed."
        sane = false
    end

    if any(cs_H .== 0)
        @warn "Sanity check: One or more columns of H are null vectors, leading to division by zero in the update of H and W. Try reducing regularization if any, removing the scaling of W or modifying the initiazation method. The problem might as well be ill-posed."
        sane = false
    end

    if any(rs_H .== 0)
        @warn "Sanity check: One or more rows of H are null vectors, leading to division by zero in the update of H and W. Try reducing regularization if any, removing the scaling of W or modifying the initiazation method. The problem might as well be ill-posed."
        sane = false
    end

    sane
end


text:
# both energy and angular momentum
# assume time only coupled to radial coordinate
# need to think of a nice way to keep this efficient
# whilst allowing metrics with other couplings

"""
    E(m::AbstractMatrix{T}, v)
    E(m::AbstractMetricParams{T}, u, v)

Compute the energy for a numerically evaluated metric, and some velocity four vector `v`,
```math
E = - p_t = - g_{t\\nu} p^\\nu.
```

For null geodesics, the velocity is the momentum ``v^\\nu = p^\\nu``. For massive geodesics,
the mass ``\\mu`` needs to be known to compute ``\\mu v^\\nu = p^\\nu``.
"""
function E(metric::AbstractMatrix{T}, v) where {T}
    T(@inbounds -(metric[1, 1] * v[1] + metric[1, 4] * v[4]))
end
E(m::AbstractMetricParams{T}, u, v) where {T} = E(metric(m, u), v)

"""
    Lz(m::AbstractMatrix{T}, v)
    Lz(m::AbstractMetricParams{T}, u, v)

Compute the angular momentum for a numerically evaluated metric, and some velocity four vector `v`.
```math
L_z = p_\\phi = - g_{\\phi\\nu} p^\\nu.
```
"""
function Lz(metric::AbstractMatrix{T}, v) where {T}
    T(@inbounds metric[4, 4] * v[4] + metric[1, 4] * v[1])
end
Lz(m::AbstractMetricParams{T}, u, v) where {T} = Lz(metric(m, u), v)


text:
#
# Correctness Tests
#

require("test.jl")

using DataFrames

my_tests = ["test/dataarray.jl",
            "test/data.jl",
            "test/index.jl",
            "test/dataframe.jl",
            "test/operators.jl",
            "test/io.jl",
            # "test/formula.jl",
            "test/datastream.jl",
            "test/datamatrix.jl",
            "test/constructors.jl",
            "test/abstractarray.jl",
            "test/booleans.jl",
            "test/indexing.jl",
            "test/indexedvector.jl",
            "test/iteration.jl",
            "test/colfuncs.jl",
            "test/duplicates.jl",
            "test/padding.jl",
            "test/tabulation.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
vec_infeasible = @stochastic_model begin
    @stage 1 begin
        @parameters begin
            lb = zeros(2)
            c  = [3., 2.]
        end
        @decision(model, x[i in 1:2] >= lb[i])
        @objective(model, Min, dot(c, x))
    end
    @stage 2 begin
        @parameters begin
            q = [-15., -12.]
            T = [-1. 0;
                  0  -1.]
            W = [3. 2.;
                 2. 5.]
        end
        @uncertain ξ₁ ξ₂
        ub = [ξ₁, ξ₂]
        lb = 0.8 * ub
        @variable(model, lb[i] <= y[i in 1:2] <= ub[i])
        @objective(model, Min, dot(q, y))
        @constraint(model, T * x + W * y in MOI.Nonpositives(2))
    end
end

ξ₁ = Scenario(ξ₁ = 6., ξ₂ = 8., probability = 0.5)
ξ₂ = Scenario(ξ₁ = 4., ξ₂ = 4., probability = 0.5)

vec_infeasible_res = SPResult([27.2,41.6], 36.4, 9.2, 27.2, Inf, 9.2, Inf)
push!(problems, (vec_infeasible, [ξ₁,ξ₂], vec_infeasible_res, "Infeasible"))


text:
module Lift

export @^

include("liftmacro.jl")

end # module


text:
import Base: Float64, Float32, Float16, Int

abstract type AbstractLogFixPoint <: AbstractFloat end

""" LogFixPoint16 is a logarithmic fixed-point number format with

1 sign bit
7 signed integer bits, with its signbit (=bit 15) flipped (i.e. excess-8)
8 fraction bits.

Exceptions
0x0000      = zero
0x8000      = Not-a-Real (NaR)."""
primitive type LogFixPoint16 <: AbstractLogFixPoint 16 end

Base.iszero(x::LogFixPoint16) = reinterpret(UInt16,x) == 0x0000
Base.isnan(x::LogFixPoint16) = reinterpret(UInt16,x) == 0x8000
Base.signbit(x::LogFixPoint16) = (reinterpret(UInt16,x) & 0x8000) == 0x8000

Base.zero(::Type{LogFixPoint16}) = reinterpret(LogFixPoint16,0x0000)
nan(::Type{LogFixPoint16}) = reinterpret(LogFixPoint16,0x8000)

Base.floatmin(::Type{LogFixPoint16}) = reinterpret(LogFixPoint16,0x0001)
Base.floatmax(::Type{LogFixPoint16}) = reinterpret(LogFixPoint16,0x7fff)

# In the absence of -Inf,Inf define typemin,typemax as floatmin,floatmax
Base.typemin(::Type{LogFixPoint16}) = floatmin(LogFixPoint16)
Base.typemax(::Type{LogFixPoint16}) = floatmax(LogFixPoint16)

Base.one(::Type{LogFixPoint16}) = reinterpret(LogFixPoint16,0x4000)

function Base.:(-)(x::LogFixPoint16)
    iszero(x) && return zero(LogFixPoint16)
    isnan(x) && return nan(LogFixPoint16)
    return reinterpret(LogFixPoint16,reinterpret(UInt16,x) ⊻ 0x8000)
end

function Base.inv(x::LogFixPoint16)
    iszero(x) && return nan(LogFixPoint16)
    ui = reinterpret(UInt16,x)
    sign = ui & 0x8000
    val = -ui & 0x7fff		# two's complement of all except sign bit
    return reinterpret(LogFixPoint16,sign | val)
end

# Default number of fraction bits
const nfrac = Ref{Int}(9)
const nint = Ref{Int}(15-nfrac[])
const scale = Ref{Int}(2^nfrac[])
const scale_over_logof2 = Ref{Float32}(scale[]/log(2f0))
const max_nfrac_supported = 11
const rounding_mode = Ref{Symbol}(:lin)

""" 32-bit value of floatmax/min are derived from

    Int32(round(scale[]*(-(2^i-2^-j)-1))) + bit15flip

where i=nint-1, j=nfrac."""
const logfixpoint16_max = Int32(2^15-1)             # val of floatmax
const logfixpoint16_halfmin = Int32(-2^nfrac[]+1)   # val of half floatmin
const bit15flip = Int32(2^14)    # flip the meaning of the sign bit of integers
const bit14flip = UInt16(2^13)
const uintbit15flip = UInt16(2^14)    # used for power2

"""Constant to convert from round-to-nearest in log to linear space.

Let
    1 = s*log2(f1)
    1.5 = s*log2(f2)
    2 = s*log2(f3)

With some scale s. Now f2 is not the arithmetic mean of f1,f3. Let f2* be that mean

    f2* = 2^(1/s-1) + 2^(2/s-1)

Then we find a constant c_b to be added before rounding by

    1.5 = c_b + 2*log2(f2*)
    => c_b = 1.5 - s*log2(f2*)
           = 1/2 - s*(log2(2^(1/s) + 1) - 1)"""
function rounding_correction(s::Int)
    return 1/2 - s*(log2(2^(1/s)+1)-1)
end

const c_b = Ref{Float32}(Float32(rounding_correction(scale[])))

"""Convert a Float32 to LogFixPoint16 via the base-2 logarithm and rounding.
Round to nearest is applied in lin-space due to the addition of c_b."""
function LogFixPoint16(f::Float32)
    iszero(f) && return zero(LogFixPoint16)
    ~isfinite(f) && return nan(LogFixPoint16)

    # scale_over_logof2 = scale/log(2) to scale by nfrac fraction bits
    # use log instead of log2 as it's faster and include 1/log(2) in the constant
    val = Int32(round(c_b[] + scale_over_logof2[]*log(abs(f)))) + bit15flip
    sign = ((reinterpret(UInt32,f) & 0x8000_0000) >> 16) % UInt16

    sign,val = under_overflow_check(sign,val)
    result = sign | (val % UInt16)
    return reinterpret(LogFixPoint16,result)
end

"""Check for under or overflow:

    - Any result beyond floatmax returns floatmax. (No-overflow rounding mode)
    - Any result smaller than floatmin/2 returns 0, any result between
        floatmin and floatmin/2 returns floatmin."""
function under_overflow_check(sign::UInt16,val::Int32)

    underflow = val < logfixpoint16_halfmin
    overflow = val > logfixpoint16_max

    val, sign = underflow ? (zero(Int32), 0x0000) : (max(val,one(Int32)), sign)
    val = overflow ? logfixpoint16_max : val

    return sign,val
end

"""Pre-calcualte look-up table for conversion from LogFixPoint16 to Float32."""
function createF32LookupTable(nint::Int,nfrac::Int)

    N = 2^(nint+nfrac)      # length of the table
    table = Array{UInt32,1}(undef,N)
    s = 2^nfrac             # scale
    bias = 2^(nint+nfrac-1) # exponent bias (=bit15flip)

    for i in 1:N
        # convert index to signed integer, with signbit flipped
        si = i-bias
        # store float32s as UInt32 for bitwise operations
        table[i] = reinterpret(UInt32,Float32(2.0 ^ (si/s)))
    end

    return table
end

const f32lookup = createF32LookupTable(nint[],nfrac[])

"""Conversion function from LogFixPoint16 to Float32 via table lookup
after branching off the special cases 0,NaR."""
function Float32(x::LogFixPoint16)
    iszero(x) && return 0f0
    isnan(x) && return NaN32

    ui = reinterpret(UInt16,x)
    val = ui & 0x7fff     # set the signbit to zero

    # set all except sign bit to zero
    # sign = either 0x8000_0000 or 0x0000_0000
    sign = ((ui & 0x8000) % UInt32) << 16

    # combine sign and table lookup
    @inbounds f = sign | f32lookup[val]
    return reinterpret(Float32,f)
end

# conversion between LogFixPoint16 and various floats
Float64(x::LogFixPoint16) = Float64(Float32(x))
Float16(x::LogFixPoint16) = Float16(Float32(x))
LogFixPoint16(x::Float64) = LogFixPoint16(Float32(x))
LogFixPoint16(x::Float16) = LogFixPoint16(Float32(x))
Int(x::LogFixPoint16) = Int(Float32(x))
LogFixPoint16(x::Int) = LogFixPoint16(Float32(x))

"""Multiplication for LogFixPoint16, equivalent to an addition of the exponents."""
function Base.:(*)(x::LogFixPoint16,y::LogFixPoint16)
    # special cases NaR*y = x*NaR = NaR
    isnan(x) | isnan(y) && return nan(LogFixPoint16)
    # if no NaN present: 0*y = x*0 = 0
    iszero(x) | iszero(y) && return zero(LogFixPoint16)

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    # mask exponent
    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    # mask the sign and cast to Int32
    xval = (uix & 0x7fff) % Int32
    yval = (uiy & 0x7fff) % Int32

    # resulting sign
    sign = xsign ⊻ ysign

    # ADD EXPONENTS
    # xval and yval both contain the bias bit15flip
    # so subtract bit15flip, such that the result is again
    # biased with (a single) + bit15flip
    val = xval + yval - bit15flip

    # check for over or underflow
    sign,val = under_overflow_check(sign,val)

    # merge sign and exponent back together
    result = sign | (val % UInt16)
    return reinterpret(LogFixPoint16,result)
end

"""Division for LogFixPoint16, equivalent to a subtraction of the exponents."""
function Base.:(/)(x::LogFixPoint16,y::LogFixPoint16)
    # special case, if either x,y =NaN, or y=0 return NaR
    isnan(x) | isnan(y) | iszero(y) && return nan(LogFixPoint16)
    # else if x=0 return 0
    iszero(x) && return zero(LogFixPoint16)

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    # extract sign
    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    # extract exponent
    xval = (uix & 0x7fff) % Int32
    yval = (uiy & 0x7fff) % Int32

    # combine to obtain the sign of the result
    sign = xsign ⊻ ysign

    # SUBTRACT EXPONENTS
    # xval and yval both contain the bias bit15flip
    # so add bit15flip, such that the result is again
    # biased with (a single) + bit15flip
    val = xval - yval + bit15flip

    sign,val = under_overflow_check(sign,val)

    # merge sign and exponent val
    result = sign | (val % UInt16)
    return reinterpret(LogFixPoint16,result)
end

function Base.sqrt(x::LogFixPoint16)
    iszero(x) && return zero(LogFixPoint16)
    signbit(x) && return nan(LogFixPoint16) #TODO throw DomainError?

    uix = reinterpret(UInt16,x)
    #TODO is this correct rounding?
    uix = (uix >> 1) + bit14flip
    return reinterpret(LogFixPoint16,uix)
end

function power2(x::LogFixPoint16)
    iszero(x) && return zero(LogFixPoint16)
    isnan(x) && return nan(LogFixPoint16)
    uix = reinterpret(UInt16,x)
    uix = (uix << 1) + uintbit15flip
    return reinterpret(LogFixPoint16,uix & 0x7fff)
end

""" Precomputes the Gaussian Logarithm as a table lookup for addition. Let

X = (-1)^signx * 2^x
Y = (-1)^signy * 2^y
Z = (-1)^signz * 2^z

then additions can be computed via the Gaussian logarithms

z = log2(||X|+|Y||) = x + log2(1+2^(y-x))

For fixed-point number exponents, interpreted as integers
including a bias c, x̂ = ax + c, where a = 2^n_frac the number of fraction bits:

ẑ = x̂ + a*log2(1+2^((ŷ-x̂)/a)     or
ẑ = ŷ - (ŷ-x̂) + a*log2(1+2^((ŷ-x̂)/a)

The last two terms are precomputed into a table lookup."""
function createAddLookup(scale::Int)

    tab = Array{Int32,1}(undef,max_table_size)

    for i in 0:max_table_size-1
        tab[i+1] = -i + Int(round(c_b[] + scale*log2(1+2^(i/scale))))
    end

    return tab
end

""" Precomputes the Gaussian Logarithm as a table lookup for subtraction.
Subtractions can be computed via the Gaussian logarithms

z = log2(||X|-|Y||) = x + log2(abs(1-2^(y-x)))

For fixed-point number exponents, interpreted as integers
including a bias c, x̂ = ax + c, where a = 2^n_frac the number of fraction bits:

ẑ = x̂ + a*log2(abs(1-2^((ŷ-x̂)/a))     or
ẑ = ŷ - (ŷ-x̂) + a*log2(abs(1-2^((ŷ-x̂)/a))

The last two terms are precomputed into a table lookup."""
function createSubLookup(scale::Int)
    tab = Array{Int32,1}(undef,max_table_size)

    # set the first entry manually to avoid -Inf due to the log2(0)
    tab[1] = -logfixpoint16_max-1

    for i in 1:max_table_size-1
        tab[i+1] = -i + Int(round(c_b[] + scale*log2(abs(1-2^(i/scale)))))
    end

    return tab
end

# table indices higher than that are 0
function find_max_diff_res(scale::Int)
    i = 1
    while (-i + round(c_b[] + scale*log2(1+2^(i/scale)))) > 0.0
        i += 1
    end
    return Int32(i)
end

const max_diff_resolvable = Ref{Int32}(find_max_diff_res(scale[]))
const max_table_size = find_max_diff_res(2^max_nfrac_supported)+1
const addTable = createAddLookup(scale[])
const subTable = createSubLookup(scale[])

function Base.:(+)(x::LogFixPoint16,y::LogFixPoint16)
    isnan(x) | isnan(y) && return nan(LogFixPoint16)

    # zero is an exceptional case and not reprsentable as (-1)^s * 2^x
    iszero(x) && return y
    iszero(y) && return x

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    xval = (uix & 0x7fff) % Int32
    yval = (uiy & 0x7fff) % Int32

    # y is always the larger value
    # resulting sign is always that of y
    # in a-a sign is 1, but case is caught in underflow: sign set to 0
    xval,yval,xsign,ysign = xval > yval ? (yval,xval,ysign,xsign) : (xval,yval,xsign,ysign)
    diff = yval - xval     # diff >= 0

    # pull Gaussian logarithms from addTable, subTable
    # update yval with the increment
    @inbounds increment = diff < max_diff_resolvable[] ?
        (xsign == ysign ? addTable[diff+1] : subTable[diff+1]) : zero(Int32)
    yval += increment

    # check for over or underflow
    ysign,yval = under_overflow_check(ysign,yval)

    result = ysign | (yval % UInt16)
    return reinterpret(LogFixPoint16,result)
end

function diff_val(x::LogFixPoint16,y::LogFixPoint16)
    isnan(x) | isnan(y) && return nan(LogFixPoint16)

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    xval = (uix & 0x7fff) % Int32
    yval = (uiy & 0x7fff) % Int32

    # y is always the larger value
    # resulting sign is always that of y
    # in a-a sign is 1, but case is caught in underflow: sign set to 0
    xval,yval,xsign,ysign = xval > yval ? (yval,xval,ysign,xsign) : (xval,yval,xsign,ysign)
    diff = yval - xval     # diff >= 0
end

"""Subtraction for LogFixPoint16 via sign switch of y as in x-y = x + (-y)."""
Base.:(-)(x::LogFixPoint16,y::LogFixPoint16) = x + (-y)

"""Rounding to integer via conversion to Float32."""
Base.round(x::LogFixPoint16, r::RoundingMode{:Up}) = Int(ceil(Float32(x)))
Base.round(x::LogFixPoint16, r::RoundingMode{:Down}) = Int(floor(Float32(x)))
Base.round(x::LogFixPoint16, r::RoundingMode{:Nearest}) = Int(round(Float32(x)))

for t in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
	@eval Base.promote_rule(::Type{LogFixPoint16}, ::Type{$t}) = LogFixPoint16
end

function Base.nextfloat(x::LogFixPoint16)
    isnan(x) && return nan(LogFixPoint16)
    (x == -floatmin(LogFixPoint16)) && return zero(LogFixPoint16)

    sign = signbit(x)
    ui = reinterpret(UInt16,x)
    sign && return reinterpret(LogFixPoint16,ui - 0x0001)
    return reinterpret(LogFixPoint16,ui + 0x0001)
end

function Base.nextfloat(x::LogFixPoint16,n::Int)
    for i in 1:n
        x = nextfloat(x)
    end
    return x
end

function Base.prevfloat(x::LogFixPoint16)
	# -floatmax is 0xffff and would otherwise prevfloat to 0=0x0000
    isnan(x) | (x == -floatmax(LogFixPoint16)) && return nan(LogFixPoint16)

    sign = signbit(x)
    ui = reinterpret(UInt16,x)
    sign && return reinterpret(LogFixPoint16,ui + 0x0001)
    return reinterpret(LogFixPoint16,ui - 0x0001)
end

function Base.prevfloat(x::LogFixPoint16,n::Int)
    for i in 1:n
        x = prevfloat(x)
    end
    return x
end

function Base.log2(x::LogFixPoint16)
	signbit(x) | iszero(x) && return nan(LogFixPoint16)
	ui =  reinterpret(UInt16,x) % Int
	ui -= Int(bit15flip)
	ui /= scale
	return LogFixPoint16(ui)
end

function Base.:(==)(x::LogFixPoint16,y::LogFixPoint16)
    isnan(x) | isnan(y) && return false
    return reinterpret(UInt16,x) == reinterpret(UInt16,y)
end

function Base.:(>)(x::LogFixPoint16,y::LogFixPoint16)
    isnan(x) | isnan(y) && return nan(LogFixPoint16)

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    xsign > ysign && return false
    xsign < ysign && return true

    return uix > uiy
end

function Base.:(<)(x::LogFixPoint16,y::LogFixPoint16)
    isnan(x) | isnan(y) && return nan(LogFixPoint16)

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    xsign > ysign && return true
    xsign < ysign && return false

    return uix < uiy
end

function Base.:(<=)(x::LogFixPoint16,y::LogFixPoint16)
    isnan(x) | isnan(y) && return nan(LogFixPoint16)

    uix = reinterpret(UInt16,x)
    uiy = reinterpret(UInt16,y)

    xsign = uix & 0x8000
    ysign = uiy & 0x8000

    xsign > ysign && return true
    xsign < ysign && return false

    return uix <= uiy
end

# Showing
function Base.show(io::IO, x::LogFixPoint16)
    if isnan(x)
        print(io, "NaR")
    else
		io2 = IOBuffer()
        print(io2,Float32(x))
        f = String(take!(io2))
        print(io,"LogFixPoint16("*f*")")
    end
end

Base.bitstring(x::LogFixPoint16) = bitstring(reinterpret(UInt16,x))

function Base.bitstring(x::LogFixPoint16,mode::Symbol)
    if mode == :split	# split into sign, integer, fraction
        s = bitstring(x)
		return "$(s[1]) $(s[2:nint[]+1]) $(s[nint[]+2:end])"
    else
        return bitstring(x)
    end
end


text:
# ### Box plot (2)

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

## This outputs BoundsError. This is a bug of StatPlots.jl.
## Boxplot and violin broken - https://github.com/JuliaPlots/StatPlots.jl/issues/198
using StatsPlots, Statistics, Random
Random.seed!(2018)

n = 30
science = rand(1:10, n)
english = rand(1:10, n)

@show median(science)
@show median(english)

boxplot([science english], label=["science" "english"])


text:
"""
VMeasure is type containing the group angle (in radians), group velocity (in m/s, or km/s), error on group angle (in radians), error on log(v).
"""
struct VMeasure
    angle::Real    # group angle (radians)
    value::Real  # group velocity
    err_angle::Real   # error on group angle
    err_value::Real  # error on log(v)
end

"""
    VMeasure(angle::Vector, value::Vector, err_angle::Vector, err_value::Vector)

Return an array of VMeasure, element-wise.
"""
function VMeasure(angle::Vector, value::Vector, err_angle::Vector, err_value::Vector)
    return [VMeasure(angle[i], value[i], err_angle[i], err_value[i]) for i in eachindex(angle)]
end

"""
    VMeasure(angle::Vector, value::Vector, err_angle::Real, err_value::Real)

Return an array of VMeasure, element-wise, using the same errors on angle and log(v).
"""
function VMeasure(angle::Vector, value::Vector, err_angle::Real, err_value::Real)
    return [VMeasure(angle[i], value[i], err_angle, err_value) for i in eachindex(angle)]
end


text:
recA{T}(::Type{T},::Chebyshev,k) = 2one(T)
recB{T}(::Type{T},::Chebyshev,::) = zero(T)
recC{T}(::Type{T},::Chebyshev,k) = one(T)   # one(T) ensures we get correct type

recα{T}(::Type{T},::Chebyshev,::) = zero(T)
recβ{T}(::Type{T},::Chebyshev,k) = ifelse(k==1,one(T),one(T)/2)   # one(T) ensures we get correct type,ifelse ensures inlining
recγ{T}(::Type{T},::Chebyshev,k) = one(T)/2   # one(T) ensures we get correct type

## Evaluation

function evaluatechebyshev{T<:Number}(n::Integer,x::T)
    if n == 1
        [one(T)]
    else
        p = zeros(T,n)
        p[1] = one(T)
        p[2] = x

        for j=2:n-1
            p[j+1] = 2x*p[j] - p[j-1]
        end

        p
    end
end

function getindex{DD<:Interval}(op::ConcreteEvaluation{Chebyshev{DD},Bool},j::Integer)
    T=eltype(op)
    if op.order == 0
        ifelse(op.x || isodd(j),  # right rule
            one(T),
            -one(T))
    else
        #TODO: Fast version
        op[j:j][1]
    end
end

function getindex{DD<:Interval}(op::ConcreteEvaluation{Chebyshev{DD},Bool},k::Range)
    T=eltype(op)
    x = op.x
    d = domain(op)
    p = op.order
    cst = T((2/(d.b-d.a))^p)
    n=length(k)

    if x
        ret = ones(T,n)
    else
        ret = Array(T,n)
        k1=1-first(k)
        @simd for j=k
            @inbounds ret[j+k1]=(-1)^(p+1)*(-one(T))^j
        end
    end

    for m=0:p-1
        k1=1-first(k)
        @simd for j=k
            @inbounds ret[j+k1] *= (j-1)^2-m^2
        end
        scal!(T(1/(2m+1)), ret)
    end

    scal!(cst,ret)
end

function getindex{DD<:Interval,M<:Real}(op::ConcreteEvaluation{Chebyshev{DD},M},j::Integer)
    if op.order == 0
        evaluatechebyshev(j,tocanonical(domain(op),op.x))[end]
    else
        error("Only zero–second order implemented")
    end
end

function getindex{DD<:Interval,M<:Real}(op::ConcreteEvaluation{Chebyshev{DD},M},k::Range)
    if op.order == 0
        evaluatechebyshev(k[end],tocanonical(domain(op),op.x))[k]
    else
        error("Only zero–second order implemented")
    end
end

# Multiplication

Base.stride{U<:Chebyshev,V<:Chebyshev}(M::ConcreteMultiplication{U,V}) =
    stride(M.f)

function chebmult_getindex(cfs::Vector,k::Integer,j::Integer)
    n=length(cfs)

    ret=zero(eltype(cfs))

    # Toeplitz part
    if k == j
        ret += cfs[1]
    elseif k > j && k-j+1 ≤ n
        ret += cfs[k-j+1]/2
    elseif k < j && j-k+1 ≤ n
        ret += cfs[j-k+1]/2
    end

    # Hankel part
    if k ≥ 2 && k+j-1 ≤ n
        ret += cfs[k+j-1]/2
    end

    ret
end

getindex{T,C<:Chebyshev}(M::ConcreteMultiplication{C,C,T},k::Integer,j::Integer) =
    chebmult_getindex(coefficients(M.f),k,j)

getindex{PS<:PolynomialSpace,T,C<:Chebyshev}(M::ConcreteMultiplication{C,PS,T},k::Integer,j::Integer) =
    M[k:k,j:j][1,1]

function Base.convert{C<:Chebyshev,V,T}(::Type{BandedMatrix},S::SubOperator{T,ConcreteMultiplication{C,C,V,T},Tuple{UnitRange{Int},UnitRange{Int}}})
    ret=bzeros(S)

    kr,jr=parentindexes(S)
    cfs=parent(S).f.coefficients

    # Toeplitz part
    sym_toeplitz_axpy!(1.0,0.5,cfs,kr,jr,ret)

    #Hankel part
    hankel_axpy!(0.5,cfs,kr,jr,ret)

    # divide first row by half
    if first(kr)==1
        if first(jr)==1
            ret[1,1]+=0.5cfs[1]
        end

        for j=1:min(1+ret.u,size(ret,2))
            ret[1,j]/=2
        end
    end

    ret
end

function Base.convert{PS<:PolynomialSpace,V,T,C<:Chebyshev}(::Type{BandedMatrix},S::SubOperator{T,ConcreteMultiplication{C,PS,V,T},
                                                                            Tuple{UnitRange{Int},UnitRange{Int}}})
    M=parent(S)
    kr,jr=parentindexes(S)

    A=bzeros(S)

    a=coefficients(M.f)

    shft=bandshift(A)

    for k=kr ∩ jr
        A[k-kr[1]+1,k-jr[1]+1]=a[1]
    end

    if length(a) > 1
        sp=M.space
        jkr=max(1,min(kr[1],jr[1])-length(a)+1):max(kr[end],jr[end])+length(a)-1

        #Multiplication is transpose
        J=Recurrence(sp)[jkr,jkr]
        C1=J

        # the sub ranges of jkr that correspond to kr, jr
        kr2,jr2=kr-jkr[1]+1,jr-jkr[1]+1

        BLAS.axpy!(a[2],@compat(view(C1,kr2,jr2)),A)
        C0=beye(size(J,1),size(J,2),0,0)

        for k=1:length(a)-2
            C1,C0=2J*C1-C0,C1
            BLAS.axpy!(a[k+2],@compat(view(C1,kr2,jr2)),A)
        end
    end

    A
end

## Derivative

Derivative{DD<:Interval}(sp::Chebyshev{DD},order::Integer) =
    ConcreteDerivative(sp,order)

rangespace{DD<:Interval}(D::ConcreteDerivative{Chebyshev{DD}}) =
    Ultraspherical(D.order,domain(D))
bandinds{DD<:Interval}(D::ConcreteDerivative{Chebyshev{DD}}) = 0,D.order
Base.stride{DD<:Interval}(D::ConcreteDerivative{Chebyshev{DD}}) = D.order

function getindex{DD<:Interval,K,T}(D::ConcreteDerivative{Chebyshev{DD},K,T},k::Integer,j::Integer)
    m=D.order
    d=domain(D)

    if j==k+m
        C=.5pochhammer(1.,m-1)*(4./(d.b-d.a)).^m
        (C*(m+k-one(T)))::T
    else
        zero(T)
    end
end

linesum{DD<:Interval}(f::Fun{Chebyshev{DD}}) =
    sum(setcanonicaldomain(f))*arclength(d)/2

## Clenshaw-Curtis functional

for (Func,Len) in ((:DefiniteIntegral,:complexlength),(:DefiniteLineIntegral,:arclength))
   ConcFunc = parse("Concrete"*string(Func))
   @eval begin
       function getindex{D<:Interval,T}(Σ::$ConcFunc{Chebyshev{D},T},k::Integer)
           d = domain(Σ)
           C = $Len(d)/2

           isodd(k) ? 2C/(k*(2-k)) : zero(T)
       end
       function getindex{D<:Interval,T}(Σ::$ConcFunc{Chebyshev{D},T},kr::Range)
           d = domain(Σ)
           C = $Len(d)/2

           promote_type(T,typeof(C))[isodd(k) ? 2C/(k*(2-k)) : zero(T) for k in kr]
       end
   end
end


text:
using TreeKnit
TreeKnit.comonicon_install()

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
struct BGGUser
    id::Int
    name::String
    country::String
    yearregistered::Int
    lastlogin::Date
end

"""
    userinfo(name::String)
    userinfo(names::AbstractArray{<:String})

Get user from name string.
"""
function userinfo(name::AbstractString; waittime=0.0f0)
    doc = get_xml("$XMLAPI2/user?name=$(name)")
    id = parse(Int, findfirst("/user", doc)["id"])
    name = findfirst("/user", doc)["name"]
    yearregistered = parse(Int, findfirst("/user/yearregistered", doc)["value"])
    lastlogin = Date(findfirst("/user/lastlogin", doc)["value"])
    country = findfirst("/user/country", doc)["value"]

    sleep(waittime)
    return BGGUser(id, name, country, yearregistered, lastlogin)
end

function userinfo(names::AbstractArray{<:String}; waittime=2.0f0) # add waittime by default
    usercount = length(names)
    users = Vector{BGGUser}(undef, usercount) # pre-allocate
    @showprogress 1 "Scraping $usercount user profiles..." for (i, n) in enumerate(names)
        users[i] = userinfo(n; waittime=waittime)
    end
    return users
end

"""
    buddies(name::String)
    buddies(names::AbstractArray{<:String})
    buddies(user::BGGUser)

Return list of user names of buddies of a user.
"""
function buddies(name::String; waittime=0.0f0)
    doc = get_xml("$XMLAPI2/user?name=$(name)&buddies=1")

    sleep(waittime)
    return [b["name"] for b in findall("/user/buddies/buddy", doc)]
end

function buddies(names::AbstractArray{<:String}; waittime=2.0f0) # add waittime by default
    return [buddies(n; waittime=waittime) for n in names]
end

buddies(user::BGGUser) = buddies(user.name)

"""
    userreviews(name::String)
    userreviews(names::AbstractArray{<:String})
    userreviews(user::BGGUser)

Return list of all board game reviews a user wrote (ignoring expansions).
"""
function userreviews(name::String; waittime=0.0f0)
    doc = get_xml(
        "$(XMLAPI2)/collection?username=$(name)&rated=1&stats=1&excludesubtype=boardgameexpansion",
    )
    games = findall("/items/item", doc)

    sleep(waittime)
    return _parse_review_from_collection.(games, name)
end

function userreviews(names::AbstractArray{<:String}; waittime=2.0f0) # add waittime by default
    return collect(Iterators.flatten(userreviews(n; waittime=waittime) for n in names))
end

userreviews(user::BGGUser) = userreviews(user.name)

function _parse_review_from_collection(n::EzXML.Node, username)
    id = parse(Int, n["objectid"])

    gamename = findfirst("name", n).content
    rating = parse(Float32, findfirst("stats/rating", n)["value"])
    _comment = findfirst("comment", n)
    comment = ""
    if !isnothing(_comment)
        comment = _comment.content
    end
    return BGGReview(id, gamename, username, rating, comment)
end

"""
    collection(name::String)
    collection(user::BGGUser)

Return game ids in collection of user.
"""
collection(name::String) = getproperty.(userreviews(name), :id)
collection(user::BGGUser) = collection(user.name)


text:
int i = 1
ans = []
@async function f1()
  i = 1 + i%3
  ans.push(i)
end


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))
using JuliaAcademyData; activate("Foundations of machine learning")

# # Representing data in a computer
#
# The core of data science and machine learning is **data**: we are interested in extracting knowledge from data.
#
# But how exactly do computers represent data? Let's find out exactly what an "artificial intelligence" has at its disposal to learn from.

#-

# ## Data is represented as arrays
#
# Let's take a look at some fruit. Using the `Images.jl` library, we can load in some images:

using Images

apple = load(datapath("data/10_100.jpg"))

#-

banana = load(datapath("data/104_100.jpg"))

# Here we have images of apples and bananas. We would eventually like to build a program that can automatically distinguish between the two. However, the computer doesn't "see" an apple or a banana; instead, it just sees numbers.
#
# An image is encoded in something called an **array**, which is like a container that has boxes or slots for individual pieces of data:
#
# An array is a bunch of numbers in connected boxes; the figure above shows a 1-dimensional array. Our images are instead 2-dimensional arrays, or matrices, of numbers, arranged something like this:

#-

# <img src="https://raw.githubusercontent.com/JuliaComputing/JuliaAcademyData.jl/master/courses/Foundations%20of%20machine%20learning/data/array2d.png" alt="attachment:array2d.png" width="500"/>

#-

# For example, `apple` is an image, consisting of a 100x100 array of numbers:

typeof(apple)

#-

size(apple)

#-

a = [ 1 2 3;4 5 6]

#-

typeof(a)

#-

size(a)

# We can grab the datum stored in the box at row `i` and column `j` by *indexing* using square brackets: `[i, j]`. For example, let's get the pixel (piece of the image) in box $(40, 60)$, i.e. in the 40th row and 60th column of the image:

apple

#-

dump(apple[40, 60])

#-

apple[18:20,29:31]

# We see that Julia displays a coloured box! Julia, via the `Colors.jl` package, is clever enough to display colours in a way that is useful to us humans!
#
# So, in fact, an image is a 2D array, in which each element of the array is an object (a collection of numbers) describing a coloured pixel.

#-

# ## Colors as numbers
#
# How, then, are these colors actually stored? Computers store colors in RGB format, that is they store a value between 0 and 1 for each of three "channels": red, green, and blue. Here, 0 means none of that color and 1 means the brightest form of that color. The overall color is a combination of those three colors.
#
# For example, we can pull out the `red` value using the function `red` applied to the color. Since internally the actual value is stored in a special format, we choose to convert it to a standard floating-point number using the `Float64` function:

Float64(red(apple[40, 60]))

#-

using Statistics

#-

[ mean(float.(c.(img))) for c = [red,green,blue], img = [apple,banana] ]

#-

using Plots
histogram(float.(green.(apple[:])),color="red",label="apple", normalize=true, nbins=25)
histogram!(float.(green.(banana[:])),color="yellow",label="banana",normalize=true, nbins=25)

#-

apple

#-

float(red(banana[50,20]))

#-

banana[50,20]

#-

pixel = apple[40, 60]

red_value   = Float64( red(pixel) )
green_value = Float64( green(pixel) )
blue_value  = Float64( blue(pixel) )

print("The RGB values are ($red_value, $green_value, $blue_value)")

# Since the red value is high while the others are low, this means that at pixel `(40, 60)`, the picture of the apple is very red. If we do the same at one of the corners of the picture, we get the following:

pixel = apple[1, 1]

red_value   = Float64( red(pixel) )
green_value = Float64( green(pixel) )
blue_value  = Float64( blue(pixel) )

print("The RGB values are ($red_value, $green_value, $blue_value)")

#-

apple

# We see that every color is bright, which corresponds to white.

#-

# ## Working on an image as a whole
#
# In Julia, to apply a function to the whole of an array, we place a `.` between the function name and the left parenthesis (`(`), so the following gives us the `red` value of every pixel in the image:

redpartofapple = Float64.(red.(apple))
mean(redpartofapple)

#-

using Plots

#-

gr()

#-

histogram(redpartofapple[:],color=:red,label="redness in the apple")

# Note that we get a 2D array (matrix) back.

#-

# Julia's [mathematical standard library](https://docs.julialang.org/en/stable/stdlib/math/#Mathematics-1) has many mathematical functions built in. One of them is the `mean` function, which computes the average value. If we apply this to our apple:

mean(Float64.(red.(apple)))

# we see that the value indicates that the average amount of red in the image is a value between the amount of red in the apple and the amount of red in the white background.
#
# *Somehow we need to teach a computer to use this information about a picture to recognize that there's an apple there!*

#-

# ## A quick riddle
#
# Here's a quick riddle. Let's check the average value of red in the image of the banana.

mean(Float64.(red.(banana)))

# Oh no, that's more red than our apple? This isn't a mistake and is actually true! Before you move onto the next exercise, examine the images of the apple and the banana very carefully and see if you can explain why this is expected.

#-

# #### Exercise 1
#
# What is the average value of blue in the banana?
#
# (To open a new box use <ESC>+b (b is for "below", what do you think a does?))

#-

# #### Solution
#
# We can calculate the average value of blue in the banana via

mean(Float64.(blue.(banana)))

# The result is approximately `0.8`.

#-

# #### Exercise 2
#
# Does the banana have more blue or more green?

#-

# #### Solution
#
# The average value of green in the banana is

mean(Float64.(green.(banana)))

# which gives approximately `0.88`. The banana has more green on average.

max_stars_repo_path: Courses/Foundations of machine learning/0100.Representing-data-in-a-computer.jl

max_stars_repo_name: fercarozzi/JuliaAcademyMaterials

max_issues_repo_name: fercarozzi/JuliaAcademyMaterials

max_forks_repo_name: fercarozzi/JuliaAcademyMaterials


text:
function fermat(phifunc_xy, forcefunc_xy, rs, tgα, N=100; reltol=1e-21, abstol=1e-21, ext=1.0)
    const C2 = C*C
    function f(x,yvec)
        (y, v) = yvec
        y_prime = v
        v_prime = 0.0

        phi = phifunc_xy(x,y)/C/C
        dphidx, dphidy = forcefunc_xy(x,y)
        dphidx /= C2
        dphidy /= C2
        v2 = v*v

        v_prime = 2*(dphidx*(v2*v+v)-dphidy*(v2+1)^2) / (1-2*phi)
        [y_prime, v_prime]
    end
    _x = linspace(0.0,1,div(N,2)).^2
    _l = [-reverse(_x)*rs ; _x*rs.*ext]
    x, u = ode45(f, [0.0; tgα], _l; reltol=reltol, abstol=abstol)
    y = map(u -> u[1], u);
    dy = map(u -> u[2], u);
    x,y
end

function fermat_y(forcefunc_y, rs, tgα, N=100; reltol=1e-21, abstol=1e-21, ext=1.0)
    const C2 = C*C
    function f(x,yvec)
        (y, v) = yvec
        y_prime = v
        v_prime = 0.0

        v_prime = -2*forcefunc_y(x,y)/C2
        [y_prime, v_prime]
    end
    _x = linspace(0.0,1,div(N,2)).^2
    _l = [-reverse(_x)*rs ; _x*rs.*ext]
    x, u = ode45(f, [0.0; tgα], _l; reltol=reltol, abstol=abstol)
    y = map(u -> u[1], u);
    dy = map(u -> u[2], u);
    x,y
end

max_stars_repo_name: skariel/Lensing

max_issues_repo_name: skariel/Lensing

max_forks_repo_name: skariel/Lensing


text:
module MUC16Analysis

include("distance_matrix.jl")
include("letter_list.jl")
include("subset_fasta.jl")
include("upgma.jl")

export sub_mat_ali
export dist_mat_gen
export letter_list
export subset_sequence
export add_missing
export minimum_cell
export combine_strings
export combine_strings_dist
export combine_df
export upgma
export create_dm_df
export fasta_to_text

end #module


text:
# Author: Ashwani Rathee
# Program to find quotient and remainder numbers

println("Enter Dividend and Divisor")
dividend=readline()
dividend=parse(Int64,dividend)
divisor=readline()
divisor=parse(Int64,divisor)
quotient=dividend/divisor
#remainder=dividend%divisor
remainder=mod(dividend,divisor)
println("The quotient is :", quotient)
println("The remainder is: ", remainder)

max_stars_repo_name: ashwani-rathee/Learning--Julia

max_issues_repo_name: ashwani-rathee/Learning--Julia

max_forks_repo_name: ashwani-rathee/Learning--Julia


text:
@testset "test_parallel_Vec.jl" begin
   if MPI.Comm_size(MPI.COMM_WORLD) == 2
        # Test serialization
        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            vec = PetscVec()
            setSize!(vec, n_local=(Int32)(3))

            vec[1:3] = [i for i in 1:3]

            assemble!(vec)

            @test serializeToZero(vec) == (Float64)[1., 2., 3., 4., 5., 6.]
        else
            vec = PetscVec()
            setSize!(vec, n_local=(Int32)(3))

            vec[4:6] = [i for i in 4:6]

            assemble!(vec)
            serializeToZero(vec)
        end

        # Test getArray and restoreArray
        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            vec = PetscVec()
            setSize!(vec, n_local=(Int32)(3))

            vec[1:3] = [i for i in 1:3]

            assemble!(vec)

            array = MiniPETSc._getArray(vec)
            MiniPETSc._restoreArray(vec, array)

            @test vec[1:3] == [i for i in 1:3]
        else
            vec = PetscVec()
            setSize!(vec, n_local=(Int32)(3))

            vec[4:6] = [i for i in 4:6]

            assemble!(vec)

            array = MiniPETSc._getArray(vec)
            MiniPETSc._restoreArray(vec, array)

            @test vec[4:6] == [i for i in 4:6]
        end
    end
end


text:
function de_bruijn_ize(str::ASCIIString, k::Integer)
    edges = Array{Tuple{ASCIIString,ASCIIString},1}()
    nodes = Set{ASCIIString}()
    for i in 1:(length(str)-k)
        push!(edges, (str[i:i+k-1],str[i+1:i+k]))
        push!(nodes, str[i:i+k-1], str[i+1:i+k])
    end
    return edges,nodes
end

type Node
    node_name::ASCIIString
    num_in::Int64
    num_out::Int64

    function Node(name::ASCIIString)
        new(name,0,0)
    end
end

isSemiBalanced(nd::Node) = abs(nd.num_in - nd.num_out) == 1
isBalanced(nd::Node) = nd.num_in == nd.num_out
name(nd::Node) = nd.node_name

pointer_from_node(nd::Node) = convert(Ptr{Node}, pointer_from_objref(nd))

function chop(str::ASCIIString, k::Int64)
    num_kmers = length(str)-k+1
    mers = Array{Tuple{ASCIIString,ASCIIString,ASCIIString},1}()
    for i in 1:num_kmers
        kmer = str[i:i+k-1]
        push!(mers,(kmer,kmer[1:end-1],kmer[2:end]))
    end
    mers
end

type DebruijnGraph
    G::Dict{Node,Array{Node,1}}
    nodes::Dict{ASCIIString,Node} # map k-1 mers to Node objects

    nsemi::Int64
    nbal::Int64
    nneither::Int64
    head::Ptr{Node}
    tail::Ptr{Node}

    function DebruijnGraph(strIter,k,circularize=false)
        G = Dict{Node,Array{Node,1}}()
        nodes = Dict{ASCIIString,Node}()
        for str in strIter
            if circularize
                str = string(str, str[1:k-1])
            end
            for (_,km1L,km1R) in chop(str,k)
                if !in(km1L, keys(nodes))
                    nodes[km1L] = Node(km1L)
                end
                if !in(km1R, keys(nodes))
                    nodes[km1R] = Node(km1R)
                end
                nodeL = nodes[km1L]
                nodeR = nodes[km1R]

                nodeL.num_in += 1
                nodeL.num_out += 1

                if !in(nodeL,keys(G))
                    G[nodeL] = Node[nodeR]
                else
                    push!(G[nodeL],nodeR)
                end
            end
        end
        nsemi,nbal,nneither = 0,0,0
        head,tail = Ptr{Node}(0),Ptr{Node}(0)

        for nd in values(nodes)
            if isBalanced(nd)
                nbal += 1
            elseif isSemiBalanced()
                if nd.num_in == nd.num_out + 1
                    tail = pointer_from_node(nd)
                else
                    head = pointer_from_node(nd)
                end
                nsemi += 1
            else
                nneither += 1
            end
        end

        new(G,nodes,nsemi,nbal,nneither,head,tail)
    end

end

nnodes(g::DebruijnGraph) = length(g.nodes)
nedges(g::DebruijnGraph) = length(g.G)
hasEulerianWalk(g::DebruijnGraph) = g.nneither ==0 && g.nsemi == 2
hasEulerianCycle(g::DebruijnGraph) = g.nneither ==0 && g.nsemi == 0
isEulerian(g::DebruijnGraph) = hasEulerianWalk(g) || hasEulerianCycle(g)

function eulerianWalkOrCyle(g::DebruijnGraph)
    @assert isEulerian(g) == true

    G = g.G
    if hasEulerianWalk(g)
       if !in(g.head, G[unsafe_load(g.tail)])
           G[unsafe_load(g.tail)] = Node[unsafe_load(g.head)]
       else
           push!(G[unsafe_load(g.tail)],unsafe_load(g.head))
       end
    end

    tour = Node[]
    src = collect(keys(G))[1]

    function _visit(nd)
        while length(G[nd])>0
            dst = pop!(G[nd])
            _visit(dst)
        end
        push!(tour,nd)
    end

    _visit(src)
    tour =  reverse(tour)[1:end-1]

    if hasEulerianWalk(G)
        ind = findfirst(tour,unsafe_load(g.head))
        tour = vcat(tour[ind:end],tour[1:ind-1])
    end

    map(name,tour)
end


text:
ch_size = 10

@testset "Basics" begin
    worker_inf = PH._launch_workers(ch_size, ch_size)

    @test PH._number_workers(worker_inf) == 1
    @test PH._is_running(worker_inf)

    my_task = @async begin
        PH._send_message(worker_inf, myid(), Ping())
        PH._retrieve_message_type(worker_inf, Ping)
    end
    @test timeout_wait(my_task)

    my_task = @async begin
        PH._abort(worker_inf)
        PH._wait_for_shutdown(worker_inf)
    end
    @test timeout_wait(my_task)

    @test !PH._is_running(worker_inf)
    @test PH._number_workers(worker_inf) == 0
end

@testset "Error Handling" begin
    worker_inf = PH._launch_workers(ch_size, ch_size)

    struct MakeWorkerError <: PH.Message end

    PH._send_message(worker_inf, myid(), MakeWorkerError())
    my_task = @async begin
        PH._wait_for_shutdown(worker_inf)
    end
    if(timeout_wait(my_task))
        @test istaskfailed(my_task)
        @test typeof(my_task.exception) <: RemoteException
    else
        error("Test timed out")
    end
end


text:
module ClusterTrees

include("trees.jl")
include("simpletrees.jl")
include("pointerbasedtrees.jl")

include("octree.jl")
include("leveltree.jl")
include("blocktree.jl")

end # module


text:
function mark_field(r)
    side=Ost
    mark_row(r,side)

    while isborder(r,side)==false
        move!(r,Nord)
        side=inverse(side)
        mark_row(r,side)
    end
 end
function mark_row(r::Robot,side::HorizonSide)
    putmarker!(r)

    while move_if_posible!(r,side) == true
        putmarker!(r)
    end
end

max_stars_repo_path: zadacha 16.jl

max_stars_repo_name: Anpilogova/-

max_issues_repo_path: zadacha 16.jl

max_issues_repo_name: Anpilogova/-

max_forks_repo_path: zadacha 16.jl

max_forks_repo_name: Anpilogova/-


text:
# Tests that the weighting scheme works for binary_classification

module TestWeights

using Test
using DelimitedFiles
using LightGBM

@testset "weightsTest.jl" begin
    try

        LGBM_PATH = ENV["LIGHTGBM_EXAMPLES_PATH"]
        LGBM_PATH = if isabspath(LGBM_PATH) LGBM_PATH else abspath(joinpath(pwd(), "..", LGBM_PATH)) end

        binary_test = readdlm(joinpath(LGBM_PATH, "examples", "binary_classification", "binary.test"), '\t')
        binary_train = readdlm(joinpath(LGBM_PATH, "examples", "binary_classification", "binary.train"), '\t')
        X_train = binary_train[:, 2:end]
        y_train = binary_train[:, 1]
        X_test = binary_test[:, 2:end]
        y_test = binary_test[:, 1]

        binary_test_weight = readdlm(joinpath(LGBM_PATH, "examples", "binary_classification", "binary.test.weight"), '\t')[:,1]
        binary_train_weight = readdlm(joinpath(LGBM_PATH, "examples", "binary_classification", "binary.train.weight"), '\t')[:,1]

        # Test binary estimator.
        estimator = LightGBM.LGBMClassification(
            objective = "binary",
            num_class = 1,
            num_iterations = 20,
            learning_rate = .1,
            early_stopping_round = 1,
            feature_fraction = .8,
            bagging_fraction = .9,
            bagging_freq = 1,
            num_leaves = 1000,
            metric = ["auc", "binary_logloss"],
            is_training_metric = true,
            max_bin = 255,
            min_sum_hessian_in_leaf = 0.,
            min_data_in_leaf = 1,
        )

        # Test fitting.
        LightGBM.fit!(estimator, X_train, y_train, verbosity = -1, weights=binary_train_weight[:,1])
        LightGBM.fit!(estimator, X_train, y_train, (X_test, y_test), verbosity = -1, weights=binary_train_weight[:,1])

        @test true
    catch
        @test false
    end
end

end # module


text:
"""
Harsanyi dividends for S ⊂ N
"""
function Δₕ(G::Game, S::Players)
    G.v(S) - ∑(Δₕ(G, T) for T in propersubsets(S))
end

function degreeHarsanyi(G::GraphGame)
    game = graphtoMyerson(G)
    ϕ = zeros(length(G.N))

    for S in subsets(G.N)
        sg, vmap = induced_subgraph(G.g, S)

        Δ = Δₕ(game, S)
        ds = degree(sg)
        dsum = sum(ds)

        if dsum == 0 continue end

        for (i, d) in enumerate(ds)
            ϕ[vmap[i]] += (Δ * d) / dsum
        end
    end

    return ϕ
end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 80cd4b88-884e-41c5-80c1-30c15fe7e7f4
md"""
## Tutorials

content for this page is coming soon ...
"""

# ╔═╡ d2005d3b-8d74-4308-898e-1a1013a946fe
md"""
* [Welcome to Genie](/docs/tutorials/Overview.html)
* [Installing Genie](/docs/tutorials/Installing-Genie.html)
* [Getting Started](/docs/tutorials/Getting-Started.html)
* [Creating a Web service](/docs/tutorials/Developing-Web-Services.html)
* [Developing MVC web applications](/docs/tutorials/Developing-MVC-Web-Apps.html)
* [Handling URI query params](/docs/tutorials/Handling-Query-Params.html)
* [Working with forms and POST payloads](/docs/tutorials/Working-with-POST-Payloads.html)
* [Using JSON payloads](/docs/tutorials/Using-JSON-Payloads.html)
* [Uploading files](/docs/tutorials/Handling-File-Uploads.html)
* [Adding your libraries into Genie](/docs/tutorials/Publishing-Your-Julia-Code-Online-With-Genie-Apps.html)
* [Loading and starting Genie apps](/docs/tutorials/Loading-Genie-Apps.html)
* [Managing Genie app's dependencies](/docs/tutorials/Managing-External-Packages.html)
* [Advanced routing](/docs/tutorials/Advanced-Routing-Techniques.html)
* [Auto-loading configuration code with initializers](/docs/tutorials/Initializers.html)
* [The secrets file](/docs/tutorials/The-Secrets-File.html)
* [Auto-loading user libraries](/docs/tutorials/The-Lib-Folder.html)
* [Using Genie with Docker](/docs/tutorials/Using-Genie-With-Docker.html)
* [Working with Web Sockets](/docs/tutorials/Working-with-Web-Sockets.html)
* [Force compiling route handlers](/docs/tutorials/Force-Compiling-Routes.html)
* [Deploying to Heroku with Buildpacks](/docs/tutorials/Deploying-With-Heroku-Buildpacks.html)
* [Deploying to Heroku with Nginx](/docs/tutorials/Deploying-Genie-Server-Apps-with-Nginx.html)
"""

# ╔═╡ Cell order:
# ╠═80cd4b88-884e-41c5-80c1-30c15fe7e7f4
# ╟─d2005d3b-8d74-4308-898e-1a1013a946fe


text:
## P O S I T I O N S   A N D   R E F E R E N C E S
"""
Wrapper for cell and range references
"""
abstract type JLXref end
"""
A cell or range reference within the same sheet
"""
abstract type LocalRef <: JLXref end # within a single unnamed sheet
"""
A cell or range reference that specifies the sheet as well as the local position
"""
abstract type GlobalRef <: JLXref end # with reference to the sheet

struct LocalPosition <: LocalRef
    row::UInt
    column::UInt
end
struct CellRange <: LocalRef
    start::LocalPosition
    stop::LocalPosition
end
struct ColumnRange <: LocalRef
    start::UInt
    stop::UInt
end
struct RowRange <: LocalRef
    start::UInt
    stop::UInt
end
const LocalRangeType = Union{CellRange,ColumnRange,RowRange}

struct GlobalRange <: GlobalRef
    sheet::String
    lref::LocalRangeType
end
struct GlobalPosition <: GlobalRef
    sheet::String
    lref::LocalPosition
end
struct SheetRef <: GlobalRef
    sheet::String
end
const Position = Union{LocalPosition,GlobalPosition}

## C E L L S
"""
    A struct that wraps a position together with the contents of the cell (value, formula, etc.)
"""
abstract type Cell end
mutable struct ConstCell <: Cell
    position::GlobalPosition
    value::Any
end
mutable struct DynamicCell <: Cell
    position::GlobalPosition
    formula::String
    func::Function
    ancestors::Array{GlobalRef}
    value::Any
end
# mutable struct CloneCell <: Cell
#     position::GlobalPosition
#     cloneOf::DynamicCell
#     ancestors::Array{GlobalPosition}
#     value::Any
#     #style::JLXstyle
# end

"""
A collection of cells forming a 2D array, belonging to a workbook.
"""
abstract type AbstractSheet end
"""
A collection of sheets, and a map of cause-and-effect
"""
abstract type AbstractWorkbook end

## S H E E T S
mutable struct Sheet <: AbstractSheet
    name::String
    wb::AbstractWorkbook # its "owner"
    map::OrderedDict{LocalPosition, Cell} # mapping positions to cells
end

## W O R K B O O K
mutable struct Workbook <: AbstractWorkbook
    sheets::Dict{String,AbstractSheet}
    map::OrderedDict{GlobalPosition, Cell} # mapping positions to cells
    graph::SimpleDiGraph # directed edge A -> B means A's value is an input to B's formula
    order::Vector{Int} # topological ordering, telling us how to cascade computation (using map)
end

## F I L E S = JWL = Julia Workbook Library
"""
A wrapper combining a workbook with a name, folder location, `Project.toml` and `include.jl` file.
"""
mutable struct JWL
    name::AbstractString
    folder::AbstractString
    wb::Workbook
end


text:
"""
InternalError()

Internal error
"""
struct InternalError <: Exception end

"""
TimeoutError()

Parsing timeout error
"""
struct TimeoutError <: Exception end


text:
"""
    plotMandel( xrange::Tuple{Float64,Float64}, yrange::Tuple{Float64,Float64}; n=100 )

Plots the results of `mandelbrot()` function using the `heatmap()` function.
- Must load `Plots` module.

**Special case**
```
plotMandelbrot( vrange::Tuple{Float64,Float64}; n=100 ) = plotMandelbrot( vrange, vrange; n )
```
"""
function plotMandelbrot( xrange::Tuple{Float64,Float64}, yrange::Tuple{Float64,Float64}; n=100 )
	M = rand(n,n);
	xval = range(xrange[1],xrange[2]; length=n );
	yval = range(yrange[1],yrange[2]; length=n );
	for i in 1:n, j in 1:n
		M[i,j] = mandelbrot( xval[i] + im*yval[j]  )
	end
	heatmap( M, c=:thermal
		,aspect_ratio=1.0
		,xaxis=nothing
		,yaxis=nothing
	)
end
#
plotMandelbrot( vrange::Tuple{Float64,Float64}; n=100 ) = plotMandelbrot( vrange, vrange; n )


text:
using Viruses, Test

ddf = 5.9
@test meltrate(-1, ddf)==0
@test meltrate(1, ddf)==ddf


text:
using Genie.Configuration, Logging

const config = Settings(
  server_port                     = 8000,
  server_host                     = "0.0.0.0",
  log_level                       = Logging.Error,
  log_to_file                     = false,
  server_handle_static_files      = true, # for best performance set up Nginx or Apache web proxies and set this to false
  path_build                      = "build",
  format_julia_builds             = false,
  format_html_output              = false
)

if config.server_handle_static_files
  @warn("For performance reasons Genie should not serve static files (.css, .js, .jpg, .png, etc) in production.
         It is recommended to set up Apache or Nginx as a reverse proxy and cache to serve static assets.")
end

ENV["JULIA_REVISE"] = "off"


text:
# The Val{i} optimizes it into a branching statement (https://discourse.julialang.org/t/meta-programming-an-if-else-statement-of-user-defined-length/53525)
function BINOP!(x::AbstractVector{T}, y::AbstractVector{T}, ::Val{i}, ::Val{clen}, options::Options) where {i,clen,T<:Real}
    op = options.binops[i]
    # broadcast!(op, x, x, y)
    @inbounds @simd for j=1:clen
        x[j] = op(x[j], y[j])
    end
end

function UNAOP!(x::AbstractVector{T}, ::Val{i}, ::Val{clen}, options::Options) where {i,clen,T<:Real}
    op = options.unaops[i]
    @inbounds @simd for j=1:clen
        x[j] = op(x[j])
    end
end

#isnan(x::AbstractFloat) = (x != x)::Bool
#isfinite(x::AbstractFloat) = x - x == 0

# Evaluate an equation over an array of datapoints
function evalTreeArray(tree::Node, cX::AbstractMatrix{T}, options::Options)::Tuple{AbstractVector{T}, Bool} where {T<:Real}
    clen = size(cX)[1]
    if tree.degree == 0
        if tree.constant #TODO: Make this done with types instead
            return (fill(convert(T, tree.val), clen), true)
        else
            return (copy(cX[:, tree.feature]), true)
        end
    end

    (cumulator, complete) = evalTreeArray(tree.l, cX, options)
    if !complete
        return (cumulator, false)
    end

    if tree.degree == 1
        op_idx = tree.op
        UNAOP!(cumulator, Val(op_idx), Val(clen), options)
    else
        (array2, complete2) = evalTreeArray(tree.r, cX, options)
        if !complete2
            return (cumulator, false)
        end
        op_idx = tree.op
        BINOP!(cumulator, array2, Val(op_idx), Val(clen), options)
    end

	# TODO - consider checking this at end.
    @inbounds @simd for i=1:clen
        if (cumulator[i] != cumulator[i])::Bool || (cumulator[i] - cumulator[i] != convert(T, 0))::Bool
            return (cumulator, false)
        end
    end

    return (cumulator, true)
end


text:
include("libavcodec_h.jl")

include("avcodec.jl")
#include("dv_profile.jl")
#include("vdpau.jl")


text:
# InverseGamma distribution

export InverseGamma

@parameterized InverseGamma(shape) ≃ Lebesgue(ℝ₊)

function logdensity(μ::InverseGamma{(:shape,)}, x)
    α = μ.shape
    xinv = 1/x

    return xlogy(α + 1, xinv) - xinv - loggamma(α)
end

Base.rand(rng::AbstractRNG, T::Type, μ::InverseGamma{(:shape,)}) = rand(rng, Dists.InverseGamma(μ.shape))

≪(::InverseGamma, ::Lebesgue{X}) where X <: Real = true

TV.as(::InverseGamma) = asℝ₊

# @μσ_methods InverseGamma(shape)


text:
test_msh = (
    msh_version = "4.1",
    raw_PhysicalNames = SubString{String}[],
    raw_Entities = SubString{String}[
        "4 4 1 0",
        "1 0 0 0 0 ",
        "2 1 0 0 0 ",
        "3 1 1 0 0 ",
        "4 0 1 0 0 ",
        "1 0 0 0 1 0 0 0 2 1 -2 ",
        "2 1 0 0 1 1 0 0 2 2 -3 ",
        "3 0 1 0 1 1 0 0 2 3 -4 ",
        "4 0 0 0 0 1 0 0 2 4 -1 ",
        "1 0 0 0 1 1 0 0 4 1 2 3 4 "
    ],
    raw_Nodes = SubString{String}[
        "9 30 1 30",
        "0 1 0 1",
        "1",
        "0 0 0",
        "0 2 0 1",
        "2",
        "1 0 0",
        "0 3 0 1",
        "3",
        "1 1 0",
        "0 4 0 1",
        "4",
        "0 1 0",
        "1 1 0 3",
        "5",
        "6",
        "7",
        "0.2499999999994121 0 0",
        "0.499999999998694 0 0",
        "0.7499999999993416 0 0",
        "1 2 0 3",
        "8",
        "9",
        "10",
        "1 0.2499999999994121 0",
        "1 0.499999999998694 0",
        "1 0.7499999999993416 0",
        "1 3 0 3",
        "11",
        "12",
        "13",
        "0.7500000000003465 1 0",
        "0.5000000000020591 1 0",
        "0.2500000000010404 1 0",
        "1 4 0 3",
        "14",
        "15",
        "16",
        "0 0.7500000000003465 0",
        "0 0.5000000000020591 0",
        "0 0.2500000000010404 0",
        "2 1 0 14",
        "17",
        "18",
        "19",
        "20",
        "21",
        "22",
        "23",
        "24",
        "25",
        "26",
        "27",
        "28",
        "29",
        "30",
        "0.3640932128839348 0.7867687832230399 0",
        "0.2110423218192967 0.3753337416957934 0",
        "0.6146003733357942 0.2544869979106268 0",
        "0.8080509614599705 0.6243689539557383 0",
        "0.1979210047255827 0.6215671767908416 0",
        "0.4308090314147045 0.5056502726999197 0",
        "0.3712961061742487 0.2167990069479849 0",
        "0.6139924015723897 0.7478815024664035 0",
        "0.8112741566735278 0.3855439771726542 0",
        "0.6557453848912773 0.5035863408410685 0",
        "0.1830127018930926 0.816987298107279 0",
        "0.816987298107459 0.8169872981071898 0",
        "0.1830127018926974 0.1830127018931337 0",
        "0.8351749060017326 0.1780061950165386 0"
    ],
    raw_Elements = SubString{String}[
        "9 62 1 62",
        "0 1 15 1",
        "1 1 ",
        "0 2 15 1",
        "2 2 ",
        "0 3 15 1",
        "3 3 ",
        "0 4 15 1",
        "4 4 ",
        "1 1 1 4",
        "5 1 5 ",
        "6 5 6 ",
        "7 6 7 ",
        "8 7 2 ",
        "1 2 1 4",
        "9 2 8 ",
        "10 8 9 ",
        "11 9 10 ",
        "12 10 3 ",
        "1 3 1 4",
        "13 3 11 ",
        "14 11 12 ",
        "15 12 13 ",
        "16 13 4 ",
        "1 4 1 4",
        "17 4 14 ",
        "18 14 15 ",
        "19 15 16 ",
        "20 16 1 ",
        "2 1 2 42",
        "21 19 22 23 ",
        "22 17 22 24 ",
        "23 22 19 26 ",
        "24 24 22 26 ",
        "25 22 18 23 ",
        "26 17 21 22 ",
        "27 20 24 26 ",
        "28 12 13 17 ",
        "29 15 16 18 ",
        "30 9 10 20 ",
        "31 6 7 19 ",
        "32 19 25 26 ",
        "33 14 15 21 ",
        "34 11 12 24 ",
        "35 5 6 23 ",
        "36 8 9 25 ",
        "37 21 18 22 ",
        "38 15 18 21 ",
        "39 12 17 24 ",
        "40 6 19 23 ",
        "41 9 20 25 ",
        "42 4 14 27 ",
        "43 3 11 28 ",
        "44 13 4 27 ",
        "45 16 1 29 ",
        "46 1 5 29 ",
        "47 10 3 28 ",
        "48 7 2 30 ",
        "49 2 8 30 ",
        "50 14 21 27 ",
        "51 11 24 28 ",
        "52 5 23 29 ",
        "53 8 25 30 ",
        "54 17 13 27 ",
        "55 20 10 28 ",
        "56 18 16 29 ",
        "57 19 7 30 ",
        "58 21 17 27 ",
        "59 23 18 29 ",
        "60 24 20 28 ",
        "61 25 19 30 ",
        "62 25 20 26 "
    ]
)