text:
using Test
using LinearAlgebra
using MathProgBase.SolverInterface

function linprogsolvertest(solver::AbstractMathProgSolver, eps = Base.rtoldefault(Float64))
    @testset "Testing LP interface with $solver" begin
        @testset "Basic interface" begin
            m = LinearQuadraticModel(solver)

            # Min -x
            # s.t. x + y <= 1
            # x, y >= 0
            loadproblem!(m, [1. 1.], [0.,0.], [Inf, Inf], [-1.,0.], [-Inf],[1.], :Min)
            @test numvar(m) == 2
            @test numconstr(m) == 1
            @test getvarLB(m) == [0.,0.]
            @test getvarUB(m)[1] >= 1e20
            @test getvarUB(m)[2] >= 1e20
            @test getconstrLB(m)[1] < -1e30
            @test getconstrUB(m)[1] == 1.0
            @test getobj(m) == [-1.,0.]
            @test getsense(m) == :Min

            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), -1, atol=eps)
            @test isapprox(norm(getsolution(m) - [1.0, 0.0]), 0.0, atol=eps)
            @test isapprox(getconstrsolution(m)[1], 1.0, atol=eps)
            @test isapprox(getconstrduals(m)[1], -1.0, atol=eps)
            @test isapprox(norm(getreducedcosts(m) - [0.0, 1.0]), 0.0, atol=eps)

            setsense!(m, :Max)
            # max x
            setobj!(m, [1.0,0.0])
            optimize!(m)

            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 1, atol=eps)
            @test isapprox(norm(getsolution(m) - [1.0, 0.0]), 0.0, atol=eps)
            @test isapprox(getconstrsolution(m)[1], 1.0, atol=eps)
            @test isapprox(getconstrduals(m)[1], 1.0, atol=eps)
            @test isapprox(norm(getreducedcosts(m) - [0.0, -1.0]), 0.0, atol=eps)

            # add new variable to get:
            # max x + 2z
            # s.t. x + y + z <= 1
            # x,y,z >= 0
            addvar!(m, [1], [1.0], 0, Inf, 2.0)

            @test numvar(m) == 3
            @test numconstr(m) == 1

            optimize!(m)

            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 2, atol=eps)
            @test isapprox(norm(getsolution(m) - [0.0, 0.0, 1.0]), 0.0, atol=eps)
            @test isapprox(getconstrsolution(m)[1], 1.0, atol=eps)
            @test isapprox(getconstrduals(m)[1], 2.0, atol=eps)
            @test isapprox(norm(getreducedcosts(m) - [-1.0, -2.0, 0.0]), 0.0, atol=eps)

            setvarLB!(m, [-1.0,0.0,0.0])
            optimize!(m)

            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 3, atol=eps)

            # fix z to zero
            setvarLB!(m, [0.0,0.0,0.0])
            setvarUB!(m, [Inf,Inf,0.0])
            optimize!(m)

            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 1, atol=eps)

            setconstrUB!(m, [2.0])
            setconstrLB!(m, [2.0])
            optimize!(m)
            @test isapprox(getobjval(m), 2, atol=eps)

            setobj!(m, [1.0, 2.0, 0.0])
            optimize!(m)
            @test isapprox(getobjval(m), 4, atol=eps)
            @test isapprox(norm(getsolution(m) - [0.0, 2.0, 0.0]), 0.0, atol=eps)

            # we now have:
            # max x+2y
            # s.t. x + y + z == 2
            # x,y >= 0, z = 0

            # add constraint x - y >= 0
            addconstr!(m, [1,2], [1.0,-1.0], 0.0, Inf)
            optimize!(m)

            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 3, atol=eps)
            @test isapprox(norm(getsolution(m) - [1.0, 1.0, 0.0]), 0.0, atol=eps)
            @test isapprox(norm(getconstrduals(m) - [1.5, -0.5]), 0.0, atol=eps)
            @test isapprox(norm(getreducedcosts(m) - [0.0, 0.0, -1.5]), 0.0, atol=eps)
        end

        @testset "addvar! interface" begin
            m = LinearQuadraticModel(solver)

            # Min -x
            # s.t. x + y <= 1
            # x, y >= 0
            addvar!(m, 0, Inf, -1)
            addvar!(m, 0, Inf, 0)
            addconstr!(m, [1, 2], [1.0, 1.0], -Inf, 1.0)
            @test numvar(m) == 2
            @test numconstr(m) == 1
            @test getvarLB(m) == [0.,0.]
            @test getvarUB(m)[1] >= 1e20
            @test getvarUB(m)[2] >= 1e20
            @test getconstrLB(m)[1] <= -1e20
            @test getconstrUB(m)[1] == 1.0
            @test getobj(m) == [-1.,0.]
            @test getsense(m) == :Min

            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), -1, atol=eps)
            @test isapprox(norm(getsolution(m) - [1.0, 0.0]), 0.0, atol=eps)
            @test isapprox(getconstrsolution(m)[1], 1.0, atol=eps)
            @test isapprox(getconstrduals(m)[1], -1.0, atol=eps)
            @test isapprox(norm(getreducedcosts(m) - [0.0, 1.0]), 0.0, atol=eps)
        end

        @testset "setconstrLB! and setconstrUB!" begin
            ####################################
            # test setconstrLB! and setconstrUB!
            # Test that:
            #   Modifying lower bound works
            #   Modifying upper bound works
            #   Setting upper and lower bound to same value works
            #

            m = LinearQuadraticModel(solver)
            # Min  x - y
            # s.t. 0.0 <= x
            #             y <= 0.0
            # x,y unbounded

            loadproblem!(m, [ 1.0 0.0 ; 0.0 1.0 ], [-Inf, -Inf], [Inf,Inf], [1.0, -1.0], [0.0,-Inf], [Inf,0.0], :Min)

            optimize!(m)

            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0, atol=eps)
            @test isapprox(norm(getsolution(m)), 0.0, atol=eps)

            # Min  x - y
            # s.t. 100.0 <= x
            #               y <= 0.0
            # x,y unbounded
            setconstrLB!(m,[100.0,-Inf])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 100.0, atol=eps)
            @test isapprox(norm(getsolution(m) - [ 100.0, 0.0 ]), 0.0, atol=eps)

            # Min  x - y
            # s.t. 100.0 <= x
            #               y <= -100.0
            # x,y unbounded
            setconstrUB!(m,[Inf,-100.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 200.0, atol=eps)
            @test isapprox(norm(getsolution(m) - [ 100.0, -100.0 ]), 0.0, atol=eps)
        end

        @testset "Issue #40 from Gurobi.jl" begin
            # Test issue #40 from Gurobi.jl
            # min  x
            # s.t. x >= 0
            #      x >= 3
            m = LinearQuadraticModel(solver)
            loadproblem!(m, [1.0 1.0]', [-Inf], [Inf], [1.0], [0.0, 3.0], [Inf, Inf], :Min)
            optimize!(m)
            for i = 1:length(getconstrLB(m))
                @test getconstrLB(m)[i] <= getconstrsolution(m)[i] + eps
                @test getconstrsolution(m)[i] <= getconstrUB(m)[i] + eps
            end

            # min  x
            # s.t. x <= 0
            #      x <= 3
            m = LinearQuadraticModel(solver)
            loadproblem!(m, [1.0 1.0]', [-Inf], [Inf], [1.0], [-Inf, -Inf], [0.0, 3.0], :Max)
            optimize!(m)
            for i = 1:length(getconstrLB(m))
                @test getconstrLB(m)[i] <= getconstrsolution(m)[i] + eps
                @test getconstrsolution(m)[i] <= getconstrUB(m)[i] + eps
            end
        end

        @testset "Change coeffs, del constr, del var" begin
            #####################################
            # Start from simple LP
            # Solve it
            # Copy and solve again
            # Chg coeff, solve, change back solve
            # del constr and solve
            # del var and solve

            #   maximize x + y
            #
            #   s.t. 2 x + 1 y <= 4
            #        1 x + 2 y <= 4
            #        x >= 0, y >= 0
            #
            #   solution: x = 1.3333333, y = 1.3333333, objv = 2.66666666

            m = LinearQuadraticModel(solver)

            loadproblem!(m, [ 2.0 1.0 ; 1.0 2.0 ], [0.0,0.0], [Inf,Inf], [1.0, 1.0], [-Inf, -Inf], [4.0,4.0], :Max)

            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 2.6666666666, atol=eps)
            @test isapprox(norm(getsolution(m) - [1.3333333333, 1.3333333333]), 0.0, atol=eps)

            # copy and solve again

            if applicable(copy, m)
                m2 = copy(m)
                optimize!(m2)
                @test status(m2) == :Optimal
                @test isapprox(getobjval(m2), 2.6666666666, atol=eps)
                @test isapprox(norm(getsolution(m2) - [1.3333333333, 1.3333333333]), 0.0, atol=eps)
            end

            # change coeff
            #   maximize x + y
            #
            #   s.t. 2 x + 2 y <= 4
            #        1 x + 2 y <= 4
            #        x >= 0, y >= 0
            #
            #   solution: x = 0, y = 2, objv = 2
            if applicable(changecoeffs!, m, [1], [2],  [2.])
                changecoeffs!(m, [1], [2],  [2.])
                optimize!(m)
                @test status(m) == :Optimal
                @test isapprox(getobjval(m), 2.0, atol=eps)
                @test isapprox(norm(getsolution(m) - [0.0, 2.0]), 0.0, atol=eps)
            end

            # delconstrs and solve
            #   maximize x + y
            #
            #   s.t. 1 x + 2 y <= 4
            #        x >= 0, y >= 0
            #
            #   solution: x = 4, y = 0, objv = 4
            if applicable(delconstrs!, m, [1])
                delconstrs!(m, [1])
                optimize!(m)
                @test status(m) == :Optimal
                @test isapprox(getobjval(m), 4.0, atol=eps)
                @test isapprox(norm(getsolution(m) - [4.0, 0.0]), 0.0, atol=eps)
            end

            # delvars and solve
            #   maximize y
            #
            #   s.t.  2 y <= 4
            #           y >= 0
            #
            #   solution: y = 2, objv = 2
            if applicable(delvars!, m, [1])
                delvars!(m, [1])
                optimize!(m)
                @test status(m) == :Optimal
                @test isapprox(getobjval(m), 2.0, atol=eps)
                @test isapprox(getsolution(m)[1], 2.0, atol=eps)
            end
        end
    end
end

function linprogsolvertestextra(solver::AbstractMathProgSolver; eps = Base.rtoldefault(Float64))
    @testset "Testing LP interface extra with $solver" begin
        ####################################
        # test setconstrLB! and setconstrUB!
        # Test that:
        #   Modifying lower bound works
        #   Modifying upper bound works
        #   Setting upper and lower bound to same value works
        #
        @testset "setconstrLB! and setcontrUB!" begin
            m = LinearQuadraticModel(solver)
            # Min  x - y
            # s.t. 0.0 <= x <= 0.0
            #      0.0 <= y <= 0.0
            # x,y unbounded
            loadproblem!(m, [ 1.0 0.0 ; 0.0 1.0 ], [-Inf, -Inf], [Inf,Inf], [1.0, -1.0], [0.0,0.0], [0.0,0.0], :Min)

            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0)
            @test isapprox(norm(getsolution(m)), 0.0)

            # Min  x - y
            # s.t. 0.0 <= x <= 100.0
            #      0.0 <= y <= 100.0
            # x,y unbounded
            setconstrUB!(m,[100.0,100.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), -100.0)
            @test isapprox(norm(getsolution(m) - [ 0.0, 100.0 ]), 0.0)

            # Min  x - y
            # s.t. -100.0 <= x <= 100.0
            #      -100.0 <= y <= 100.0
            # x,y unbounded
            setconstrLB!(m,[-100.0,-100.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), -200.0)
            @test isapprox(norm(getsolution(m) - [ -100.0, 100.0 ]), 0.0)

            # Min  x - y
            # s.t. -100.0 <= x <= 100.0
            #      -100.0 <= y <= 100.0
            # x,y unbounded
            setconstrLB!(m,[10.0,10.0])
            setconstrUB!(m,[10.0,10.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0)
            @test isapprox(norm(getsolution(m) - [ 10.0, 10.0 ]), 0.0)

            # Min  x - y
            # s.t. 0.0  <= x <= Inf
            #      -Inf <= y <= 0.0
            # x,y unbounded
            setconstrLB!(m,[0.0,-Inf])
            setconstrUB!(m,[Inf,0.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0)
            @test isapprox(norm(getsolution(m)), 0.0)

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= 0.0
            # x,y unbounded
            setconstrLB!(m,[-Inf,-Inf])
            setconstrUB!(m,[Inf,0.0])
            optimize!(m)
            @test status(m) == :Unbounded
            @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8

            # Min  x - y
            # s.t. 0.0 <= x <= Inf
            #      -Inf <= y <= Inf
            # x,y unbounded
            setconstrLB!(m,[0.0,-Inf])
            setconstrUB!(m,[Inf,Inf])
            optimize!(m)
            @test status(m) == :Unbounded
            @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8
        end

        ####################################
        # test setvarLB! and setvarUB!
        # Test that:
        #   Modifying lower bound works
        #   Modifying upper bound works
        #   Setting upper and lower bound to same value works
        @testset "setvarLB! and setvarUB!" begin
            m = LinearQuadraticModel(solver)
            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            # 0.0 <= x <= 0.0
            # 0.0 <= y <= 0.0
            loadproblem!(m, [ 1.0 0.0 ; 0.0 1.0 ], [0.0,0.0], [0.0,0.0], [1.0, -1.0], [-Inf, -Inf], [Inf,Inf], :Min)

            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0)
            @test isapprox(norm(getsolution(m)), 0.0)

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            # 0.0 <= x <= 100.0
            # 0.0 <= y <= 100.0

            # x,y unbounded
            setvarUB!(m,[100.0,100.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), -100.0)
            @test isapprox(norm(getsolution(m) - [ 0.0, 100.0 ]), 0.0)

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            # -100.0 <= x <= 100.0
            # -100.0 <= y <= 100.0
            setvarLB!(m,[-100.0,-100.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), -200.0)
            @test isapprox(norm(getsolution(m) - [ -100.0, 100.0 ]), 0.0)

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            # -100.0 <= x <= 100.0
            # -100.0 <= y <= 100.0
            setvarLB!(m,[10.0,10.0])
            setvarUB!(m,[10.0,10.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0)
            @test isapprox(norm(getsolution(m) - [ 10.0, 10.0 ]), 0.0)

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            # 0.0  <= x <= Inf
            # -Inf <= y <= 0.0
            # x,y unbounded
            setvarLB!(m,[0.0,-Inf])
            setvarUB!(m,[Inf,0.0])
            optimize!(m)
            @test status(m) == :Optimal
            @test isapprox(getobjval(m), 0.0)
            @test isapprox(norm(getsolution(m)), 0.0)

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            # -Inf <= x <= Inf
            # -Inf <= y <= 0.0
            setvarLB!(m,[-Inf,-Inf])
            setvarUB!(m,[Inf,0.0])
            optimize!(m)
            @test status(m) == :Unbounded
            @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8

            # Min  x - y
            # s.t. -Inf <= x <= Inf
            #      -Inf <= y <= Inf
            #  0.0 <= x <= Inf
            # -Inf <= y <= Inf
            # x,y unbounded
            setvarLB!(m,[0.0,-Inf])
            setvarUB!(m,[Inf,Inf])
            optimize!(m)
            @test status(m) == :Unbounded
            @test dot(getunboundedray(m),[1.0,-1.0]) < 1e-8
        end
    end
end


text:
using MocosSim
using CodecZlib
using FileIO
using JLD2
using DataFrames
using FixedPointNumbers
using Random
using Profile

individuals_df = load( (@__DIR__) * "/wroclaw_v4.jld2", "individuals_df")
rng = Random.MersenneTwister(0);

@info "creating state"
@time state = MocosSim.SimState(nrow(individuals_df))

@info "loading data"
@time params = MocosSim.load_params(
  population=individuals_df,
  mild_detection_prob=0.4,
  backward_tracing_prob=0.2,
  forward_tracing_prob=0.2,
  constant_kernel_param=0.7425,
  household_kernel_param=0.07,
  british_strain_multiplier=1.7,

  age_coupling_thresholds=[0, 18, 40],
  age_coupling_weights=[1.1 1.2 1.3;
                        2.1 2.2 2.3;
                        3.1 3.2 3.3],
  age_coupling_use_genders=false,
  age_coupling_param=1.0,

  infection_modulation_name="TanhModulation",
  infection_modulation_params=(
    scale=12500,
    loc=300000,
    weight_detected=1,
    weight_deaths=0,
    limit_value=0.545455
  ),

  forward_detection_delay=1.75,
  backward_detection_delay=1.75,
  testing_time=3.0
);

mutable struct Callback
  max_time::Float64
  num_events::Int
end

function (cb::Callback)(event::MocosSim.Event, ::MocosSim.SimState, ::MocosSim.SimParams)::Bool
  cb.num_events += 1
  MocosSim.time(event) < cb.max_time
end

@info "initializing state"
@time MocosSim.reset!(state, MersenneTwister(0))
@time MocosSim.InstantOutsideCases(;num_infections=100)(state, params)

@info "warm-up"
cb = Callback(50, 0)
@time MocosSim.simulate!(state, params, cb)
@info "events executed = $(cb.num_events), simtime=$(state.time)"

@info "measuring allocations"
cb = Callback(100, 0)
Profile.clear_malloc_data()
@time MocosSim.simulate!(state, params, cb)
@info "events executed = $(cb.num_events)"

ext: jl


text:
# UMLS
# ====
#
# Julia utilities to interact with the Unified Medical Language System (UMLS) REST API
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioServices.jl/blob/master/LICENSE.md

"""
    BioServices.UMLS

Julia interface to [Unified Medical Language REST API](https://github.com/HHS/uts-rest-api)

The UMLS provides basic functionality to authenticate and query the UMLS API.
The following functions are exported from this module:
1. `get_tgt`: Get Ticket-Granting Ticket
2. `search_umls`: Search the UMLS
3. `best_match_cui`: Get the best matching CUI for a term
4. `get_cui`: Get Concept ID
5. `get_semantic_types`: Get UMLS semantic type

See "UMLS REST API Home Page"
(https://documentation.uts.nlm.nih.gov/rest/home.html) for more details.
"""
module UMLS

export get_tgt,
       search_umls,
       best_match_cui,
       get_cui,
       get_semantic_types

using Gumbo
import HTTP
import JSON
using Dates

#------------- Endpoints -----------------------
const uri = "https://utslogin.nlm.nih.gov"
const service = "http://umlsks.nlm.nih.gov"
const rest_uri = "https://uts-ws.nlm.nih.gov"
#----------------------------------------------

#-------------- Common Utils--------------------
struct BadResponseException <: Exception
    code::Int64
end

function Base.showerror(io::IO, e::BadResponseException)
    print(io, "BadResponseException with code: ",STATUS_CODES[e.code])
end

#-------------- Authentication -------------------
"""
    time_to_last_save(file)

Get how many hours since last save
"""
function time_to_last_save(file)
    #unix time is GMT
    time_diff = Dates.value(now() - Dates.unix2datetime(mtime(file)))/ (1000 * 60 * 60)
    return time_diff
end

"""
    tgt_exists(; tgt_file = "UTS_TGT.txt", hours = 1)

Check if there is a TGT in disk and if it has not expired
"""
function tgt_exists(; tgt_file = "UTS_TGT.txt", hours = 1)
    if isfile(tgt_file)
        #check time
        time_elapsed = time_to_last_save(tgt_file)
        # Expiration time should be 8 hours - but I tend to expirience bad TGT after few hours
        if time_elapsed > hours
            @info "UTS TGT Expired"
            rm(tgt_file)
            return false
        else
            return true
        end
    end
    return false
end

"""
    get_tgt(; force_new::Bool = false, kwargs...)

Retrieve a ticket granting ticket (TGT) using

1. UTS username and password OR
2. apikey

A tgt is valid for 8 hours. Therefore, look for UTS_TGT.txt in the local
directory to see if it has been recently stored. One can force getting a
new ticket by passing keyword argument `force_new=true`

####Examples

```julia
tgt = get_tgt(username = "myuser", password = "mypass")
```

```julia
tgt = get_tgt(apikey = "mykey")
```
"""
function get_tgt(; force_new::Bool = false, kwargs...)
    params = Dict(kwargs)
    body = HTTP.escapeuri(params)
    auth_endpoint = "/cas/v1/tickets/"

    if !haskey(params, :apikey)
        if !haskey(params, :username) && !haskey(params, :password)
            error("UMLSRestAPI: Connection requieres passing keword arguments: apikey or username and password")
        end
    else
        auth_endpoint = "/cas/v1/api-key"
    end

    # Check if there is a valid ticket on disk
    tgt_file = "UTS_TGT.txt"
    if tgt_exists(tgt_file=tgt_file) && !force_new
        @info "UTS: Reading TGT from file"
        return readline(tgt_file)
    end

    @info "UTS: Requesting new TGT"
    headers = Dict("Content-type"=> "application/x-www-form-urlencoded",
    "Accept"=> "text/plain", "User-Agent"=>"julia" )
    r = HTTP.request("POST", uri*auth_endpoint, body=body, headers=headers)
    ascii_r = String(r.body)

    doc = parsehtml(ascii_r)
    #for now - harcoded
    #TO DO:: parse and check
    ticket = ""
    try
        ticket = getattr(doc.root.children[2].children[2], "action")
    catch
        error("Could not get TGT: Unexpected structure of UTS response")
    end

    open(tgt_file, "w") do f
        write(f, ticket)
    end

    return ticket
end

"""
    get_ticket(tgt)

Retrieve a single-use Service Ticket using TGT
"""
function get_ticket(tgt)
    params = Dict("service"=> service)
    body = HTTP.escapeuri(params)

    headers = Dict("Content-type"=> "application/x-www-form-urlencoded",
    "Accept"=> "text/plain", "User-Agent"=>"JuliaBioServices" )
    r = HTTP.Response(503)
    try
        r = HTTP.request("POST", tgt; body=body, headers=headers, retry_non_idempotent=true)
    catch
        isdefined(r, :code) ? error("UMLS GET error: ", r.code) : error("UMLS COULD NOT GET")
    end
    return String(r.body)
end

#------------------Search ---------------------
"""
    search_umls(tgt, query)

Search UMLS Rest API. For more info see
[UMLS_API](https://documentation.uts.nlm.nih.gov/rest/search/)

####Arguments

- `tgt`: Ticket Granting Ticket
- `query`: UMLS query containing the search term
- `version:` Optional - defaults to current

####Output

- `result_pages`: Array, where each entry is a dictionary containing a page of
results. e.g
```julia
Dict{AbstractString,Any} with 3 entries:
"pageSize"   => 25
"pageNumber" => 1
"result"     => Dict{AbstractString,Any}("classType"=>"searchResults","result…
```

####Examples

```julia
credentials = Credentials(user, psswd)
tgt = get_tgt(credentials)
term = "obesity"
query = Dict("string"=>term, "searchType"=>"exact" )
all_results= search_umls(tgt, query)
```
"""
function search_umls(tgt, query; version::String="current", timeout=1)
    page = 0
    content_endpoint = "/rest/search/current"

    #each page of results is appended to the output list
    #where each entry is a dictionary containing that pages's results
    # e.g
    # Dict{AbstractString,Any} with 3 entries:
    #   "pageSize"   => 25
    #   "pageNumber" => 1
    #   "result"     => Dict{AbstractString,Any}("classType"=>"searchResults","result…
    result_pages = Array{Any,1}()

    while true

        #get a new ticket per page if necessary
        ticket = ""
        ticket = get_ticket(tgt)

        page +=1
        #append ticket to query
        query["ticket"]= ticket
        query["pageNumber"]= string(page)

        r = HTTP.request("GET", rest_uri*content_endpoint, query=query, timeout=timeout)

        if r.status != 200
            error("Bad HTTP status $(r.status)")
        end

        json_response = JSON.parse(String(r.body))
        # println("No Results ", length(json_response["result"]["results"]))
        if json_response["result"]["results"][1]["ui"]=="NONE"
            break
        end
        push!(result_pages,json_response)
    end

    return result_pages
end

"""
    best_match_cui(result_pages)

Retrieve the best match from array of all result pages

####Example

```julia
cui = best_match_cui(all_results)
```
"""
function best_match_cui(result_pages)
    return result_pages[1]["result"]["results"][1]["ui"]
end

#------------------Content--------------------------------
"""
    get_cui(tgt,cui)

Retrieve information (name, semantic types, number of atoms, etc) for a known CUI
from latest UMLS version or a specific release.

Returns UTS json response

See: https://documentation.uts.nlm.nih.gov/rest/concept

####Example

```julia
tgt = get_tgt(apikey = "mykey")
cui = "C0028754"
concept = get_cui(tgt, cui)
```
"""
function get_cui(tgt,cui; version="current")
    content_endpoint = "/rest/content/$(version)/CUI/"*cui
    #get a new ticket
    ticket = ""
    try
        ticket = get_ticket(tgt)
    catch err
        rethrow(err)
    end

    r = HTTP.Response(503)
    try
        r = HTTP.request("GET", rest_uri*content_endpoint,query=Dict("ticket"=> ticket))
    catch
        isdefined(r, :code) ? error("UMLS GET error: ", r.code) : error("UMLS COULD NOT GET")
    end

    return JSON.parse(String(r.body))
end

"""
    get_semantic_types(c::Credentials, cui)

Return an array of the semantic types associated with a cui

####Example

```julia
tgt = get_tgt(apikey = "mykey")
cui = "C0028754"
sm = get_semantic_types(tgt, cui)
```
"""
function get_semantic_types(tgt, cui; version="current")
    json_response = get_cui(tgt,cui; version=version)
    st = json_response["result"]["semanticTypes"]
    return [String(concept["name"]) for concept in st]
end

end

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
function init_xmins(data::AbstractArray,xmins::AbstractArray,xmax::Int64)
  sorted_data = sort(data)
  bins_data = array_bins(sorted_data)
  if (xmins == [])
    xmins = sort(unique(sorted_data))
    if (xmins[end] > xmax)
      max_indx = findlast(x -> x <= xmax, xmins)
      xmins = xmins[1:max_indx]
    end
    xmins = xmins[1:end]
    if xmins[1] <= 0
      println("removing elements smaller than 0")
      xmins = xmins[findfirst(x -> x > 0,xmins), end]
    end
  else
    real_xmins = Array(Bool,length(xmins))
    for i=1:length(xmins)
      if (haskey(bins_data,xmins[i]))
        real_xmins[i] = true
      else
        real_xmins[i] = false
      end
    end
    xmins = xmins[real_xmins]
  end
  return sorted_data,bins_data,xmins
end

function _estimate_xmin(sorted_data::AbstractArray,bins_data::Dict{Float64,Int64},distribution::DataType,xmins::AbstractArray,xmax::Int64)
  if (length(xmins) == 0)
    println("No xmins")
    return Union{}
  end

  min_dist = Inf
  best_fit = Union{}

  for xmin in xmins
    fit_data = sorted_data[bins_data[xmin]:end]
    f = fit(distribution,fit_data)

    negloglike(alpha) = begin
      d = distribution(alpha[1],f.θ)
      r = -sum(logpdf(d,fit_data))
      if(Inf == r || -Inf == r)
        r = 1e12
      end
      return r
    end
    try
      opt_alfa = fminbox(DifferentiableFunction(negloglike),[f.α],[1.0],[Inf])
      f = distribution(opt_alfa.minimum[1],f.θ)
    catch
      #if fminbox throws error it means that function cannot be optimized
    end
    if (f.α == Inf)
      continue
    end
    d = Kolmogorov_smirnov_test(fit_data,f,xmin,xmax)

    if ((min_dist > d))
      best_fit = f
      min_dist = d
    end
  end
  return best_fit,min_dist
end

function estimate_xmin(data::AbstractArray,distribution::Type{con_powerlaw};xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5))

  sorted_data,bins_data,xmins = init_xmins(data,xmins,xmax)

  _estimate_xmin(sorted_data,bins_data,distribution,xmins,xmax)
end

function estimate_xmin(data::AbstractArray,distribution::Type{dis_powerlaw};xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5))
  if (data != floor(data))
    println("Data should be discreate. Use round or floor function.")
    return Union{}
  end

  sorted_data,bins_data,xmins = init_xmins(data,xmins,xmax)
  if xmins[1] < 1
    xmins = xmins[findfirst(x -> x >= 1), end]
  end

  _estimate_xmin(sorted_data,bins_data,distribution,xmins,xmax)
end

function estimate_xmin(data::AbstractArray,distribution::DataType;xmins::AbstractArray = [],xmax::Int64 = round(Int,1e5))
  min_dist = Inf
  best_fit = Union{}
  if ((distribution <: DiscreteUnivariateDistribution)&&(data != floor(data)))
    println("Data should be discreate. Use round or floor function.")
    return Union{}
  end

  sorted_data,bins_data,xmins = init_xmins(data,xmins,xmax)
  if ((distribution <: DiscreteUnivariateDistribution)&&(xmins[1] < 1))
    xmins = xmins[findfirst(x -> x >= 1), end]
  end

  if (length(xmins) == 0)
    println("No xmins")
    return Union{}
  end

  for xmin in xmins
    fit_data = sorted_data[bins_data[xmin]:end]
    f = fit(distribution,fit_data)
    d = Kolmogorov_smirnov_test(fit_data,f,xmin,xmax)

    if ((min_dist > d))
      best_fit = f
      min_dist = d
    end
  end
  return best_fit,min_dist
end

#helper function  return dict{Int,Int} where key is number and value is where in array it starts
#input array must be sorted
function array_bins(arr::AbstractArray)
  num = arr[1]
  bins = Dict{Float64,Int64}()
  bins[num] = 1
  for i=1: length(arr)
    if num != arr[i]
      num = arr[i]
      bins[num] = i
    end
  end
  return bins
end


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

module TestContext

using MXNet
using Test

function test_num_gpus()
  @info "Context::num_gpus"

  @test num_gpus() >= 0
end

function test_context_macro()
  @info "Context::@context"

  @context mx.CPU 42 begin
    ctx = mx.current_context()
    @test ctx.device_type == mx.CPU
    @test ctx.device_id   == 42

    @context mx.GPU 24 begin
      ctx = mx.current_context()
      @test ctx.device_type == mx.GPU
      @test ctx.device_id   == 24
    end

    ctx = mx.current_context()
    @test ctx.device_type == mx.CPU
    @test ctx.device_id   == 42
  end

  function f()
    ctx = mx.current_context()
    @test ctx.device_type == mx.GPU
    @test ctx.device_id   == 123
  end

  @context mx.GPU 123 begin
    f()
  end

  @context mx.GPU begin
    ctx = mx.current_context()
    @test ctx.device_type == mx.GPU
    @test ctx.device_id   == 0
  end

  @context mx.CPU begin
    ctx = mx.current_context()
    @test ctx.device_type == mx.CPU
    @test ctx.device_id   == 0
  end

  @info "Context::@gpu"
  @gpu 123 f()
  @gpu begin
    ctx = mx.current_context()
    @test ctx.device_type == mx.GPU
    @test ctx.device_id   == 0
  end
  let n = 321
    @gpu n begin
      ctx = mx.current_context()
      @test ctx.device_type == mx.GPU
      @test ctx.device_id   == 321
    end
  end

  @info "Context::@cpu"
  @cpu 123 begin
    ctx = mx.current_context()
    @test ctx.device_type == mx.CPU
    @test ctx.device_id   == 123
  end
  @cpu begin
    ctx = mx.current_context()
    @test ctx.device_type == mx.CPU
    @test ctx.device_id   == 0
  end
  let n = 321
    @cpu n begin
      ctx = mx.current_context()
      @test ctx.device_type == mx.CPU
      @test ctx.device_id   == 321
    end
  end
end

@testset "Context Test" begin
  test_num_gpus()
  test_context_macro()
end

end  # module TestContext


text:
using Expectations, Distributions
using Test, Random, Statistics, LinearAlgebra

@testset "Iterable distributions" begin include("iterable.jl") end
@testset "Mixture models" begin include("mixturemodels.jl") end


text:
fnames = [
          "area",
          "ftc",

          "substitution",
          "integration_by_parts",
          "partial_fractions", # XX add in trig integrals (cos()sin() stuff? mx or ^m... XXX
          "improper_integrals", ##

          "mean_value_theorem",
          "area_between_curves",
          "center_of_mass",
          "volumes_slice",
          #"volumes_shell", ## XXX add this in if needed, but not really that excited to now XXX
          "arc_length",
          "surface_area"
]

function process_file(nm, twice=false)
    include("$nm.jl")
    mmd_to_md("$nm.mmd")
    markdownToHTML("$nm.md")
    twice && markdownToHTML("$nm.md")
end

process_files(twice=false) = [process_file(nm, twice) for nm in fnames]

"""
## TODO integrals

* add in volumes shell???
* mean value theorem is light?
* could add surface area problems

"""


text:
module Raft
using ProtoBuf

include("./common_types.jl")

include("./leader.jl")
include("./candidate.jl")
include("./follower.jl")

end


text:
#####
##### Schema
#####

const ALLOWED_SCHEMA_NAME_CHARACTERS = Char['-', '.', 'a':'z'..., '0':'9'...]

"""
    Legolas.is_valid_schema_name(x::AbstractString)

Return `true` if `x` is a valid schema name, return `false` otherwise.

Valid schema names are lowercase, alphanumeric, and may contain hyphens or periods.
"""
is_valid_schema_name(x::AbstractString) = all(i -> i in ALLOWED_SCHEMA_NAME_CHARACTERS, x)

"""
    Legolas.Schema{name,version}

A type representing the schema of a [`Legolas.Row`](@ref). The `name` (a `Symbol`) and `version` (an `Integer`)
are surfaced as type parameters, allowing them to be utilized for dispatch.

For more details and examples, please see `Legolas.jl/examples/tour.jl` and the "Tips for Schema Authors"
section of the Legolas.jl documentation.

See also: [`schema_name`](@ref), [`schema_version`](@ref), [`schema_parent`](@ref)
"""
struct Schema{name,version} end

"""
    Legolas.Schema(name::AbstractString, version::Integer)

Return `Legolas.Schema{Symbol(name),version}()`. This constructor will throw an `ArgumentError` if `name` is
not a valid schema name.

Prefer using this constructor over `Legolas.Schema{Symbol(name),version}()` directly.
"""
function Schema(name::AbstractString, version::Integer)
    version >= 0 || throw(ArgumentError("`Legolas.Schema` version must be non-negative, recieved: $version"))
    is_valid_schema_name(name) || throw(ArgumentError("argument is not a valid `Legolas.Schema` name: \"$name\""))
    return Schema{Symbol(name),version}()
end

"""
    Legolas.Schema(s::AbstractString)

Return `Legolas.Schema(name, n)` where `s` is a valid schema identifier of the form `"name@n"`.

`s` may also be a fully qualified schema identifier of the form `"name@n>...>..."`.
"""
function Schema(s::AbstractString)
    x = split(first(split(s, '>', limit=2)), '@')
    if length(x) == 2
        name, version = x
        version = tryparse(Int, version)
        version isa Int && return Schema(name, version)
    end
    throw(ArgumentError("argument is not a valid `Legolas.Schema` string: \"$s\""))
end

"""
    schema_name(::Type{<:Legolas.Schema{name}})
    schema_name(::Legolas.Schema{name})

Return `name`.
"""
@inline schema_name(::Type{<:Schema{name}}) where {name} = name
@inline schema_name(schema::Schema) = schema_name(typeof(schema))

"""
    schema_version(::Type{Legolas.Schema{name,version}})
    schema_version(::Legolas.Schema{name,version})

Return `version`.
"""
@inline schema_version(::Type{<:Schema{name,version}}) where {name,version} = version
@inline schema_version(schema::Schema) = schema_version(typeof(schema))

"""
    schema_parent(::Type{Legolas.Schema{name,version}})
    schema_parent(::Legolas.Schema{name,version})

Return the `Legolas.Schema` instance that corresponds to the parent of the given `Legolas.Schema`.
"""
@inline schema_parent(::Type{<:Schema}) = nothing
@inline schema_parent(schema::Schema) = schema_parent(typeof(schema))

Base.show(io::IO, schema::Schema) = print(io, "Schema(\"$(schema_name(schema))@$(schema_version(schema))\")")

#####
##### methods overloaded by `@row`
#####

struct UnknownSchemaError <: Exception
    schema::Legolas.Schema
end

function Base.showerror(io::IO, e::UnknownSchemaError)
    print(io, """
              encountered unknown `Legolas.Schema` type: $(e.schema)

              This generally indicates that this schema has not been defined (i.e.
              the schema's corresponding `@row` statement has not been executed) in
              the current Julia session.

              In practice, this can arise if you try to read a Legolas table with a
              prescribed schema, but haven't actually loaded the schema definition
              (or commonly, haven't loaded the dependency that contains the schema
              definition - check the versions of loaded packages/modules to confirm
              your environment is as expected).

              Note that if you're in this particular situation, you can still load
              the raw Arrow table as-is via `Arrow.Table`.
              """)
    return nothing
end

# Note that there exist very clean generic implementations of `transform`/`validate`:
#
#    function transform(schema::Schema; fields...)
#        parent = schema_parent(schema)
#        parent isa Schema && (fields = transform(parent; fields...))
#        return _transform(schema; fields...)
#    end
#
#    function validate(tables_schema::Tables.Schema, legolas_schema::Schema)
#        parent = schema_parent(legolas_schema)
#        parent isa Schema && validate(parent, tables_schema)
#        _validate(tables_schema, legolas_schema)
#        return nothing
#    end
#
# However, basic benchmarking demonstrates that the above versions can allocate
# unnecessarily for schemas with a few ancestors, while the "hardcoded" versions
# generated by the current implementation of the `@row` macro (see below) do not.

transform(s::Legolas.Schema; fields...) = throw(UnknownSchemaError(s))

function _transform end

"""
    Legolas.validate(tables_schema::Tables.Schema, legolas_schema::Legolas.Schema)

Throws an `ArgumentError` if `tables_schema` does comply with `legolas_schema`, otherwise
returns `nothing`.

Specifically, `tables_schema` is considered to comply with `legolas_schema` if:

- every non-`>:Missing` field required by `legolas_schema` is present in `tables_schema`.
- `T <: S` for each field `f::T` in `tables_schema` that matches a required `legolas_schema` field `f::S`.
"""
validate(::Tables.Schema, s::Legolas.Schema) = throw(UnknownSchemaError(s))

function _validate end

function validate_expected_field(schema::Tables.Schema, name::Symbol, ::Type{T}) where {T}
    i = findfirst(==(name), schema.names)
    if isnothing(i)
        Missing <: T || throw(ArgumentError("could not find expected field `$name::$T` in $schema"))
    else
        schema.types[i] <: T || throw(ArgumentError("field `$name` has unexpected type; expected <:$T, found $(schema.types[i])"))
    end
    return nothing
end

"""
    schema_qualified_string(::Legolas.Schema{name,version})

Return this `Legolas.Schema`'s fully qualified schema identifier string. This string is
serialized as the `\"$LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY\"`` field value in table
metadata for table written via [`Legolas.write`](@ref).
"""
schema_qualified_string(s::Legolas.Schema) = throw(UnknownSchemaError(s))

#####
##### Row
#####

"""
    Legolas.Row(schema::Schema; fields...)
    Legolas.Row(schema::Schema, row)

Return a `Legolas.Row <: Tables.AbstractRow` instance whose fields are the provided `fields`
(or the fields of `row`) validated/transformed in accordance with provided `schema`.

For more details and examples, please see `Legolas.jl/examples/tour.jl`.
"""
struct Row{S<:Schema,F} <: Tables.AbstractRow
    schema::S
    fields::F
    function Row(schema::Schema; fields...)
        fields = transform(schema; fields...)
        return new{typeof(schema),typeof(fields)}(schema, fields)
    end
end

Row{S}(args...; kwargs...) where {S} = Row(S(), args...; kwargs...)

Row(schema::Schema, fields) = Row(schema, NamedTuple(Tables.Row(fields)))
Row(schema::Schema, fields::Row) = Row(schema, getfield(fields, :fields))
Row(schema::Schema, fields::NamedTuple) = Row(schema; fields...)

Base.propertynames(row::Row) = propertynames(getfield(row, :fields))
Base.getproperty(row::Row, name::Symbol) = getproperty(getfield(row, :fields), name)

Tables.getcolumn(row::Row, i::Int) = Tables.getcolumn(getfield(row, :fields), i)
Tables.getcolumn(row::Row, nm::Symbol) = Tables.getcolumn(getfield(row, :fields), nm)
Tables.columnnames(row::Row) = Tables.columnnames(getfield(row, :fields))

Base.:(==)(a::Row, b::Row) = getfield(a, :schema) == getfield(b, :schema) && getfield(a, :fields) == getfield(b, :fields)
Base.isequal(a::Row, b::Row) = isequal(getfield(a, :schema), getfield(b, :schema)) && isequal(getfield(a, :fields), getfield(b, :fields))
Base.hash(a::Row, h::UInt) = hash(Row, hash(getfield(a, :schema), hash(getfield(a, :fields), h)))

function Base.show(io::IO, row::Row)
    print(io, "Row($(getfield(row, :schema)), ")
    show(io, getfield(row, :fields))
    print(io, ")")
    return nothing
end

function _parse_schema_expr(x)
    if x isa Expr && x.head == :call && x.args[1] == :> && length(x.args) == 3
        child, _ = _parse_schema_expr(x.args[2])
        parent, _ = _parse_schema_expr(x.args[3])
        return child, parent
    end
    return nothing, nothing
end

_parse_schema_expr(str::AbstractString) = Schema(str), nothing

"""
    @row("name@version", field_expressions...)
    @row("name@version" > "parent_name@parent_version", field_expressions...)

Define a new `Legolas.Schema{name,version}` whose required fields are specified by `field_expressions`.
Returns `Legolas.Row{Legolas.Schema{name,version}}` which can be conveniently aliased to the caller's
preferred binding for a row constructor associated with `Legolas.Schema{name,version}`.

Each element of `field_expression` defines a required field for `Legolas.Schema{name,version}`, and is
an expression of the form `field::F = rhs` where:

- `field` is the corresponding field's name
- `::F` denotes the field's type constraint (if elided, defaults to `::Any`).
- `rhs` is the expression which produces `field::F` (if elided, defaults to `field`).

As implied above, the following alternative forms are also allowed:

- `field::F` (interpreted as `field::F = field`)
- `field = rhs` (interpreted as `field::Any = rhs`)
- `field` (interpreted as `field::Any = field`)

For more details and examples, please see `Legolas.jl/examples/tour.jl` and the "Tips for Schema Authors"
section of the Legolas.jl documentation.
"""
macro row(schema_expr, fields...)
    schema, parent = _parse_schema_expr(schema_expr)
    isnothing(schema) && throw(ArgumentError("`@row` schema argument must be of the form `\"name@X\"` or `\"name@X\" > \"parent@Y\"`. Received: $schema_expr"))
    fields = map(fields) do f
        original_f = f
        f isa Symbol && (f = Expr(:(::), f, :Any))
        f.head == :(::) && (f = Expr(:(=), f, f.args[1]))
        f.head == :(=) && f.args[1] isa Symbol && (f.args[1] = Expr(:(::), f.args[1], :Any))
        f.head == :(=) && f.args[1].head == :(::) || throw(ArgumentError("malformed `@row` field expression: $original_f"))
        return f
    end
    validate_fields = map(fields) do f
        name, type = f.args[1].args
        return :(validate_expected_field(tables_schema, $(Base.Meta.quot(name)), $(esc(type))))
    end
    field_names = [esc(f.args[1].args[1]) for f in fields]
    schema_type = Base.Meta.quot(typeof(schema))
    quoted_parent = Base.Meta.quot(parent)
    schema_qualified_string = string(schema_name(schema), '@', schema_version(schema))
    parent_transform = nothing
    parent_validate = nothing
    if !isnothing(parent)
        schema_qualified_string = :(string($schema_qualified_string, '>', Legolas.schema_qualified_string($quoted_parent)))
        parent_transform = :(fields = transform($quoted_parent; fields...))
        parent_validate = :(validate(tables_schema, $quoted_parent))
    end
    return quote
        Legolas.schema_qualified_string(::$schema_type) = $schema_qualified_string

        Legolas.schema_parent(::Type{<:$schema_type}) = $quoted_parent

        function Legolas._transform(::$schema_type; $([Expr(:kw, f, :missing) for f in field_names]...), other...)
            $(map(esc, fields)...)
            return (; $([Expr(:kw, f, f) for f in field_names]...), other...)
        end

        function Legolas._validate(tables_schema::Tables.Schema, legolas_schema::$schema_type)
            $(validate_fields...)
            return nothing
        end

        function Legolas.transform(schema::$schema_type; fields...)
            $parent_transform
            return _transform(schema; fields...)
        end

        function Legolas.validate(tables_schema::Tables.Schema, legolas_schema::$schema_type)
            $parent_validate
            return _validate(tables_schema, legolas_schema)
        end

        Legolas.Row{$schema_type}
    end
end

# Support (de)serialization as an Arrow column value via Arrow.ArrowTypes overloads.
#
# Note that this only really works in relatively simple cases; rely on this at your own peril.
# See https://github.com/JuliaData/Arrow.jl/issues/230 for more details.
#
# Note also that the limited support here that DOES work participates in SemVer,
# e.g. if we break this in future Legolas versions we should treat it as a breaking
# change and bump version numbers accordingly.

const LEGOLAS_ROW_ARROW_NAME = Symbol("JuliaLang.Legolas.Row")
Arrow.ArrowTypes.arrowname(::Type{<:Legolas.Row}) = LEGOLAS_ROW_ARROW_NAME
Arrow.ArrowTypes.ArrowType(::Type{Legolas.Row{_,F}}) where {_,F} = Tuple{String,Int,F}
Arrow.ArrowTypes.toarrow(row::Legolas.Row{S}) where {S} = (String(Legolas.schema_name(S)), Legolas.schema_version(S), getfield(row, :fields))
Arrow.ArrowTypes.JuliaType(::Val{LEGOLAS_ROW_ARROW_NAME}, ::Any) = Legolas.Row
Arrow.ArrowTypes.fromarrow(::Type{<:Legolas.Row}, name, version, fields) = Legolas.Row(Legolas.Schema(name, version), fields)


text:
# Automatic differentiation on sn.jl:
include("transit_poly.jl")

using ForwardDiff
using DiffResults

function transit_poly_grad(r::T,b::T,u_n::Array{T,1}) where {T <: Real}
  # Computes the derivative of transit_poly(r,b,u_n) with respect to r, b, u_n.
  # Create a vector for use with ForwardDiff:
  x=[r;b;u_n]
  # Now, define a wrapper of transit_poly for use with ForwardDiff:
  function diff_transit_poly(x::Array{T,1}) where {T <: Real}
  # x should be a n+2-element vector with values [r; b; u_n]
  return transit_poly(x[1],x[2],x[3:length(x)])
  end
  dfdp = zeros(typeof(r),length(u_n)+2)
  # Set up a type to store s_n and it's Jacobian with respect to x:
  out = DiffResults.GradientResult(x)
  # Compute the Jacobian (and value):
  out = ForwardDiff.gradient!(out,diff_transit_poly,x)
  # Place the value in the s_n vector:
  flux = DiffResults.value(out)
  # And, place the Jacobian in an array:
  dfdp = DiffResults.gradient(out)
return flux,dfdp
end

function transit_poly_grad_c(r::T,b::T,c_n::Array{T,1}) where {T <: Real}
  # Computes the derivative of transit_poly_c(r,b,c_n) with respect to r, b, c_n.
  # Create a vector for use with ForwardDiff:
  x=[r;b;c_n]
  # Now, define a wrapper of transit_poly for use with ForwardDiff:
  function diff_transit_poly_c(x::Array{T,1}) where {T <: Real}
  # x should be a n+2-element vector with values [r; b; c_n]
  return transit_poly_c(x[1],x[2],x[3:length(x)])
  end
  dfdp = zeros(typeof(r),length(c_n)+2)
  # Set up a type to store s_n and it's Jacobian with respect to x:
  out = DiffResults.GradientResult(x)
  # Compute the Jacobian (and value):
  out = ForwardDiff.gradient!(out,diff_transit_poly_c,x)
  # Place the value in the s_n vector:
  flux = DiffResults.value(out)
  # And, place the Jacobian in an array:
  dfdp = DiffResults.gradient(out)
return flux,dfdp
end

# Here are some example lines of how to use it:
#r = 100.0; b=100.5

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
using OmegaCore, Test, Distributions

function test()
  function f(ω)
    a = (1 ~ Normal(0, 1))(ω)
    b = (2 ~ Normal(0, 1))(ω)
    a + b
  end

  function g(ω)
    x = 0.0
    for i = 1:100
      x += (i <| f)(ω)
    end
    x
  end

  function h(ω)
    α = 1.0
    for i = 1:10
      α *= (i <| f)(ω)
    end
    α
  end

  i(ω) = (1 <| g)(ω) + (2 <| h)(ω)

  ω = defω()
  i(ω)
  @test length(keys(ω)) == 220
end

function normal_normals_test()
  μ = 1 ~ Normal(0, 1)
  x = 2 ~ Normal(μ, 1)
end

function blr()
  m = 1 ~ Normal(0, 1)
  c = 2 ~ Normal(0, 1)
  n = 100
  data = rand(n)
  f(x, m_, c_) = m_*x + c_
  obs(id, ω) = f(x, m_(ω), c_(ω)) + id ~ Normal(ω, 0, 0.1)
  obs(ω) = [f(x, m_(ω), c_(ω)) + i ~ Normal(ω, 0, 0.1) for i in data]
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
import .Convex
using .GLPK

"""
    ot_convex(D, P1, P2)

Solve the optimal transport problem using Convex.jl. This function is only available if `using Convex.jl, GLPK`.

# Arguments:
- `D`: Distance matrix
- `P1`: Weight vector 1
- `P2`: Weight vector 2
"""
function ot_convex(D::AbstractMatrix{T}, P1, P2; kwargs...) where T
    n = length(P1)
    @assert all(isfinite, P1) "Got nonfinite P1"
    @assert all(isfinite, P2) "Got nonfinite P2"
    @assert sum(P1) ≈ 1 "Sum(P1) ≠ 1: $(sum(P1))"
    @assert sum(P2) ≈ 1 "Sum(P2) ≠ 1: $(sum(P2))"
    γ = Convex.Variable(n^2)
    Γ = reshape(γ,n,n)
    con1 = sum(Γ, dims=1) == P2'
    con2 = sum(Γ, dims=2) == P1
    problem = Convex.minimize(γ'vec(D), γ >= 0, con1, con2)
    Convex.solve!(problem, GLPK.Optimizer)
    if Int(problem.status) != 1
        @error problem.status
    end
    α, β = -con2.dual[:], -con1.dual[:]
    α .-= mean(α)
    β .-= mean(β)
    reshape(γ.value,n,n), α, β
end


text:
module HM
export HMT, TVar, mk_tcstate, UN, VerboseUN
export Var, Nom, Bound, Var, Arrow, App, Tup, Forall, Implicit
export gTrans, gTransCtx, gCheck
export tvar_of_int, int_of_tvar, IllFormedType
export ⪯, less_than_under_evidences

using MLFS: @typedIO
import MLFS
using MLStyle
import Base

const VerboseUN = Ref(false)

mutable struct UN # unique name
    name :: Symbol
    ts :: UInt64
end
const pool = Dict{Tuple{UInt64, UInt64}, UN}()

function UN(n::Symbol)
    un = UN(n, time_ns())
    pool[(un.ts, objectid(un))] = un
    un
end

MLFS.TypedIO.toVec(::Type{UN}, un::UN) = [
    MLFS.TypedIO.toVec(Symbol, un.name), string(objectid(un)), string(un.ts)]
function MLFS.TypedIO.fromVec(::Type{UN}, v::Vector)
    n = Symbol(v[1])
    objectid = parse(UInt64, v[2])
    ts = parse(UInt64, v[3])
    key = (ts, objectid)

    Some(get!(pool, key) do
        un = UN(n, ts)
    end)

end

const noImplicits = Val(false)

Base.show(io::IO, unique_name::UN) =
    if VerboseUN[]
        Base.print(io, unique_name.name, "@", objectid(unique_name))
    else
        Base.print(io, unique_name.name)
    end

@data HMT begin
    Var(id::UInt)
    Nom(n::Symbol)
    Bound(n::UN)
    App(f::HMT, arg::HMT)
    Arrow(from::HMT, to::HMT)
    Tup(Vector{HMT})
    Forall(Vector{UN}, HMT)

    Implicit(HMT)
end

⊤ = let a = UN(:a)
    Forall([a], Bound(a))
end

@typedIO Nom(Symbol)

@typedIO HMT = begin
    Nom(Symbol)
    Bound(UN)
    App(HMT, HMT)
    Arrow(HMT, HMT)
    Tup(Vector{HMT})
    Forall(Vector{UN}, HMT)

    Implicit(HMT)
end

need_parens(hmt::HMT) =
    @match hmt begin
        ::App || ::Arrow || ::Forall => true
        _ => false
    end

with_parens(f::Function, io::IO, need_parens::Bool) =
    if need_parens
        print(io, '(')
        f()
        print(io, ')')
    else
        f()
    end

Base.show(io::IO, hmt::HMT) =
    @switch hmt begin
    @case Var(v)
        print(io, '\'')
        print(io, v)
        return
    @case Nom(n)
        print(io, n)
        return

    @case Bound(n)
        print(io, n)
        return

    @case App(f, a)
        Base.show(io, f)
        print(io, " ")
        with_parens(io, need_parens(a)) do
            Base.show(io, a)
        end
        return
    @case Arrow(a, r)
        with_parens(io, need_parens(a)) do
            Base.show(io, a)
        end
        print(io, "->")
        Base.show(io, r)
        return
    @case Tup(xs)
        print(io, "{")
        print(io, join(repr.(xs), ", "))
        print(io, "}")
        return
    @case Forall(ns, t)
        print(io, "forall ")
        print(io, join(string.(ns), " "), ".")
        Base.show(io, t)
        return
    @case Implicit(t)
        print(io, "implicit[$t]")
    end

function gTrans(self::Function, root::HMT)
    @match root begin
        Var(_) || Nom(_) || Bound(_) => root
        App(f, arg) => App(self(f), self(arg))
        Arrow(arg, ret) => Arrow(self(arg), self(ret))
        Tup(xs) => Tup(HMT[self(x) for x in xs])
        Forall(ns, p) => Forall(ns, self(p))
        Implicit(x) => Implicit(self(x))
    end
end

function gTransCtx(self′::Function, ctx::Ctx, root::HMT) where Ctx
    @inline self(root::HMT) = self′(ctx, root)
    @match root begin
        Var(_) || Nom(_) || Bound(_) => root
        App(f, arg) => App(self(f), self(arg))
        Arrow(arg, ret) => Arrow(self(arg), self(ret))
        Tup(xs) => Tup(HMT[self(x) for x in xs])
        Forall(ns, p) => Forall(ns, self(p))
        Implicit(x) => Implicit(self(x))
    end
end

function gCheck(self::Function, root::HMT)
    @match root begin
        Var(_) || Nom(_) || Bound(_) => true
        App(f, arg) => self(f) && self(arg)
        Arrow(arg, ret) => self(arg) && self(ret)
        Tup(xs) => all(self(x) for x in xs)
        Forall(_, p) => self(p)
        Implicit(x) => self(x)
    end
end

function fresh(substitutions::Dict{UN, V}, root::HMT) where V <: HMT
    substitutions = Dict{Bound, V}(Bound(k) => v for (k, v) in substitutions)
    fresh(substitutions, root)
end

function fresh(substitutions::Dict{Bound, V}, root::HMT) where V <: HMT
    function freshRec(root::HMT)
        @match root, substitutions begin
        (Bound(_), Dict(root => var)) => var
        # (Forall(ns, p), _) => let restore = Pair{Bound, HMT}[]
        #         for n in ns
        #             k = Bound(n)
        #             v = pop!(substitutions, k, nothing)
        #             if v !== nothing
        #                 push!(restore, k => v)
        #             end
        #         end
        #         r = isempty(substitutions) ? p : gTrans(freshRec, p)
        #         for (k, v) in restore
        #             substitutions[k] = v
        #         end
        #         Forall(ns, r)
        #     end
        (t, _) => gTrans(freshRec, t)
        end
    end
    freshRec(root)
end

function subst(substitutions::Dict{A, B}, root::HMT) where {A <: HMT, B <: HMT}
    @match substitutions begin
        Dict(root => replaced) => replaced
        _ => gTransCtx(subst, substitutions, root)
    end
end
struct IllFormedType <: Exception
    msg::String
end

function tvar_of_int(i::Integer)
    Var(i)
end

function int_of_tvar(x::Var)
    x.id
end

function un_of_bound(x::Bound)
    x.n
end

function bound_of_un(x::UN)
    Bound(x)
end

function mk_tcstate(tctx::Vector{HMT})
    bound_links = Dict{UN, Set{UInt}}()

    function unlink(un::Union{UN, Bound})
        unlink(un) do _::UInt
            nothing
        end
    end

    function unlink(f::Function, b::Bound)
        unlink(f, un_of_bound(b))
    end

    function unlink(f::Function, un::UN)
        rels = pop!(bound_links, un, nothing)
        rels === nothing && return
        for rel in rels
            tctx[rel] = Var(rel)
            f(rel)
        end
        empty!(rels)
        nothing
    end

    function new_tvar()::HMT
        vid = UInt(length(tctx) + 1)
        tvar = tvar_of_int(vid)
        push!(tctx, tvar)
        tvar
    end

    function new_tvar_id()::UInt
        vid = UInt(length(tctx) + 1)
        tvar = tvar_of_int(vid)
        push!(tctx, tvar)
        vid
    end

    function occur_in(i::UInt, ty::HMT)
        @switch ty begin
            @case Var(i′) && if i′ === i end
                return false
            @case _
                visit_func(x) =
                    @match x begin
                        Var(i′) && if i′ === i end => false
                        _ => gCheck(visit_func, x)
                    end
                return !visit_func(ty)
        end
    end

    function link!(typeref::UInt, b::Bound)
        un = un_of_bound(b)
        links =
            get!(bound_links, un) do
                Set{UInt}(UInt[typeref])
            end
        typeref in links || begin
            push!(links, typeref)
            tctx[typeref] = b
        end
        b
    end

    function link!(typeref::UInt, t::HMT)
        tctx[typeref] = t
        t
    end

    function prune(a::HMT)
        @match a begin
            Var(i) =>
                @match tctx[i] begin
                    Var(i′) && if i′ === i end => a
                    t => link!(i, prune(t))
                end
            a => gTrans(prune, a)
        end
    end

    function pruneWithVarCheck(f::Function, a::HMT)
        function prune(a::HMT)
            v = @match a begin
                Var(i) =>
                    @match tctx[i] begin
                        Var(i′) && if i′ === i end => a

                        t => link!(i, prune(t))
                    end
                a => gTrans(prune, a)
            end
            if v isa Var
                f(v)
            end
            v
        end
        prune(a)
    end

    function instantiate(x::HMT)
        @match x begin
            Forall(ns, ty) =>
                let freshmap = Dict{Bound, Var}(Bound(a) => new_tvar() for a in ns)
                    subst(freshmap, ty)
                end
            _ => x
        end
    end

    function unifyImplicits!(lhs::HMT, rhs::HMT, implicits::Vector{HMT})
        lhs = prune(lhs)
        rhs = prune(rhs)
        lhs === rhs && return true

        @match lhs, rhs begin
            (Forall(_, lhs), _) => unifyImplicits!(lhs, rhs, implicits)
            (Var(_), _) ||
            (_, Var(_)) =>
                unifyINST(lhs, rhs)

            (_, Forall()) =>
                begin
                    unifyImplicits!(lhs, instantiate(rhs), implicits)
                end

            (_, Arrow(Implicit(im), rhs)) =>
                begin # for debugger
                    push!(implicits, im)
                    unifyImplicits!(lhs, rhs, implicits)
                end

            _ => unifyINST(lhs, rhs)
        end
    end

    function unifyINST(lhs::HMT, rhs::HMT)
        lhs = prune(lhs)
        rhs = prune(rhs)
        lhs === rhs && return true
        # @info :INST lhs rhs

        @match lhs, rhs begin
            (Forall(_, lhs), _) => unifyINST(lhs, rhs)

            (Var(i && ai), b) ||
            (b, Var(i && ai)) =>
                if occur_in(ai, b)
                    throw(IllFormedType("a = a -> b"))
                else
                    link!(i, b)
                    true
                end

            (_, Forall()) =>
                begin # for debugger
                    unifyINST(lhs, instantiate(rhs))
                end

            (_, Implicit(rhs)) => unifyINST(lhs, rhs)
            (Implicit(lhs), _) => unifyINST(lhs, rhs)
            (Arrow(a1, r1), Arrow(a2, r2)) =>
                begin # for debugger
                    unifyINST(a2, a1) && unifyINST(r1, r2)
                end

            (App(f1, a1), App(f2, a2)) =>
                begin # for debugger
                    unifyINST(f1, f2) && unifyINST(a1, a2)
                end

            (Tup(xs1), Tup(xs2)) =>
                all(zip(xs1, xs2)) do (lhs, rhs)
                    unifyINST(lhs, rhs)
                end
            _ => false
        end
    end

    function unify(lhs::HMT, rhs::HMT)
        lhs = prune(lhs)
        rhs = prune(rhs)
        lhs === rhs && return true

        @match lhs, rhs begin
            (Forall(ns1, p1), Forall(ns2, p2)) =>
                (begin
                    N1, N2 = length(ns1), length(ns2)
                    N1 !== N2 && return false

                    subst1 = Dict{UN, Var}(a => new_tvar() for a in ns1)
                    generic = p2
                    unify(fresh(subst1, p1), generic) || return false

                    generic = prune(generic)
                    remap2 = Bound[Bound(n) for n in ns2]
                    remap1 = Dict{Bound, Bound}()

                    for (k, v) in subst1
                        v = prune(v)
                        v in remap2 || return false
                        haskey(remap1, v) && return false
                        remap1[v] = Bound(k)
                    end

                    backmap1 = Dict{Bound, Bound}()
                    backmap2 = Dict{Bound, Bound}()

                    for un in ns2
                        unlink(un) do i::UInt
                            backmap1[Var(i)] = remap1[Bound(un)]
                            backmap2[Var(i)] = Bound(un)
                        end
                    end

                    unify(subst(backmap2, p2), p2) &&
                    unify(subst(backmap1, p1), p1)
                end)

            (Var(i && ai), b) ||
            (b, Var(i && ai)) =>

                if occur_in(ai, b)
                    throw(IllFormedType("a = a -> b"))
                else
                    link!(i, b)
                    true
                end

            (Arrow(a1, r1), Arrow(a2, r2)) =>
                unify(a1, a2) && unify(r1, r2)

            (App(f1, a1), App(f2, a2)) =>
                unify(f1, f2) && unify(a1, a2)

            (Tup(xs1), Tup(xs2)) =>
                all(zip(xs1, xs2)) do (lhs, rhs)
                    unify(lhs, rhs)
                end
            _ => false
        end
    end

    (unifyINST = unifyINST,
     unifyImplicits! = unifyImplicits!,
     unify = unify,
     tctx = tctx,
     new_tvar = new_tvar,
     new_tvar_id = new_tvar_id,
     unlink = unlink,
     occur_in = occur_in,
     instantiate = instantiate,
     prune = prune,
     pruneWithVarCheck = pruneWithVarCheck)

end

(lhs::HMT) ⪯ (rhs::HMT)  = begin
    small_tc = mk_tcstate(HMT[])
    subst_table = Dict{UInt, Var}()
    function subst(root::HMT)
        @match root begin
            Var(i) =>
                get!(subst_table, i) do
                    small_tc.new_tvar()
                end

            _ => gTrans(subst, root)
        end
    end
    small_tc.unifyImplicits(subst(lhs), subst(rhs), HMT[])
end

function less_than_under_evidences(lhs::HMT, rhs ::HMT)
    evidences = HMT[]
    small_tc = mk_tcstate(HMT[])
    subst_table = Dict{UInt, Var}()
    function subst(root::HMT)
        @match root begin
            Var(i) =>
                get!(subst_table, i) do
                    small_tc.new_tvar()
                end

            _ => gTrans(subst, root)
        end
    end
    small_tc.unifyImplicits!(subst(lhs), subst(rhs), evidences), evidences
end

end


text:
module SatelliteAnalysis

using Dates
using LinearAlgebra
using PrettyTables
using ReferenceFrameRotations
using StaticArrays
using SatelliteToolbox
using Statistics
using TerminalPager

################################################################################
#                                   Includes
################################################################################

include("./beta_angle.jl")
include("./eclipse_time.jl")
include("./lighting_condition.jl")
include("./misc/find_crossing.jl")

end # module


text:
## CHull.jl
## (c) 2013 David Al van Leeuwen
## A Julia wrapper around a python wrapper around the qhull Convex Hull library

## This code is licensed under the GNU General Public License, version 2
## See the file LICENSE in this distribution

module QHull

export Chull, chull

using PyCall
const spatial = PyNULL()

function __init__()
    copy!(spatial, pyimport_conda("scipy.spatial", "scipy"))
end

mutable struct Chull{T<:Real}
    points::Matrix{T}
    vertices::Vector{Int}
    simplices::Vector{Vector{Int}}
    facets::Matrix{T}
    area::Float64
    volume::Float64
end

## helper for base-0 / base-1 difference
incone(x) = for i in 1:length(x)
    x[i] += 1
end

function chull(x::Matrix{T}) where T<:Real
    py = spatial.ConvexHull(x)
    points = convert(Matrix{T}, py."points")
    vertices = convert(Vector{Int}, py."vertices")
    incone(vertices)
    simplices = convert(Vector{Vector{Int}}, py."simplices")
    for simplex in simplices
        incone(simplex)
    end
    facets = convert(Matrix{T}, py."equations")
    area = convert(Float64, py."area")
    volume = convert(Float64, py."volume")
    Chull(points, vertices, simplices, facets, area, volume)
end

include("polyhedron.jl")

function Base.show(io::IO, ::MIME"text/plain", ch::Chull)
    println(io, string("Convex Hull of ", size(ch.points, 1), " points in ", size(ch.points, 2), " dimensions"))
    println(io, "Hull segment vertex indices:")
    println(io, ch.vertices)
    println(io, "Points on convex hull in original order:\n")
    println(io, ch.points[sort(ch.vertices[:, 1]), :])
end

end


text:
# Test iswithinrectangle function

using RANSAC: iswithinrectangle

using RegionTrees: HyperRectangle, vertices
using RegionTrees: Cell, findleaf, adaptivesampling!

rectangle = HyperRectangle(SVector(0.0,0,0), SVector(1.0,1,1))

@testset "corner points" begin
    # false
    v000 = SVector(0.0,0,0)
    v100 = SVector(1.0,0,0)
    v010 = SVector(0.0,1,0)
    v001 = SVector(0.0,0,1)
    v110 = SVector(1.0,1,0)
    v101 = SVector(1.0,0,1)
    v011 = SVector(0.0,1,1)

    # true
    v111 = SVector(1.0,1,1)

    @test iswithinrectangle(rectangle, v000) == false
    @test iswithinrectangle(rectangle, v100) == false
    @test iswithinrectangle(rectangle, v010) == false
    @test iswithinrectangle(rectangle, v001) == false
    @test iswithinrectangle(rectangle, v110) == false
    @test iswithinrectangle(rectangle, v101) == false
    @test iswithinrectangle(rectangle, v011) == false

    @test iswithinrectangle(rectangle, v111) == true
end

@testset "edge midpoints" begin
    # false
    v500 = SVector(0.5,0,0)
    v501 = SVector(0.5,0,1)
    v050 = SVector(0,0.5,0)
    v051 = SVector(0,0.5,1)
    v150 = SVector(1,0.5,0)
    v510 = SVector(0.5,1,0)
    v105 = SVector(1,0,0.5)
    v005 = SVector(0,0,0.5)
    v015 = SVector(0,1,0.5)

    # true
    v151 = SVector(1,0.5,1)
    v511 = SVector(0.5,1,1)
    v115 = SVector(1,1,0.5)

    @test iswithinrectangle(rectangle, v500) == false
    @test iswithinrectangle(rectangle, v501) == false
    @test iswithinrectangle(rectangle, v050) == false
    @test iswithinrectangle(rectangle, v051) == false
    @test iswithinrectangle(rectangle, v150) == false
    @test iswithinrectangle(rectangle, v510) == false
    @test iswithinrectangle(rectangle, v105) == false
    @test iswithinrectangle(rectangle, v005) == false
    @test iswithinrectangle(rectangle, v015) == false

    @test iswithinrectangle(rectangle, v151) == true
    @test iswithinrectangle(rectangle, v511) == true
    @test iswithinrectangle(rectangle, v115) == true
end

@testset "face midpoints" begin
    # false
    v550 = SVector(0.5,0.5,0)
    v505 = SVector(0.5,0,0.5)
    v055 = SVector(0,0.5,0.5)

    # true
    v551 = SVector(0.5,0.5,1)
    v515 = SVector(0.5,1,0.5)
    v155 = SVector(1,0.5,0.5)

    @test iswithinrectangle(rectangle, v550) == false
    @test iswithinrectangle(rectangle, v505) == false
    @test iswithinrectangle(rectangle, v055) == false

    @test iswithinrectangle(rectangle, v551) == true
    @test iswithinrectangle(rectangle, v515) == true
    @test iswithinrectangle(rectangle, v155) == true
end

@testset "inside points" begin
    # like facepoints but inside
    v5501 = SVector(0.5,0.5,0.1)
    v5015 = SVector(0.5,0.1,0.5)
    v0155 = SVector(0.1,0.5,0.5)
    v5509 = SVector(0.5,0.5,0.9)
    v5095 = SVector(0.5,0.9,0.5)
    v0955 = SVector(0.9,0.5,0.5)

    vin = [v5501, v5015, v0155, v5509, v5095, v0955]
    for vi in vin
        @test iswithinrectangle(rectangle, vi)
    end
end

@testset "outside points" begin
    # like facepoints but outside
    v5501m = SVector(0.5,0.5,-0.1)
    v5015m = SVector(0.5,-0.1,0.5)
    v0155m = SVector(-0.1,0.5,0.5)
    v5509p = SVector(0.5,0.5,1.1)
    v5095p = SVector(0.5,1.1,0.5)
    v0955p = SVector(1.1,0.5,0.5)

    vout = [v5501m, v5015m, v0155m, v5509p, v5095p, v0955p]
    for vo in vout
        @test iswithinrectangle(rectangle, vo) == false
    end
end

@testset "getnthcell tests" begin
    ps = [SVector(i, j, k)/3 for i in 0:5 for j in 0:5 for k in 0:5]
    pc = RANSACCloud(ps, ps, 1)

    l = findleaf(pc.octree, ps[117])

    @test l.data.depth == 3
    @test RANSAC.getnthcell(l, l.data.depth) == l
    @test RANSAC.getnthcell(l, 3) == l
    @test RANSAC.getnthcell(l, 2) == parent(l)
    @test RANSAC.getnthcell(l, 1) == parent(parent(l))
    @test RANSAC.getnthcell(l, 1) == pc.octree

    @test RANSAC.getnthcell(l, -1) === nothing
    @test RANSAC.getnthcell(l, 0) === nothing
    @test RANSAC.getnthcell(l, 4) === nothing
    @test RANSAC.getnthcell(l, 5) === nothing
end

@testset "octree depth" begin
    ps = [SVector(i, j, k)/3 for i in 0:5 for j in 0:5 for k in 0:5]
    pc = RANSACCloud(ps, ps, 1)
    @test RANSAC.octreedepth(pc.octree) == 3
    @test RANSAC.octreedepth(pc) == 3
end


text:
@testset "$TEST $G" begin

if !@isdefined(is_nonbacktracking)
    #is_nonbacktracking: a predicate that tests if a walk is nonbacktracking.
    #    That is that for no index i walk[i+2]==walk[i]
    function is_nonbacktracking(walk)
        n = length(walk)
        if n < 3
           return true
        end
        for i in 1:n-2
            if walk[i+2] == walk[i]
                return false
            end
        end
        return true
    end

    #test_nbw: check that a walk is nonbacktracking and print it if is isn't.
    #Used only for testing and debugging.
    function test_nbw(g, start, len)
        w = nonbacktracking_randomwalk(g, start, len)
        if is_nonbacktracking(w)
            return true
        else
            print("walk was:\n  $w")
        end
        return false
    end
end

g = PathDiGraph(10, DG)
@test randomwalk(g, 1, 5) == [1:5;]
@test randomwalk(g, 2, 100) == [2:10;]
@test_throws BoundsError randomwalk(g, 20, 20)

g = PathGraph(10, G)
@test self_avoiding_randomwalk(g, 1, 20) == [1:10;]
@test_throws BoundsError self_avoiding_randomwalk(g, 20, 20)

g = PathGraph(10, G)
@test nonbacktracking_randomwalk(g, 1, 20) == [1:10;]
@test_throws BoundsError nonbacktracking_randomwalk(g, 20, 20)

g = digraph(PathGraph(10, G))
@test nonbacktracking_randomwalk(g, 1, 20) == [1:10;]
@test_throws BoundsError nonbacktracking_randomwalk(g, 20, 20)

g = PathDiGraph(10, DG)
@test nonbacktracking_randomwalk(g, 10, 20) == [10]

g = PathDiGraph(10, DG)
@test nonbacktracking_randomwalk(g, 1, 20) == [1:10;]

g = CycleGraph(10, G)
visited = nonbacktracking_randomwalk(g, 1, 20)
@test visited == [1:10; 1:10;] || visited == [1; 10:-1:1; 10:-1:2;]

g = CycleDiGraph(10, DG)
@test nonbacktracking_randomwalk(g, 1, 20) == [1:10; 1:10;]

n = 10
g = CycleGraph(n, G)
for k = 3:n-1
    add_edge!(g, 1, k)
end

for len = 1:3*n
    @test test_nbw(g,1,len)
    @test test_nbw(g,2,len)
end
#test to make sure it works with self loops.
add_edge!(g, 1, 1)
for len = 1:3*n
    @test test_nbw(g,1,len) == true
    @test test_nbw(g,2,len) == true
end

end # testset


text:
module Solution
  import LinearAlgebra
  predictAge(age...) = LinearAlgebra.norm(age) ÷ 2
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
using StableArrays
using StableArrays: StableNumber, unstabilize, logview
using LinearAlgebra
using Test

@testset "StableArray" begin
    a = rand(3,3)
    as = stabilize!(a)
    @test StableArrays.base(as') == StableArrays.base(as)'

    bs = stabilize!(rand(3,3))
    cs = as * bs
    @test cs isa StableArray
    @test maximum(abs, StableArrays.base(cs)) ≈ 1
    @test unstabilize(cs) ≈ unstabilize(as) * unstabilize(bs)

    @testset "zeros" begin
        a = zeros(4,4)
        as = stabilize(a)
        @test all(StableArrays.base(as) .== 0)
        @test unstabilize(as) == a
    end
end

@testset "StableNumber" begin
    a = stabilize!(rand(3,3))
    @test tr(a) isa StableNumber
    @test maximum(a) isa StableNumber
    x = StableNumber(1.0, 10.1)
    y = StableNumber(9.0, 5.1)

    @testset "plus" begin
        z = x + y
        @test z isa StableNumber
        @test StableArrays.exponent(z) == StableArrays.exponent(x)
        @test unstabilize(z) ≈ unstabilize(x) + unstabilize(y)
    end

    @testset "minus" begin
        z = x - y
        @test z isa StableNumber
        @test StableArrays.exponent(z) == StableArrays.exponent(x)
        @test unstabilize(z) ≈ unstabilize(x) - unstabilize(y)
    end

    @testset "times" begin
        z = x * y
        @test z isa StableNumber
        @test StableArrays.exponent(z) == StableArrays.exponent(x) + StableArrays.exponent(y)
        @test StableArrays.base(z) == StableArrays.base(x) * StableArrays.base(y)
        @test log(z) ≈ log(x) + log(y)
    end

    @testset "division" begin
        z = x / y
        @test z isa StableNumber
        @test StableArrays.exponent(z) == StableArrays.exponent(x) - StableArrays.exponent(y)
        @test StableArrays.base(z) == StableArrays.base(x) / StableArrays.base(y)
        @test log(z) ≈ log(x) - log(y)

        z = x / 2
        @test z isa StableNumber
        @test log(z) ≈ log(x) - log(2)
    end

    @testset "less than" begin
        @test 50000 > StableNumber(2, 10) > 0.9
        @test StableNumber(4, 7) > StableNumber(2, 7) > StableNumber(2, 5)
        @test StableNumber(-7, 3) < 0
    end

    @testset "abs" begin
        a = StableNumber(-2,3)
        b = abs(a)
        @test StableArrays.base(b) == 2
        @test StableArrays.exponent(b) == 3
    end

    @testset "isnan" begin
        @test isnan(StableNumber(NaN, 8))
        @test isnan(StableNumber(9, NaN))
        @test isnan(StableNumber(NaN, NaN))
        @test isnan(StableNumber(0, 8) / StableNumber(0, 2))
    end

    @testset "logview" begin
        a = StableNumber(-2, 3)
        @test unstabilize(logview(a)) ≈ unstabilize(a)
    end
end


text:
A = rand(1:4,3,3)

x = fill(1.0, (3))

b = A*x

Asym = A + A'

Apd = A'A

A\b

Atall = rand(3, 2)
display(Atall)
Atall\b

v = randn(3)
rankdef = [v v]

rankdef\b

Ashort = rand(2, 3)
display(Ashort)
Ashort\b[1:2]

max_stars_repo_path: es-es/intro-to-julia-ES/plutos/10. Algebra lineal basica.jl

max_stars_repo_name: miguelraz/JuliaTutorials

max_issues_repo_name: miguelraz/JuliaTutorials

max_forks_repo_name: miguelraz/JuliaTutorials


text:
using DataDeps

register(DataDep(
  "haberman-survival",
  "http://archive.ics.uci.edu/ml/datasets/Haberman%27s+Survival",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/haberman/haberman.data",
  "b4b7a32586a5668f9f4d6dc8be9d1bc8cd4822523affb1f6b5bfc350681ef3e2",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=4,
        feature_indices=1:3,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data