text:
# Compatibility functions that are particular to Symata.

# The following are deprecated in v0.7 to small routines.
symrsearch(s::AbstractString, c::Union{Tuple{Vararg{Char}}, AbstractVector{Char}, Set{Char}},
           pos) = something(findprev((in)(c), s, pos), 0)
symrsearch(s::AbstractString, c::Char,
           pos) = something(findprev((in)(c), s, pos), 0)
#                                              something(findprev((in)(c), s, pos), 0)
symsearch(s::AbstractString, c::Char) = something(findfirst(isequal(c), s), 0)
symsearchindex(s::AbstractString, t::AbstractString) = first(something(findfirst(t, s), 0:-1))
symsearch(s::AbstractString, c::Char, i::Integer) = something(findnext(isequal(c), s, i), 0)
symreplace(s::AbstractString, pat, f) = replace(s, pat => f)
## Warning! The deprecation says this is a drop-in replacement.
## But, it fails in some cases. Eg. in collect_times
symfindin(a, b) = findall((in)(b), a)

# as of at least
# v"0.5.0-dev+3868"
# 'symbol' is deprecated in favor of 'Symbol'
# The compatibility code here works at least as far back as
# v"0.4.6-pre+7"
# const have_new_Symbol =
#     try
#         Symbol("a","b")
#         true
#     catch
#         false
#     end

# if ! have_new_Symbol
#     Symbol(args...) = symbol(args...)
# end

#### readstring

# if ! @isdefined(readstring)
#     function readstring(fname::AbstractString)
#        stream = open(fname, "r")
#        str = readall(stream)
#        close(stream)
#        str
#    end
# end

#### view

# if ! @isdefined(view)
#     view(args...) = slice(args...)
# end

#### bytestring_beforecursor

# for v0.7-beta
# if isdefined(Base.REPL, :bytestring_beforecursor)
#     symata_beforecursor(args...) = Base.REPL.bytestring_beforecursor(args...)
# else
import REPL  # try to avoid warning
symata_beforecursor(args...) = REPL.beforecursor(args...)
#end

### denominator

## We can probably use Compat here

# FIXME: replace for v0.7-beta
# if ! @isdefined(denominator)
#     denominator(args...) = den(args...)
#     symdenominator(args...) = den(args...)
# else
symdenominator(args...) = denominator(args...)
#end

# if ! @isdefined(numerator)
#     numerator(args...) = num(args...)
#     symnumerator(args...) = num(args...)
# else
symnumerator(args...) = numerator(args...)
#end

if isdefined(Base, :invokelatest)
    invokelatest(args...) = Base.invokelatest(args...)
else
    invokelatest(args...) = args[1](args[2:end]...)
end

## The suggested replacement for findin fails
## in coefficient_times, algebra.jl.
## So, we reproduce the version from v0.6.3
## We continue to get deprecation warnings even thought this is
## Not exported to Base or Main. So, we rename it to oldfindin.
"""
    Symata.oldfindin(a, b)

Returns the indices of elements in collection `a` that appear in collection `b`.

# Examples
```jldoctest
julia> a = collect(1:3:15)
5-element Array{Int64,1}:
  1
  4
  7
 10
 13

julia> b = collect(2:4:10)
3-element Array{Int64,1}:
  2
  6
 10

julia> oldfindin(a,b) # 10 is the only common element
1-element Array{Int64,1}:
 4
```
"""
function oldfindin(a, b)
    ind = Array{Int,1}(undef, 0)
    bset = Set(b)
    @inbounds for (i,ai) in enumerate(a)
        ai in bset && push!(ind, i)
    end
    ind
end

## for v0.6, we used Base.syntax_deprecation_warnings
## 1) This function is not present in v0.7
## 2) The deprecated syntax will probably be removed soon.
## So, we make this a noop for the time being
symata_syntax_deprecation_warnings(f, args...) = f()


text:
export PoreData

@doc raw"""
`PoreData` is a collection of physical parameters for coupled geomechanics and flow simulation

- `M`: Biot modulus
- `b`: Biot coefficient
- `ρb`: Bulk density
- `ρf`: Fluid density
- `kp`: Permeability
- `E`: Young modulus
- `ν`: Poisson ratio
- `μ`: Fluid viscosity
- `Pi`: Initial pressure
- `Bf`: formation volume, $B_f=\frac{\rho_{f,0}}{\rho_f}$
- `g`: Gravity acceleration
"""
@with_kw mutable struct PoreData
    M::Float64 = 0.1
    b::Float64 = 1.0
    ρb::Float64 = 2400.
    ρf::Float64 = 1000.
    kp::Float64 = 50e6 * 9.869233e-13 # Darcy
    E::Float64 = 350e6
    ν::Float64 = 0.0
    μ::Float64 = 8.9e-4 # viscosity of water
    Pi::Float64 = 2.125e6
    Bf::Float64 = 1.0
    g::Float64 = 9.807
end

ext: jl


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

"""
    lglowand!(state::SmartsParserState) -> Union{Pair,Nothing}

LogicalLowAnd <- Or (';' Or)*

The argument `func` is a parser function which has a parser state as an argument,
process tokens found in the given text, and returns nothing if no valid tokens were found.
"""
function lglowand!(state::SmartsParserState, func)
    fmls = Set{QueryFormula}()
    fml = lgor!(state, func)
    fml === nothing && return
    while fml !== nothing
        push!(fmls, fml)
        if read(state) == ';'
            forward!(state)
            fml = lgor!(state, func)
        else
            break
        end
    end
    @assert !isempty(fmls) "(lglowand!) invalid AND(;) operation"
    if length(fmls) == 1
        return collect(fmls)[1]
    else
        return QueryFormula(:and, fmls)
    end
end

"""
    lgor!(state::SmartsParserState) -> Union{Pair,Nothing}

Or <- And (',' And)*

The argument `func` is a parser function which has a parser state as an argument,
process tokens found in the given text, and returns nothing if no valid tokens were found.
"""
function lgor!(state::SmartsParserState, func)
    fmls = Set{QueryFormula}()
    fml = lghighand!(state, func)
    fml === nothing && return
    while fml !== nothing
        push!(fmls, fml)
        if read(state) == ','
            forward!(state)
            fml = lghighand!(state, func)
        else
            break
        end
    end
    @assert !isempty(fmls) "(lgor!) invalid OR(,) operation"
    if length(fmls) == 1
        return collect(fmls)[1]
    else
        return QueryFormula(:or, fmls)
    end
end

"""
    lghighand!(state::SmartsParserState) -> Union{Pair,Nothing}

And <- Not ('&'? Not)*

The argument `func` is a parser function which has a parser state as an argument,
process tokens found in the given text, and returns nothing if no valid tokens were found.
"""
function lghighand!(state::SmartsParserState, func)
    fmls = Set{QueryFormula}()
    fml = lgnot!(state, func)
    fml === nothing && return
    while fml !== nothing
        push!(fmls, fml)
        if read(state) == '&'
            forward!(state)
        end
        fml = lgnot!(state, func)
    end
    @assert !isempty(fmls) "(lghighand!) invalid AND(&) operation"
    setdiff!(fmls, [QueryFormula(:any, true)])
    if isempty(fmls)
        return QueryFormula(:any, true)
    elseif length(fmls) == 1
        return collect(fmls)[1]
    else
        return QueryFormula(:and, fmls)
    end
end

"""
    lgnot!(state::SmartsParserState) -> Union{Pair,Nothing}

Not <- '!'? Element

The argument `func` is a parser function which has a parser state as an argument,
process tokens found in the given text, and returns nothing if no valid tokens were found.
"""
function lgnot!(state::SmartsParserState, func)
    if read(state) == '!'
        forward!(state)
        fml = func(state)
        @assert fml !== nothing "(lgnot!) invalid NOT(!) operation"
        return QueryFormula(:not, fml)
    end
    return func(state)  # can be Nothing if the parser get stop token
end


text:
density(q::Quadrature, fs::P) where {N <: Int, P <: AbstractArray{<:Real, N}} =
    sum(fs, dims = N)
density(q::Quadrature, f::P) where {P <: AbstractVector{<:Real}} = sum(f)

function velocity!(
    q::Quadrature,
    f::P,
    ρ::T,
    u::VT,
) where {T <: Real, VT <: AbstractVector{T}, P <: AbstractVector{T}}
    @inbounds for d in 1:dimension(q)
        u[d] = zero(T)
        for idx in 1:length(q.weights)
            u[d] += f[idx] * q.abscissae[d, idx]
        end
        u[d] /= ρ
    end
    return
end

function pressure(
    q::Quadrature,
    f::P,
    ρ::T,
    u::VT,
) where {T <: Real, VT <: AbstractVector{T}, P <: AbstractVector{T}}
    a_2 =
        sum(f[idx] * hermite(Val{2}, q.abscissae[:, idx], q) for idx in 1:length(q.weights))
    D = dimension(q)

    p = (tr(a_2) - ρ * (u[1]^2 + u[2]^2 - D)) / D
    return p
end

function momentum_flux(
    q::Quadrature,
    f::Population,
    ρ::T,
    u::VT,
) where {T <: Real, VT <: AbstractVector{T}, Population <: AbstractVector{T}}
    D = dimension(q)
    P = zeros(D, D)
    @inbounds for x_idx in 1:D, y_idx in 1:D
        P[x_idx, y_idx] = sum(
            # Pressure tensor
            f[f_idx] *
            (q.abscissae[x_idx, f_idx] - u[x_idx]) *
            (q.abscissae[y_idx, f_idx] - u[y_idx])

            # Stress tensor
            # f[f_idx] * (q.abscissae[x_idx, f_idx]) * (q.abscissae[y_idx, f_idx])
            for f_idx in 1:length(f)
        ) #- ρ * u[x_idx] * u[y_idx]
    end

    return q.speed_of_sound_squared * P #- I * pressure(q, f, ρ, u)

    E = 0.0
    @inbounds for idx in 1:length(f)
        E += f[idx] * (q.abscissae[1, idx]^2 + q.abscissae[2, idx]^2)
    end

    p = q.speed_of_sound_squared * (E - ρ * (u[1]^2 + u[2]^2)) / D
    return p
end

function temperature(
    q::Quadrature,
    f::P,
    ρ::T,
    u::VT,
) where {T <: Real, VT <: AbstractVector{T}, P <: AbstractVector{T}}
    return pressure(q, f, ρ, u) ./ ρ
end

"""
Computes the deviatoric tensor σ

τ is the relaxation time such that ν = cs^2 τ
"""
function deviatoric_tensor(
    q::Quadrature,
    τ::T,
    f::P,
    ρ::T,
    u::VT,
) where {T <: Real, VT <: AbstractVector{T}, P <: AbstractVector{T}}
    D = dimension(q)

    a_bar_2 = sum([
        f[idx] * hermite(Val{2}, q.abscissae[:, idx], q) for idx in 1:length(q.weights)
    ])
    a_eq_2 = equilibrium_coefficient(Val{2}, q, ρ, u, 1.0)
    σ = (a_bar_2 - a_eq_2) / (1 + 1 / (2 * τ))
    return σ - I * tr(σ) / D
end

"""
Compute the temperature from hermite coefficients
"""
function temperature(
    q::Q,
    f::VT,
    a_0::VT,
    a_1::VT,
    a_2::MT,
) where {Q <: Quadrature, T <: Real, VT <: AbstractVector{T}, MT <: AbstractMatrix{T}}
    D = dimension(q)
    P = Array{T}(undef, D, D)
    ρ = a_0
    u = a_1 / ρ
    P = q.speed_of_sound_squared * a_2 - ρ * (u * u' - I(2))

    return tr(P) / (D * ρ)
end


text:
using Documenter

makedocs()

open_in_new_tab( joinpath(HOMEDIR,"../docs/build/index.md") )


text:
md"""
# Элементы функционального программирования

## Функции - как объекты 1-го класса
"""
x = 0

s = sin  # <: Function

md"""
Пусть, например, мы хотим перемещать робота в заданном направлении до тех пор, пока не наступит некоторое событие,
выраженное соответствующим условием.

Это реализуется следующей функцией:
"""

function movements!(condition::Function, r, side)
    n=0
    while !condition(r)
        move!(r,side)
        n+=1
    end
    return n
end

#Теперь эту функцию можно вызывать, например, так:

movements!(isborder, r, side)

md"""
## Анонимные функции
Анонимные функции представляют собой функциональные значения (<: Function) без имени
 (подобно тому, как бывают значения типа Int без имени)

Например:

condition = robot -> !isborder(robot)

Здесь переменной condition присвоено некоторое функциональное значение (значение анонимной функции).

Другие римеры:

f = x -> 2x

f(2)
4

(x -> 2x)(3)
6

f2= (x,y) -> x^2+y^2

f2(2,3)
13

g= _ -> 1

g(1)
1

g(2)
1

g(3)
1

g()
ERROR

g(1,2)
ERROR

g2 = (_,_) -> 1

g2(1,2)
1

Анонимные функции часто также называю ЛЯМБДА- функциями
(это название пришло из соответствующего раздела математической логики, из так называемого лямбда-исчислниия)

В языке Python, например, определение фанкции

x -> 2x

выглядело бы так:

lambda x: 2*x

При этом, однако, в отличие от языка Python, в языке Julia нет тех сильных ограничений, которые есть в языке Python на вид выражения,
определяющее анонимную функцию.

## Функции высшего порядка

Функция называется  функцией высшего порядка, если она имеет хотя бы один аргумент функционального типа
или если она возвращает значение функционального типа.

Пример встроенной функции высшего порядка - map

julia> map(x->x^2, [1,2,3,4])
4-element Vector{Int64}:
  1
  4
  9
 16

julia> map(x->x^2, (1,2,3,4))
(1, 4, 9, 16)

julia> map(sin, (1,2,3,4))
(0.8414709848078965, 0.9092974268256817, 0.1411200080598672, -0.7568024953079282)

## do-синтаксис

julia> map([1,2,3,4,5]) do x
    x^2
end
5-element Vector{Int64}:
1
4
9
16
25

julia> map([1,2,3,4,5]) do _
    0
end
5-element Vector{Int64}:
0
0
0
0
0

do-синтаксис применим не только к встоеннной функции map, но и к любой другой, у котрой первый аргумент
имеет функциональный тип.

Этот синтаксис просто позволяет определить фактическое значение этого параметра в этой компактной конструкци,
без того, чтобы определять эту функцию отдельно стандартным образом.

## Пример функции высшего порядка, возвращающей функциональное значение. Понятие замыкания.

julia> function g(a)
           return x -> a*x # возвращается функциональное значение (функция)
       end
g (generic function with 1 method)

julia> h = g(2)
#33 (generic function with 1 method)

julia> h(2)
4

julia> h(3)
6

julia> h = g(3)
#33 (generic function with 1 method)

julia> h(3)
9

## Примеры практического использования функций высшего порядка.

Ранее были у на  такие функции

    movements!(robot, side) # перемещает робота "до упора" в заданном направлении и возвращает число сделанных "шагов"

    movements!(robot, side, num_steps) # делает num_steps шагов в заданном напралении

Их можно было бы дополнить ещё двумя методами - двумя функциями высшего порядка:

    movements!(condition::Function, robot, side)

и, соответственно,

    movements!(condition::Function, robot, side, num_steps)

где в в обоих случаях через функциональный параметр condition в функцию передаётся условие,
только и только при выполнении которого робот должен перемещаться.
При этом во втором случае через параметр num_steps передается еще наибольшее число шагов, которое робот может
сделать в заданном направлении (даже если условие, определяемое condition все время выполняется).

Функции высшего порядка movements! (новые методы функции movements!, в терминологии языка julia) могут быть
реализованы так, что функция condition будет иметь один аргумент (robot), или, в другом варианте реализации,
вообще не будут иметь аргументов.

### Варианты реализации функции movements!(condition::Function, robot, side)

1-вариант
"""
"""
movements!(condition::Function, robot, side)

-- возвращает число шагов, сделанных "до тех пор, пока выпоняется условие condition"

-- condition - функция одного аргумента (типа робот), возвращающая логическое значение
"""
function movements!(condition::Function, r, side)

    while condition(r) == true
        move!(r,side)
    end
end

# Примеры использования этой функции:

r=Ribot() # тип робота мог бы быть и другим

movements!(r->!isborder(r,Ost), r, Ost) # - это вызов функции

movements!(r, Ost) do robot
    robot->!isborder(robot, Ost)
end # - это вызов функции с использованием do-синтаксиса

movements!(ismarker, r, Ost)

movements!(robot -> ismarker(robot), r, Ost)

md"""
2-ой вариант реализации той же функции movements!(condition::Function, r, side)
(наверное, более предпочтительный)
"""

"""
movements!(condition::Function, robot, side)

-- возвращает число шагов, сделанных "до тех пор, пока выпоняется условие condition"

-- condition - функция без аргументов, возвращающая логическое значение
"""
function movements!(condition::Function, robot, side)
    #  число шагов, сделанных "до тех пор, пока выпоняется условие condition"
    while condition() == true
        move!(robot, side)
    end
end

# Пример использования так определенной функции movements!:

r=Ribot() # тип робота мог бы быть и другим

movements!(() -> !isborder(r, Ost), r, Ost) # - это вызов функции

movements!(r, Ost) do
    robot->!isborder(r, Ost)
end # - это вызов функции с использованием do-синтаксиса

md"""
### Реализация функции movements!(condition::Function, robot, side)

Определим эту функцию (обобщенную), так, чтобы она перемещала робота в заданном направлении не более чем на
num_steps шагов, пока выпоняется условие condition()==true". При этом она будет возвращать число сделанных роботом
"шагов" до своей остановки.
"""

"""
movements!(condition::Function, robot, side, num_steps)

-- делает не более чем num_steps шагов в заданном направлении, пока выпоняется условие condition()==true",
и возвращает число сделанных шагов

-- condition - функция без аргументов, возвращающая логическое значение
"""
function movements!(condition::Function, r, side, num_steps)
    n = 0
    for _ in 1:num_steps
        if condition() == true
            move!(r, side)
            n += 1
        else
            return n
        end
    end
    return n
end

md"""
Такая функция была бы удобной, например, для решения задачи, в которой требуется
на неограниченном поле нати маркер (была у нас такая задача).

Так же её можно было бы использовать для реализации функции try_move!(robot, side::HorizonSide).

Рассмотрим ещё примеры некоторых других полезных (обобщенных) функций высшего порядка, которыми стоило бы
поплнить нашу библиотеку функций.

1. Функция

    shuttle(condtion::Function, robot, side)

перемещающая робота взад-вперед (челноком, as shuttle) с увелисением амплитуды перемещений
до тех пор, пока не будет выполено заданное условие,
где side - начальное направление,

2. Функция, перемещающая робота змейкой

    snake(condition::Function, robot, side::NTuple{HorizonSide})

до тех пор, пока не будет выполено заданное условие condition()==true,
где side - пара направлений, определяющих траекторию змейки.

3. Функция, перемещающая робота по спирали

    spiral(condition::Function, robot, side)

до тех пор, пока не будет выполено заданное условие condition()==true,
где side - начальное направление.
"""


text:
using FCSFiles
using FileIO
using Test

flowrun = load("testdata/BD-FACS-Aria-II.fcs")

@test length(flowrun["SSC-A"]) == 100000


text:
 SEC1  = 1
 SEC2  = IS1
 SEC3  = IS2
 SEC4  = IS3
 SEC5  = IS4
 SEC6  = IS5
 SEC7  = IS6
 SEC8  = IS7
 SEC9  = IS8
 SEC10  = IS9
 SEC11  = IS10
 SEC12  = IS11
 SEC13  = IS12
 SEC14  = IS13
 SEC15  = IS14
 SEC16  = IS15
 SEC17  = IS1 *IS1
 SEC18  = IS1 *IS2
 SEC19  = IS2 *IS2
 SEC20  = IS16
 SEC21  = IS17
 SEC22  = IS18
 SEC23  = IS19
 SEC24  = IS20
 SEC25  = IS21
 SEC26  = IS22
 SEC27  = IS23
 SEC28  = IS24
 SEC29  = IS25
 SEC30  = IS26
 SEC31  = IS27
 SEC32  = IS1 *IS3
 SEC33  = IS2 *IS3
 SEC34  = IS1 *IS4
 SEC35  = IS2 *IS4
 SEC36  = IS1 *IS5
 SEC37  = IS2 *IS5
 SEC38  = IS1 *IS6
 SEC39  = IS2 *IS6
 SEC40  = IS1 *IS7
 SEC41  = IS2 *IS7
 SEC42  = IS28
 SEC43  = IS29
 SEC44  = IS30
 SEC45  = IS31
 SEC46  = IS32
 SEC47  = IS33
 SEC48  = IS34
 SEC49  = IS35
 SEC50  = IS36
 SEC51  = IS37
 SEC52  = IS38
 SEC53  = IS39
 SEC54  = IS40
 SEC55  = IS4 *IS4
 SEC56  = IS3 *IS5
 SEC57  = IS4 *IS5
 SEC58  = IS5 *IS5
 SEC59  = IS3 *IS6
 SEC60  = IS4 *IS6
 SEC61  = IS5 *IS6
 SEC62  = IS6 *IS6
 SEC63  = IS3 *IS7
 SEC64  = IS4 *IS7
 SEC65  = IS5 *IS7
 SEC66  = IS6 *IS7
 SEC67  = IS7 *IS7
 SEC68  = IS1 *IS8
 SEC69  = IS2 *IS8
 SEC70  = IS1 *IS9
 SEC71  = IS2 *IS9
 SEC72  = IS1 *IS10
 SEC73  = IS2 *IS10
 SEC74  = IS1 *IS11
 SEC75  = IS2 *IS11
 SEC76  = IS1 *IS12
 SEC77  = IS2 *IS12
 SEC78  = IS1 *IS13
 SEC79  = IS2 *IS13
 SEC80  = IS1 *IS14
 SEC81  = IS2 *IS14
 SEC82  = IS1 *IS15
 SEC83  = IS2 *IS15
 SEC84  = IS41
 SEC85  = IS42
 SEC86  = IS43
 SEC87  = IS44
 SEC88  = IS3 *IS12
 SEC89  = IS4 *IS12
 SEC90  = IS5 *IS12
 SEC91  = IS6 *IS12
 SEC92  = IS7 *IS12
 SEC93  = IS3 *IS13
 SEC94  = IS4 *IS13
 SEC95  = IS5 *IS13
 SEC96  = IS6 *IS13
 SEC97  = IS7 *IS13
 SEC98  = IS3 *IS14
 SEC99  = IS4 *IS14
 SEC100  = IS5 *IS14
 SEC101  = IS6 *IS14
 SEC102  = IS7 *IS14
 SEC103  = IS3 *IS15
 SEC104  = IS4 *IS15
 SEC105  = IS5 *IS15
 SEC106  = IS6 *IS15
 SEC107  = IS7 *IS15
 SEC108  = IS1 *IS16
 SEC109  = IS2 *IS16
 SEC110  = IS1 *IS17
 SEC111  = IS2 *IS17
 SEC112  = IS1 *IS18
 SEC113  = IS2 *IS18
 SEC114  = IS1 *IS19
 SEC115  = IS2 *IS19
 SEC116  = IS1 *IS20
 SEC117  = IS2 *IS20
 SEC118  = IS1 *IS21
 SEC119  = IS2 *IS21
 SEC120  = IS1 *IS22
 SEC121  = IS2 *IS22
 SEC122  = IS1 *IS23
 SEC123  = IS2 *IS23
 SEC124  = IS1 *IS24
 SEC125  = IS2 *IS24
 SEC126  = IS1 *IS25
 SEC127  = IS2 *IS25
 SEC128  = IS1 *IS26
 SEC129  = IS2 *IS26
 SEC130  = IS1 *IS27
 SEC131  = IS2 *IS27
 SEC132  = IS45
 SEC133  = IS46
 SEC134  = IS3 *IS21
 SEC135  = IS4 *IS21
 SEC136  = IS5 *IS21
 SEC137  = IS6 *IS21
 SEC138  = IS7 *IS21
 SEC139  = IS3 *IS22
 SEC140  = IS4 *IS22
 SEC141  = IS5 *IS22
 SEC142  = IS6 *IS22
 SEC143  = IS7 *IS22
 SEC144  = IS3 *IS23
 SEC145  = IS4 *IS23
 SEC146  = IS5 *IS23
 SEC147  = IS6 *IS23
 SEC148  = IS7 *IS23
 SEC149  = IS3 *IS24
 SEC150  = IS4 *IS24
 SEC151  = IS5 *IS24
 SEC152  = IS6 *IS24
 SEC153  = IS7 *IS24
 SEC154  = IS3 *IS25
 SEC155  = IS4 *IS25
 SEC156  = IS5 *IS25
 SEC157  = IS6 *IS25
 SEC158  = IS7 *IS25
 SEC159  = IS3 *IS26
 SEC160  = IS4 *IS26
 SEC161  = IS5 *IS26
 SEC162  = IS6 *IS26
 SEC163  = IS7 *IS26
 SEC164  = IS3 *IS27
 SEC165  = IS4 *IS27
 SEC166  = IS5 *IS27
 SEC167  = IS6 *IS27
 SEC168  = IS7 *IS27
 SEC169  = IS1 *IS28
 SEC170  = IS2 *IS28
 SEC171  = IS1 *IS29
 SEC172  = IS2 *IS29
 SEC173  = IS1 *IS30
 SEC174  = IS2 *IS30
 SEC175  = IS1 *IS31
 SEC176  = IS2 *IS31
 SEC177  = IS1 *IS32
 SEC178  = IS2 *IS32
 SEC179  = IS1 *IS33
 SEC180  = IS2 *IS33
 SEC181  = IS1 *IS34
 SEC182  = IS2 *IS34
 SEC183  = IS1 *IS35
 SEC184  = IS2 *IS35
 SEC185  = IS1 *IS36
 SEC186  = IS2 *IS36
 SEC187  = IS1 *IS37
 SEC188  = IS2 *IS37
 SEC189  = IS1 *IS38
 SEC190  = IS2 *IS38
 SEC191  = IS1 *IS39
 SEC192  = IS2 *IS39
 SEC193  = IS1 *IS40
 SEC194  = IS2 *IS40

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
module MonteCarloStudy

import Distributed
using Distributed: pmap

import PartiallyObservedInverseGames.InversePreSolve
using PartiallyObservedInverseGames.ForwardGame: solve_game
using PartiallyObservedInverseGames.InverseGames:
    InverseKKTConstraintSolver, InverseKKTResidualSolver, solve_inverse_game

import Random
import Statistics
import LinearAlgebra
import Distances
import VegaLite

import PartiallyObservedInverseGames.CostUtils
import PartiallyObservedInverseGames.DynamicsModelInterface
using PartiallyObservedInverseGames.TrajectoryVisualization: visualize_trajectory, visualize_trajectory_batch
using ProgressMeter: @showprogress
using VegaLite: @vlplot

#======================================== Generate Dataset =========================================#

include("dataset_generation.jl")
export generate_dataset

include("estimation.jl")
export estimate

include("forward_augmentation.jl")
export augment_with_forward_solution

include("estimator_statistics.jl")
export estimator_statistics

include("visualization.jl")
export visualize_paramerr, visualize_poserr
end


text:
using StaticArrays, Distributions, ProgressMeter, DataFrames

export neutral_growth_clones, neutral_growth_clones!

function neutral_growth_clones!(tumor::Vector{Haplotype}, clonesizes::Vector{Int}, obs::Vector{Vector{Int}}, endsize::Int, T=Inf; b, d, μ, t=0., Nthresh=endsize, Tthresh = T, showprogress=true )

    mutID = length(tumor)-1
    F = cumsum(getfield.(tumor,:n))
    N = last(F)

    if showprogress prog = ProgressUnknown("Progress: N ") end
    while 0 < N < endsize && t < T
        μp = N < Nthresh && t < Tthresh ? μ : 0.

        p = rand()*(b+d)
        m = rand(1:N)
        row = findfirst(i-> i >= m, F)
        parent = tumor[row]

        t += rand(Exponential(1.))/((b+d)*N)
        if p > d

            mnew =  (rand()<μp/2) + (rand()<μp/2)
            for n=1:mnew
                mutID += 1
                push!(tumor, Haplotype(push(parent.mutations, mutID), 1, t, row) )
                push!(F, N+n)

                push!(clonesizes, 1)

                push!(obs, clonesizes[parent.mutations])
            end
            if mnew == 0
                parent.n += 1
                F[row:end] .+= 1
            elseif mnew == 2
                parent.n -= 1
                F[row:end] .-= 1
            end
            clonesizes[parent.mutations] .+= 1
            N += 1

        else
            parent.n -= 1
            F[row:end] .-= 1
            N -= 1

            clonesizes[parent.mutations] .-= 1
        end

        if showprogress ProgressMeter.update!(prog, N) end
    end

    return (tumor=tumor, clonesizes=clonesizes, obs=obs, N=N, t=t)
end

function neutral_growth_clones(endsize, T=Inf; args... )
    mutID = 0
    tumor = [Haplotype(SVector{0, Int}(), 1, 0., 0 )]
    obs = Vector{Vector{Int}}()
    clonesizes = Int[]

    out = neutral_growth_clones!(tumor, clonesizes, obs, endsize, T; args... )
    return iszero(out[:N]) ? neutral_growth_clones(endsize, T; args... ) : out
end

# out = neutral_growth(100; b=1., d=0.2, μ=0.3)
#
# out.clonesizes
#
# tumor = out.tumor |> DataFrame

max_stars_repo_name: aangaji/Turnover

max_issues_repo_name: aangaji/Turnover

max_forks_repo_name: aangaji/Turnover


text:
if !haskey(ENV, "BUILDKITE")
    import Pkg
    Pkg.develop(Pkg.PackageSpec(; path = dirname(@__DIR__)))
end
import Thermodynamics
import RootSolvers
import CLIMAParameters
import ReportMetrics
const RM = ReportMetrics
td_dir = RM.mod_dir(Thermodynamics)

dirs_to_monitor = [
    ".",
    RM.mod_dir(Thermodynamics),
    RM.mod_dir(RootSolvers),
    RM.mod_dir(CLIMAParameters),
]

for constructor in ["ρeq", "pθq", "pTq"]
    ENV["ALLOCATION_CONSTRUCTOR"] = constructor
    RM.report_allocs(;
        job_name = constructor,
        run_cmd = `$(Base.julia_cmd()) --project=perf/ --track-allocation=all perf/alloc_per_constructor.jl`,
        dirs_to_monitor = dirs_to_monitor,
        process_filename = x -> last(split(x, "packages/")),
    )
end


text:
using HDF5

function generate_files(path,N)
    hdf5_file = path*"cumsum.hdf5"
    csv_file = path*"cumsum.csv"
    if isfile(hdf5_file)
        rm(hdf5_file)
    end
    if isfile(csv_file)
        rm(csv_file)
    end
    # uniform ids in a range for load balance among nodes
    id = collect(1:N)
    x = rand(N)
    y = rand(N)
    h5write(hdf5_file,"/id", id)
    h5write(hdf5_file,"/x", x)
    h5write(hdf5_file,"/y", y)
    writecsv(csv_file, zip(id,x,y))
end

#N = 256*10^6
N = 256
path = ENV["HOME"]*"/tmp/"
generate_files(path, N)

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
include("plotter_setup.jl")
using LsqFit
data_folder_name = "nonmoralizing_qsw"

data_goe = log10.(npzread("$data_folder_name/goe_loc_ham.npz")[:,1])
data_gue = log10.(npzread("$data_folder_name/gue_loc_ham.npz")[:,1])
n = length(data_goe)
@assert length(data_gue) == n

cla()
figure(figsize=[4, 2.4])

hist(data_goe, label="GOE", color="g", bins=70, alpha=.4, ec="k", histtype="stepfilled")
hist(data_gue, label="GUE", color="r", bins=70, alpha=.7, ec="k", histtype="stepfilled")
ylabel("frequency")

ylims = ylim()
vlines(minimum(data_goe), ylims..., color="b", linestyle="-", linewidth=.5)
vlines(minimum(data_gue), ylims..., color="r", linestyle="--", linewidth=.5)
ylim(ylims...)
xlabel(L"\log_{10} (|\lambda|)")
leg = legend(loc=2)

savefig("plots/nonmoralizing_qsw/loc_ham_analysis.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))

##

data_nonsym = npzread("$data_folder_name/probdist_non_symmetric.npz")
data_sym = npzread("$data_folder_name/probdist_symmetric.npz")
n = length(data_sym)
#

positions = -(div(n, 2)):(div(n, 2))
markersize = 3
mew = .5
linewidth = .3

cla()
figure(figsize=[2.8, 2])

plot(positions, data_nonsym, "-+r", label="NGQSW", linewidth=linewidth, markersize=markersize+1, mew=mew)
plot(positions, reverse(data_nonsym), "-xk", label="NGQSW rev.", linewidth=linewidth, markersize=markersize, mew=mew)
vlines(0, 0, .06, linestyles="--", colors="k",linewidth=.5)

legend(loc=1, fontsize=6)
ylim(-.005, 0.065)
xlabel("position")
ylabel("measurement probability")
savefig("plots/nonmoralizing_qsw/nonsymmetric.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))

cla()
figure(figsize=[2.8, 2])

plot(positions, data_sym, "-+r", label="NGQSW", linewidth=linewidth, markersize=markersize+1, mew=mew)
plot(positions, reverse(data_sym), "-xk", label="NGQSW rev.", linewidth=linewidth, markersize=markersize, mew=mew)
vlines(0, 0, .06, linestyles="--", colors="k",linewidth=.5)

legend(loc=1, fontsize=6)
ylim(-.005, 0.065)
yticks(visible=false)
xlabel("position")

savefig("plots/nonmoralizing_qsw/symmetric.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))

## scaling new
ω_select = [1,2,4,6]
ωs = (0.5:0.1:1.)[ω_select]
max_n = 300
t_scales = 0.1:0.1:6
t_scales_cont = 0.1:0.1:6
batch = 5

data = npzread("$data_folder_name/scaling_data.npz")[ω_select,:]

data_t = collect(t_scales .* max_n)[1:size(data,2)]
data_t_cont = collect(t_scales_cont .* max_n)
println(data_t[[1,end]])

cla()
figure(figsize=[5, 2.4])
b_half = div(5, 2)

#linetyles = ["--", ":", "-.", "-"]
style = ["+", ".", "d", "x"]
marksizes = [5, 4.5, 3, 4]

for (i, ω) = enumerate(ωs)
    _, data_y = scaling_exponent_generator(log.(data_t), data[i,:], batch)

    plot(data_t[b_half+1:end-b_half], data_y, "$(style[i])", label=latexstring("\\omega = $ω"), markersize=marksizes[i])
    if i < 4
        @. model(x, p) = p[1] - p[2]/(x-p[3])^p[4]
        lb = [1.5, 10. , -Inf, 0.]
        ub = [2.5, Inf, 20, 1]
        p0 = [2., 15, 1, 1]
        fit = curve_fit(model, data_t[b_half+20:end-b_half], data_y[20:end], p0, lower=lb, upper=ub)
        @show length(data_y[20:end])
        c1, c2, c3, α = fit.param
        @show ω, fit.param
        @show stderror(fit)

        plot(data_t_cont, c1 .- c2 ./ (data_t_cont .- c3) .^ α, "--k", linewidth=.8)
        println()
    end
end

legend(loc=4, fontsize=7, bbox_to_anchor=(1., 0.15))
hlines([1,2], 0, maximum(data_t_cont), linestyles="-", linewidth=.6)

xlabel("time \$t\$")
ylabel("scaling exponent \$\\alpha_t \$")
savefig("plots/nonmoralizing_qsw/scaling.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))


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

using Markdown
using InteractiveUtils

# ╔═╡ 58d05c3e-8774-11eb-2f17-3daba3edae13
begin
  using Pkg; Pkg.activate("MLJ_env", shared=true)
  using PlutoUI

    using CategoricalArrays
    import MLJLinearModels
    import DataFrames
    import CSV
    import DecisionTree
    using MLJ
    import MLJLinearModels
    import MultivariateStats
    import MLJFlux
    using Plots
    using Test
end

# ╔═╡ c3f6ddc0-8782-11eb-0197-25fdd80a50fd
begin
  using Tables

  data = (a = [1, 2, 3, 4],
    b = rand(4),
    c = rand(4),
    d = coerce(["male", "female", "female", "male"], OrderedFactor)
  )
  pretty(data)

  y₂, X₂, w₂ = unpack(data, ==(:a),
    name -> elscitype(Tables.getcolumn(data, name)) == Continuous,
    name -> true);

end

# ╔═╡ 2bb6648a-8774-11eb-17b7-0f6fef114c61
md"""
## A workshop introducing the machine learning toolbox MLJ - Part2 - Selecting, Training and Evaluating Models

  - source:[MLJ](https://alan-turing-institute.github.io/MLJ.jl/stable/)
  - based on [Machine Learning in Julia using MLJ, JuliaCon2020](https://github.com/ablaom/MachineLearningInJulia2020)
"""

# ╔═╡ 7240757e-8774-11eb-3c51-eb10a4c46db3
md"""
##### General resources

  - [List of methods introduced in this tutorial](methods.md)
  - [MLJ Cheatsheet](https://alan-turing-institute.github.io/MLJ.jl/dev/mlj_cheatsheet/)
  - [Common MLJ Workflows](https://alan-turing-institute.github.io/MLJ.jl/dev/common_mlj_workflows/)
  - [MLJ manual](https://alan-turing-institute.github.io/MLJ.jl/dev/)
  - [Data Science Tutorials in Julia](https://alan-turing-institute.github.io/DataScienceTutorials.jl/)
"""

# ╔═╡ 992cc480-8785-11eb-0572-c33c1c988efc
html"""

"""

# ╔═╡ 72108d16-8774-11eb-0003-49693e396f19
md"""
##### TOC
  - [Selecting, Training and Evaluating Models](#selecting-training-and-evaluating-models)
  - [Step 1](#step1)
  - [Step 2](#step2)
  - [Step 3](#step3)
  - [Step 4](#step4)
  - [On learning curves](#learning_curves)
  - [Resources or this part](#resources)
  - [Exercises for Part 2](#exercises)
"""

# ╔═╡ 71d64138-8774-11eb-00ff-a3d20a55ccc8
md"""
 **Goals:**
 1. Search MLJ's database of model metadata to identify model candidates for a supervised learning task.
 2. Evaluate the performance of a model on a holdout set using basic `fit!`/`predict` work-flow.
 3. Inspect the outcomes of training and save these to a file.
 3. Evaluate performance using other resampling strategies, such as cross-validation, in one line, using `evaluate!`
 4. Plot a "learning curve", to inspect performance as a function of some model hyper-parameter, such as an iteration parameter

The "Hello World!" of machine learning is to classify Fisher's
famous iris data set. This time, we'll grab the data from
[OpenML](https://www.openml.org):
"""

# ╔═╡ e1fa26aa-8774-11eb-204a-a3c824321016
begin
  iris = OpenML.load(61); # a row table
  iris = DataFrames.DataFrame(iris);
  first(iris, 4)
end

# ╔═╡ 4f4601e2-8785-11eb-2401-e5a978807dc6
html"""

  back to TOC

"""

# ╔═╡ e1e48b42-8774-11eb-0e4e-07edf1692567
md"""
**Main goal.** To build and evaluate models for predicting the  `:class` variable, given the four remaining measurement variables.

##### Step 1. Inspect and fix scientific types
"""

# ╔═╡ e1c675c6-8774-11eb-0cc8-9fcd3ec74dff
schema(iris)

# ╔═╡ e1ad7eae-8774-11eb-0040-57d9e0fb1236
begin
  coerce!(iris, :class => Multiclass);
  schema(iris)
end

# ╔═╡ 720cb48c-8785-11eb-0649-1f4163b49de6
html"""

  back to TOC

"""

# ╔═╡ e1916f20-8774-11eb-3503-ef998f66e2bd
md"""
##### Step 2. Split data into input and target parts

Here's how we split the data into target and input features, which is needed for MLJ supervised models. We randomize the data at the same time:
"""

# ╔═╡ 2f9e782a-8775-11eb-1272-efb3abb86766
begin
  yᵢ, Xᵢ = unpack(iris, ==(:class), name->true; rng=123);
  scitype(yᵢ)
end

# ╔═╡ 2f6944fa-8775-11eb-1111-df2b2b267f84
md"""
###### On searching for a model

Here's how to see *all* models (not immediately useful):
"""

# ╔═╡ 2f50d958-8775-11eb-096e-0b8998321174
all_models = models()

# ╔═╡ 2f37d368-8775-11eb-2f0c-fd13ca4fe29d
meta = all_models[3]

# ╔═╡ 2f1ec62a-8775-11eb-26e2-85c032960d6a
meta_xgboost = all_models[153]

# ╔═╡ 2f0eb12c-8775-11eb-17f9-235038e0b0a7
targetscitype = meta.target_scitype

# ╔═╡ 2ef084b8-8775-11eb-2bc0-0b949d966e90
scitype(yᵢ) <: targetscitype
# So this model won't do. Let's  find all pure julia classifiers...

# ╔═╡ 2ed935ba-8775-11eb-08e9-d74a578c22fa
begin
  filter_julia_classifiers(meta) = AbstractVector{Finite} <: meta.target_scitype &&  meta.is_pure_julia

  models(filter_julia_classifiers)
end

# ╔═╡ 2ec20a02-8775-11eb-2449-5749cbbb60e1
md"""
Find all models with "Classifier" in `name` (or `docstring`):
"""

# ╔═╡ 2ea8ed54-8775-11eb-18d6-d1d30c96d908
models("Classifier")

# ╔═╡ e1782812-8774-11eb-1b1c-3faefea1001c
md"""
Find all (supervised) models that match my data!
"""

# ╔═╡ e109c3ae-8774-11eb-3f53-1dc77c89fc61
models(matching(Xᵢ, yᵢ))

# ╔═╡ 77d08a42-8785-11eb-0e54-fda94fc25cbe
html"""

  back to TOC

"""

# ╔═╡ 7e2f4de8-8775-11eb-3cc3-510cbf8258b6
md"""
##### Step 3. Select and instantiate a model

To load the code defining a new model type we use the `@load` macro, which returns an *instance* of the type, with default hyperparameters:
"""

# ╔═╡ 7e17a8fa-8775-11eb-267a-f9732276336b
begin
  _model = @load NeuralNetworkClassifier pkg=MLJFlux
  modelᵢ = _model(epochs=15)
  info(modelᵢ)
end

# ╔═╡ 7df43b38-8775-11eb-1160-ab24fe881df7
md"""
In MLJ a *model* is just a struct containing hyper-parameters, and that's all. A model does not store *learned* parameters. Models are mutable:
"""

# ╔═╡ 7dc2a634-8775-11eb-1a2a-2b64699a3175
@test MLJFlux.NeuralNetworkClassifier(epochs=15) == modelᵢ

# ╔═╡ b14a0508-8775-11eb-2ea5-09289a0926b7
md"""
In MLJ a model and training/validation data are typically bound together in a *machine* as follows:
"""

# ╔═╡ 7da2d596-8775-11eb-22bd-bf94ffdff871
machᵢ = machine(modelᵢ, Xᵢ, yᵢ);

# ╔═╡ 7d8e181a-8775-11eb-1c85-8d0179791c4d
md"""
A machine stores *learned* parameters, among other things. We'll train this machine on 70% of the data and evaluate on a 30% holdout set. Let's start by dividing all row indices into `train` and `test` subsets:
"""

# ╔═╡ 7d70204e-8775-11eb-1bd3-111d1dd44c45
begin
  trainᵢ, testᵢ = partition(eachindex(yᵢ), 0.7)
  size(trainᵢ), size(testᵢ)
end

# ╔═╡ 7d5a596c-8775-11eb-1059-c38021f064e2
md"""
And now we can `fit` the model...
"""

# ╔═╡ 7d3c5bce-8775-11eb-242f-af3448ec4f6c
fit!(machᵢ, rows=trainᵢ, verbosity=2);

# [ Info: Loss is 1.098
# [ Info: Loss is 1.095
# [ Info: Loss is 1.091
# [ Info: Loss is 1.085
# [ Info: Loss is 1.079
# [ Info: Loss is 1.072
# [ Info: Loss is 1.067
# [ Info: Loss is 1.059
# [ Info: Loss is 1.052
# [ Info: Loss is 1.043
# [ Info: Loss is 1.033
# [ Info: Loss is 1.026
# [ Info: Loss is 1.015
# [ Info: Loss is 1.008
# [ Info: Loss is 0.9971

# ╔═╡ e0397248-8775-11eb-31f1-4745600b4927
md"""
... and `predict`:
"""

# ╔═╡ e020dba0-8775-11eb-04c6-5791fa04a0d0
predict(machᵢ, rows=testᵢ)  # or `predict(machᵢ, Xnew)`

# ╔═╡ dfebbeea-8775-11eb-0e03-19641479d64d
md"""
After training, one can inspect the learned parameters:
"""

# ╔═╡ dfd1df34-8775-11eb-083d-1f134983ca77
fitted_params(machᵢ)

# ╔═╡ f74b80a2-8775-11eb-3733-7f27c7d20f3d
md"""
Everything else the user might be interested in is accessed from the  training *report*:
"""

# ╔═╡ f7265138-8775-11eb-2873-6b37109e57df
report(machᵢ)

# ╔═╡ f70c5bca-8775-11eb-07c1-dfff7ad49f4a
md"""
You save a machine like this:
"""

# ╔═╡ 09a19f7a-8776-11eb-1203-f77b73622b87
MLJ.save("neural_net.jlso", machᵢ)

# ╔═╡ f6ef9698-8775-11eb-1400-7bc68f3dbe0f
md"""
And retrieve it like this:
"""

# ╔═╡ 15b16a2a-8776-11eb-1592-3dce20d4a188
begin
  mach₂ = machine("neural_net.jlso")
  predict(mach₂, Xᵢ)[1:3]
end

# ╔═╡ 15935300-8776-11eb-0c72-0f19462b9b10
md"""
If we want to fit a retrieved model, we will need to bind some data to it:
"""

# ╔═╡ 1576a3fe-8776-11eb-2e19-f7266559cf3f
begin
  mach₃ = machine("neural_net.jlso", Xᵢ, yᵢ)
  fit!(mach₃);
end

# ╔═╡ 2f56cae2-8776-11eb-1244-1f49372a694b
md"""
Machines remember the last set of hyper-parameters used during fit, which, in the case of iterative models, allows for a warm restart of computations in the case that only the iteration parameter is increased:
"""

# ╔═╡ 2f3979ec-8776-11eb-023e-c54951130c1a
begin
  modelᵢ.epochs = modelᵢ.epochs + 3
  fit!(machᵢ, rows=trainᵢ, verbosity=2);
end

# [ Info: Loss is 0.9844
# [ Info: Loss is 0.9706
# [ Info: Loss is 0.9555

# ╔═╡ 2f1fe682-8776-11eb-163d-67fa9dec2081
md"""
By default (for this particular model) we can also increase `:learning_rate` without triggering a cold restart:
"""

# ╔═╡ dfb95914-8775-11eb-0da8-bf42ce7e7350
begin
  modelᵢ.epochs = modelᵢ.epochs + 4
  modelᵢ.optimiser.eta = 10 * modelᵢ.optimiser.eta
  fit!(machᵢ, rows=trainᵢ, verbosity=2);
end

# [ Info: Loss is 0.7822
# [ Info: Loss is 0.7047
# [ Info: Loss is 0.6378
# [ Info: Loss is 0.5993

# ╔═╡ dfa1167e-8775-11eb-2245-a3e41d6512ac
md"""
However, change any other parameter and training will restart from scratch:
"""

# ╔═╡ 6af7b84a-8776-11eb-260a-8f04dee6e55a
begin
  modelᵢ.lambda = 0.001
  fit!(machᵢ, rows=trainᵢ, verbosity=2);
end

# [ Info: Loss is 1.055
# [ Info: Loss is 0.9297
# [ Info: Loss is 0.7862
# [ Info: Loss is 0.7138
# [ Info: Loss is 0.6547
# [ Info: Loss is 0.6223
# [ Info: Loss is 0.5939
# [ Info: Loss is 0.5621
# [ Info: Loss is 0.5432
# [ Info: Loss is 0.5187
# [ Info: Loss is 0.5065
# [ Info: Loss is 0.4927
# [ Info: Loss is 0.4777
# [ Info: Loss is 0.4196
# [ Info: Loss is 0.4481
# [ Info: Loss is 0.4189
# [ Info: Loss is 0.4034
# [ Info: Loss is 0.4599
# [ Info: Loss is 0.3912
# [ Info: Loss is 0.4244
# [ Info: Loss is 0.412
# [ Info: Loss is 0.4474
# [ Info: Loss is 0.4246
# [ Info: Loss is 0.389
# [ Info: Loss is 0.3978

# ╔═╡ 6adb6884-8776-11eb-12b2-8dc0829f1bc7
md"""
Let's train silently for a total of 50 epochs, and look at a prediction:
"""

# ╔═╡ 6aa68b5a-8776-11eb-2cbc-258c26de524c
begin
  modelᵢ.epochs = 50
  fit!(machᵢ, rows=trainᵢ);
  ŷᵢ = predict(machᵢ, Xᵢ[testᵢ, :]); # or predict(mach, rows=test)
  ŷᵢ[1]
end
# Optimising neural net:100%[=========================] Time: 0:00:00

# ╔═╡ df8118d8-8775-11eb-234c-a35cf41d7430
md"""
What's going on here?
"""

# ╔═╡ 989e2202-8776-11eb-03a1-49fa4508f71a
info(modelᵢ).prediction_type

# ╔═╡ 98800d12-8776-11eb-1c7c-457f310a3a4d
md"""
**Important**:
  - In MLJ, a model that can predict probabilities (and not just point values) will do so by default. (These models have supertype `Probabilistic`, while point-estimate predictors have supertype `Deterministic`.)
  - For most probabilistic predictors, the predicted object is a `Distributions.Distribution` object, supporting the `Distributions.jl` [API](https://juliastats.org/Distributions.jl/latest/extends/#Create-a-Distribution-1) for such objects. In particular, the methods `rand`,  `pdf`, `mode`, `median` and `mean` will apply, where appropriate.

So, to obtain the probability of "Iris-virginica" in the first test prediction, we do
"""

# ╔═╡ 9868dd68-8776-11eb-3633-afcbd37b82b6
pdf(ŷᵢ[1], "Iris-virginica")

# ╔═╡ 9850ae82-8776-11eb-33bd-4bfbb7eac41b
md"""
To get the most likely observation, we do
"""

# ╔═╡ 98393bc6-8776-11eb-0659-f7225d5839f0
mode(ŷᵢ[1])

# ╔═╡ 981f0ce2-8776-11eb-1546-83d6d1ec5727
md"""
These can be broadcast over multiple predictions in the usual way:
"""

# ╔═╡ be74448e-8776-11eb-3018-c722a9182899
broadcast(pdf, ŷᵢ[1:4], "Iris-versicolor")

# ╔═╡ be58b4d0-8776-11eb-3eea-3350fe77c4f2
mode.(ŷᵢ[1:4])

# ╔═╡ be3edd94-8776-11eb-31c6-634f25bacbc2
md"""
Or, alternatively, you can use the `predict_mode` operation instead of `predict`:
"""

# ╔═╡ be24b926-8776-11eb-09c2-335929294063
predict_mode(machᵢ, Xᵢ[testᵢ,:])[1:4] # or predict_mode(mach, rows=test)[1:4]

# ╔═╡ be14985e-8776-11eb-341f-f9c6c418f697
md"""
For a more conventional matrix of probabilities you can do this:
"""

# ╔═╡ bdf3586a-8776-11eb-0ac3-ddc972c36518
begin
  L = levels(yᵢ)
  pdf(ŷᵢ, L)[1:4, :]
end

# ╔═╡ bddcc460-8776-11eb-32a8-45a1bda54031
md"""
However, in a typical MLJ work-flow, this is not as useful as you might imagine. In particular, all probabilistic performance measures in MLJ expect distribution objects in their first slot:
"""

# ╔═╡ e700cc38-8776-11eb-1c82-35ed07df74af
cross_entropy(ŷᵢ, yᵢ[testᵢ]) |> mean

# ╔═╡ e6e1d9c2-8776-11eb-243c-4db077b59c05
md"""
To apply a deterministic measure, we first need to obtain point-estimates:
"""

# ╔═╡ e6c6f486-8776-11eb-20ed-4140d2fdfc59
misclassification_rate(mode.(ŷᵢ), yᵢ[testᵢ])

# ╔═╡ e6ad53c6-8776-11eb-37e9-f521bde89017
md"""
We note in passing that there is also a search tool for measures analogous to `models`:
"""

# ╔═╡ 98060a3a-8776-11eb-1385-ffc7eac03387
measures()

# ╔═╡ 03a66230-8777-11eb-0a4c-bda208768889
# measures(matching(y)) # experimental!

# ╔═╡ 861ed130-8785-11eb-207d-1b60a935f028
html"""

  back to TOC

"""

# ╔═╡ 0388dd00-8777-11eb-11b6-23a5968841f7
md"""
##### Step 4. Evaluate the model performance

Naturally, MLJ provides boilerplate code for carrying out a model evaluation with a lot less fuss. Let's repeat the performance evaluation above and add an extra measure, `brier_score`:

"""

# ╔═╡ 1544115e-8777-11eb-28df-3167da05d7cb
evaluate!(machᵢ,
  resampling=Holdout(fraction_train=0.7),
  measures=[cross_entropy, brier_score])

# ╔═╡ 0376892a-8777-11eb-10fe-ebd9a30329b5
md"""
Or applying cross-validation instead:
"""

# ╔═╡ 03503764-8777-11eb-0793-598490c04eba
evaluate!(machᵢ,
  resampling=CV(nfolds=6),
  measures=[cross_entropy, brier_score])

# Evaluating over 6 folds: 100%[=========================] Time: 0:00:12

# ╔═╡ da97b186-8781-11eb-0591-757d7c34d247
md"""
Or, Monte Carlo cross-validation (cross-validation repeated randomized folds)
"""

# ╔═╡ da7f0320-8781-11eb-3047-832718ddd4d0
e = evaluate!(machᵢ, resampling=CV(nfolds=6, rng=123),
  repeats=3,
  measures=[cross_entropy, brier_score])

# Evaluating over 18 folds: 100%[=========================] Time: 0:00:35

# ╔═╡ da634a7c-8781-11eb-34ce-21239b737208
md"""
One can access the following properties of the output `e` of an evaluation: `measure`, `measurement`, `per_fold` (measurement for each fold) and `per_observation` (measurement per observation, if reported).
"""

# ╔═╡ da499c0a-8781-11eb-1ecc-c78a2894425e
e.measure

# ╔═╡ da312b96-8781-11eb-1a04-e709749917af
e.per_fold

# ╔═╡ da19bd4e-8781-11eb-2316-a77ff448d985
md"""
Finally note that we can restrict the rows of observations from which train and test folds are drawn, by specifying `rows=...`.

For  example, imagining the last 30% of target observations are `missing` you might have a work-flow like this:
"""

# ╔═╡ 03fb3db0-8782-11eb-22a3-abdefb466c7b
begin
  trainₓ, testₓ = partition(eachindex(yᵢ), 0.7)
  machₓ = machine(modelᵢ, Xᵢ, yᵢ)
  evaluate!(machₓ, resampling=CV(nfolds=6),
    measures=[cross_entropy, brier_score],
    rows=trainₓ)            # cv estimate, resampling from `train`
  fit!(machₓ, rows=trainₓ)    # re-train using all of `train` observations
  predict(machₓ, rows=testₓ); # and predict missing targets
end

# ╔═╡ f980ec9a-8784-11eb-3c1f-cbd12e90eb12
html"""

"""

# ╔═╡ 03e5ccf8-8782-11eb-3b06-3fd79d29ece7
md"""
##### On learning curves

Since our model is an iterative one, we might want to inspect the out-of-sample performance as a function of the iteration parameter. For this we can use the `learning_curve` function (which, incidentally can be applied to any model hyper-parameter). This  starts by defining a one-dimensional range object for the parameter  (more on this when we discuss tuning in Part 4):
"""

# ╔═╡ 03ca941a-8782-11eb-3d2c-992431e1eaf2
begin
  r = range(modelᵢ, :epochs, lower=1, upper=50, scale=:log)

  curve = learning_curve(machₓ, range=r,
    resampling=Holdout(fraction_train=0.7), # (default)
    measure=cross_entropy)
end

# ╔═╡ 03b234ea-8782-11eb-259a-0566fa974775
plot(curve.parameter_values, curve.measurements, lw=3,
  xlabel="epochs",
  ylabel="cross entropy on holdout set")

# ╔═╡ 70bb451e-8774-11eb-18df-512352f9b934
html"""

  back to TOC

"""

# ╔═╡ 704ad2e8-8774-11eb-2083-812bdb1f4ffe
md"""
##### Resources for Part 2

  - From the MLJ manual:
    - [Getting Started](https://alan-turing-institute.github.io/MLJ.jl/dev/getting_started/)
    - [Model Search](https://alan-turing-institute.github.io/MLJ.jl/dev/model_search/)
    - [Evaluating Performance](https://alan-turing-institute.github.io/MLJ.jl/dev/evaluating_model_performance/) (using `evaluate!`)
    - [Learning Curves](https://alan-turing-institute.github.io/MLJ.jl/dev/learning_curves/)
    - [Performance Measures](https://alan-turing-institute.github.io/MLJ.jl/dev/performance_measures/) (loss functions, scores, etc)
- From Data Science Tutorials:
    - [Choosing and evaluating a model](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/choosing-a-model/)
    - [Fit, predict, transform](https://alan-turing-institute.github.io/DataScienceTutorials.jl/getting-started/fit-and-predict/)

"""

# ╔═╡ fbd65648-8774-11eb-1556-135f1f66e11a
html"""

  back to TOC

"""

# ╔═╡ b8924b1c-8774-11eb-12f0-75a0c9247554
md"""
#### Exercises for Part 2

##### Exercise 1

(a) Identify all supervised MLJ models that can be applied (without type coercion or one-hot encoding) to a supervised learning problem with input features `X4` and target `y4` defined below:
"""

# ╔═╡ c2ca6d30-8774-11eb-2f9f-853ae8843263
begin
  import Distributions
  poisson = Distributions.Poisson

  age = 18 .+ 60 * rand(10);
  # do not use symbol => crash!
  salary = coerce(rand(["big", "huge", "small"], 10), OrderedFactor);
  levels!(salary, ["small", "big", "huge"]);

  X₁ = DataFrames.DataFrame(age=age, salary=salary)
end

# ╔═╡ 87ccbfea-8782-11eb-28c5-1fd0c3ca976e
begin
  n_devices(salary) = salary > "small" ? rand(poisson(1.3)) : rand(poisson(2.9))
  y₁ = [n_devices(row.salary) for row in eachrow(X₁)];
end

# ╔═╡ 877b6546-8782-11eb-0787-7f71094de4b0
schema(X₁)

# ╔═╡ 869d4842-8782-11eb-0863-77d366291f93
models(matching(X₁, y₁))

# ╔═╡ a810df34-8782-11eb-1729-19864c77e6f0
md"""
##### Exercise 2 (unpack)

After evaluating the following ...
"""

# ╔═╡ c2a292b0-8774-11eb-11d7-b5a05aa9f616
md"""
...attempt to guess the evaluations of the following:
"""

# ╔═╡ c27462d2-8774-11eb-102e-135655b15bc4
begin
  @test y₂ == data.a
  @test typeof(y₂) <: AbstractVector
end

# ╔═╡ d824285c-8782-11eb-09be-7bfa9e2f63dc
# Tables.columnnames(data)
typeof(X₂) # Tables.getcolumn(data, [:a, :b])

# ╔═╡ d809473a-8782-11eb-2331-19eccaff712a
begin
  # Xₙ contains only the continuous attributres from data
  @test X₂ == (b=data.b, c=data.c)
end

# ╔═╡ c2563758-8774-11eb-0903-dd227f7618d1
with_terminal() do
  pretty(X₂)
end

# ╔═╡ 105a8fe8-8783-11eb-2f0b-8fb550d77372
begin
  @test w₂ == data.d
  @test typeof(w₂) <: AbstractVector
end

# ╔═╡ 103d8df8-8783-11eb-2c3d-eff502ac833c
md"""
##### Exercise 3 (first steps in modeling Horse Colic)

(a) Suppose we want to use predict the `:outcome` variable in the
Horse Colic study introduced in Part 1, based on the remaining
variables that are `Continuous` (one-hot encoding categorical
variables is discussed later in Part 3) *while ignoring the others*.

Extract from the `horse` data set (defined in Part 1) appropriate
input features `X` and target variable `y`. (Do not, however,
randomize the observations.)
"""

# ╔═╡ 3de19b84-8783-11eb-139d-9d80c29134d4
begin
  const DIR = "."
  hfile = CSV.File(joinpath(DIR, "data", "horse.csv"));
  horse = DataFrames.DataFrame(hfile); # convert to data frame without copying columns
  first(horse, 4)
end

# ╔═╡ 1025137e-8783-11eb-365c-ab4ec44ed329
DataFrames.describe(horse, :eltype, :first => first)

# ╔═╡ af1ff426-8783-11eb-239f-1b70da86ce69
md"""
Let's make some changes for the scientific types
"""

# ╔═╡ c5c80b78-8783-11eb-2687-1b29f016d826
begin
  coerce!(horse, autotype(horse));
  coerce!(horse, Count => Continuous);
  coerce!(horse,
        :surgery               => Multiclass,
        :age                   => Multiclass,
        :mucous_membranes      => Multiclass,
        :capillary_refill_time => Multiclass,
        :outcome               => Multiclass,
        :cp_data               => Multiclass);
  schema(horse)
end

# ╔═╡ 100e01fc-8783-11eb-1434-bb5bf50dd24d
begin
  yₕ, Xₕ = unpack(horse, ==(:outcome),
    name -> elscitype(Tables.getcolumn(horse, name)) == Continuous);
  scitype(yₕ)
end

# ╔═╡ 7b44603a-8783-11eb-1c85-37ad050e6e26
md"""
(b) Create a 70:30 `train`/`test` split of the data and train a
`LogisticClassifier` model, from the `MLJLinearModels` package, on
the `train` rows. Use `lambda=100` and default values for the
other hyper-parameters. (Although one would normally standardize
(whiten) the continuous features for this model, do not do so here.)
After training:

- (i) Recalling that a logistic classifier (aka logistic regressor) is
  a linear-based model learning a *vector* of coefficients for each
  feature (one coefficient for each target class), use the
  `fitted_params` method to find this vector of coefficients in the
  case of the `:pulse` feature. (You can convert a vector of pairs `v =
  [x1 => y1, x2 => y2, ...]` into a dictionary with `Dict(v)`.)

- (ii) Evaluate the `cross_entropy` performance on the `test`
  observations.

- (iii) In how many `test` observations does the predicted
  probability of the observed class exceed 50%?

- (iv) Find the `misclassification_rate` in the `test` set.
  (*Hint.* As this measure is deterministic, you will either  need to broadcast `mode` or use `predict_mode` instead of
  `predict`.)

"""

# ╔═╡ 7b2c41ce-8783-11eb-0e4e-77c7c1f7ca6e
begin
  trainₕ, testₕ = partition(eachindex(yₕ), 0.7)
  size(trainₕ), size(testₕ)
end

# ╔═╡ 7b13911c-8783-11eb-000a-69067847d3ca
begin
  @load LogisticClassifier pkg=MLJLinearModels
  modelₕ = MLJLinearModels.LogisticClassifier(lambda=100)
  info(modelₕ)
end

# ╔═╡ 7aeb0e5e-8783-11eb-1560-0bd4641490ab
begin
  machₕ = machine(modelₕ, Xₕ, yₕ);
  fit!(machₕ, rows=trainₕ);
end

# ╔═╡ 0ff1eba2-8783-11eb-345b-712a8df94a74
Dict(fitted_params(machₕ)[:coefs])[:pulse]

# ╔═╡ 09ddfafc-8784-11eb-2bcc-27382c2cd9b9
begin
  ŷₕ = predict(machₕ, Xₕ[testₕ,:])
  cross_entropy(ŷₕ, yₕ[testₕ]) |> mean
end

# ╔═╡ 09c5248c-8784-11eb-03b7-09739daa1528
info(modelₕ).prediction_type

# ╔═╡ 09a98e2a-8784-11eb-27ef-2f37c8ce8a05
begin
  # The predicted probs of the actual obs. in the test are given by
  pₕ = broadcast(pdf, ŷₕ, yₕ[testₕ]);

  # The number of times this probability exceeds 50% is:
  filter(p -> p > 0.5, pₕ) |> length
end

# ╔═╡ 0974d978-8784-11eb-13a2-bf91e30d3b3c
misclassification_rate(mode.(ŷₕ), yₕ[testₕ])

# ╔═╡ 2a62df22-8784-11eb-1225-9f45e08c433d
md"""
(c) Instead use a `RandomForestClassifier` model from the
    `DecisionTree` package and:

- (i) Generate an appropriate learning curve to convince yourself
  that out-of-sample estimates of the `cross_entropy` loss do not
  substantially improve for `n_trees > 50`. Use default values for
  all other hyper-parameters, and feel free to use all available
  data to generate the curve.

- (ii) Fix `n_trees=90` and use `evaluate!` to obtain a 9-fold
  cross-validation estimate of the `cross_entropy`, restricting
  sub-sampling to the `train` observations.

- (iii) Now use *all* available data but set
  `resampling=Holdout(fraction_train=0.7)` to obtain a score you can
  compare with the `KNNClassifier` in part (b)(iii). Which model is
  better?
"""

# ╔═╡ 2a366eba-8784-11eb-192b-af01d875e165
begin
  @load RandomForestClassifier pkg=DecisionTree
  model_rf = DecisionTree.RandomForestClassifier(n_trees=50)
end

# ╔═╡ 2a193a7a-8784-11eb-110d-4f0a37c0ce40
begin
  mach_rf =  machine(model_rf, Xₕ, yₕ)
  evaluate!(mach_rf, resampling=CV(nfolds=6), measure=cross_entropy)
end

# ╔═╡ 24fc87e0-8784-11eb-018f-87fee51f7013
#TODO ...

# ╔═╡ Cell order:
# ╟─2bb6648a-8774-11eb-17b7-0f6fef114c61
# ╠═58d05c3e-8774-11eb-2f17-3daba3edae13
# ╟─7240757e-8774-11eb-3c51-eb10a4c46db3
# ╟─992cc480-8785-11eb-0572-c33c1c988efc
# ╟─72108d16-8774-11eb-0003-49693e396f19
# ╟─71d64138-8774-11eb-00ff-a3d20a55ccc8
# ╠═e1fa26aa-8774-11eb-204a-a3c824321016
# ╟─4f4601e2-8785-11eb-2401-e5a978807dc6
# ╟─e1e48b42-8774-11eb-0e4e-07edf1692567
# ╠═e1c675c6-8774-11eb-0cc8-9fcd3ec74dff
# ╠═e1ad7eae-8774-11eb-0040-57d9e0fb1236
# ╟─720cb48c-8785-11eb-0649-1f4163b49de6
# ╟─e1916f20-8774-11eb-3503-ef998f66e2bd
# ╠═2f9e782a-8775-11eb-1272-efb3abb86766
# ╟─2f6944fa-8775-11eb-1111-df2b2b267f84
# ╠═2f50d958-8775-11eb-096e-0b8998321174
# ╠═2f37d368-8775-11eb-2f0c-fd13ca4fe29d
# ╠═2f1ec62a-8775-11eb-26e2-85c032960d6a
# ╠═2f0eb12c-8775-11eb-17f9-235038e0b0a7
# ╠═2ef084b8-8775-11eb-2bc0-0b949d966e90
# ╠═2ed935ba-8775-11eb-08e9-d74a578c22fa
# ╟─2ec20a02-8775-11eb-2449-5749cbbb60e1
# ╠═2ea8ed54-8775-11eb-18d6-d1d30c96d908
# ╠═e1782812-8774-11eb-1b1c-3faefea1001c
# ╠═e109c3ae-8774-11eb-3f53-1dc77c89fc61
# ╟─77d08a42-8785-11eb-0e54-fda94fc25cbe
# ╟─7e2f4de8-8775-11eb-3cc3-510cbf8258b6
# ╠═7e17a8fa-8775-11eb-267a-f9732276336b
# ╟─7df43b38-8775-11eb-1160-ab24fe881df7
# ╠═7dc2a634-8775-11eb-1a2a-2b64699a3175
# ╟─b14a0508-8775-11eb-2ea5-09289a0926b7
# ╠═7da2d596-8775-11eb-22bd-bf94ffdff871
# ╟─7d8e181a-8775-11eb-1c85-8d0179791c4d
# ╠═7d70204e-8775-11eb-1bd3-111d1dd44c45
# ╟─7d5a596c-8775-11eb-1059-c38021f064e2
# ╠═7d3c5bce-8775-11eb-242f-af3448ec4f6c
# ╟─e0397248-8775-11eb-31f1-4745600b4927
# ╠═e020dba0-8775-11eb-04c6-5791fa04a0d0
# ╟─dfebbeea-8775-11eb-0e03-19641479d64d
# ╠═dfd1df34-8775-11eb-083d-1f134983ca77
# ╟─f74b80a2-8775-11eb-3733-7f27c7d20f3d
# ╠═f7265138-8775-11eb-2873-6b37109e57df
# ╟─f70c5bca-8775-11eb-07c1-dfff7ad49f4a
# ╠═09a19f7a-8776-11eb-1203-f77b73622b87
# ╟─f6ef9698-8775-11eb-1400-7bc68f3dbe0f
# ╠═15b16a2a-8776-11eb-1592-3dce20d4a188
# ╟─15935300-8776-11eb-0c72-0f19462b9b10
# ╠═1576a3fe-8776-11eb-2e19-f7266559cf3f
# ╟─2f56cae2-8776-11eb-1244-1f49372a694b
# ╠═2f3979ec-8776-11eb-023e-c54951130c1a
# ╟─2f1fe682-8776-11eb-163d-67fa9dec2081
# ╠═dfb95914-8775-11eb-0da8-bf42ce7e7350
# ╟─dfa1167e-8775-11eb-2245-a3e41d6512ac
# ╠═6af7b84a-8776-11eb-260a-8f04dee6e55a
# ╟─6adb6884-8776-11eb-12b2-8dc0829f1bc7
# ╠═6aa68b5a-8776-11eb-2cbc-258c26de524c
# ╟─df8118d8-8775-11eb-234c-a35cf41d7430
# ╠═989e2202-8776-11eb-03a1-49fa4508f71a
# ╟─98800d12-8776-11eb-1c7c-457f310a3a4d
# ╠═9868dd68-8776-11eb-3633-afcbd37b82b6
# ╟─9850ae82-8776-11eb-33bd-4bfbb7eac41b
# ╠═98393bc6-8776-11eb-0659-f7225d5839f0
# ╟─981f0ce2-8776-11eb-1546-83d6d1ec5727
# ╠═be74448e-8776-11eb-3018-c722a9182899
# ╠═be58b4d0-8776-11eb-3eea-3350fe77c4f2
# ╟─be3edd94-8776-11eb-31c6-634f25bacbc2
# ╠═be24b926-8776-11eb-09c2-335929294063
# ╟─be14985e-8776-11eb-341f-f9c6c418f697
# ╠═bdf3586a-8776-11eb-0ac3-ddc972c36518
# ╟─bddcc460-8776-11eb-32a8-45a1bda54031
# ╠═e700cc38-8776-11eb-1c82-35ed07df74af
# ╟─e6e1d9c2-8776-11eb-243c-4db077b59c05
# ╠═e6c6f486-8776-11eb-20ed-4140d2fdfc59
# ╟─e6ad53c6-8776-11eb-37e9-f521bde89017
# ╠═98060a3a-8776-11eb-1385-ffc7eac03387
# ╠═03a66230-8777-11eb-0a4c-bda208768889
# ╟─861ed130-8785-11eb-207d-1b60a935f028
# ╟─0388dd00-8777-11eb-11b6-23a5968841f7
# ╠═1544115e-8777-11eb-28df-3167da05d7cb
# ╟─0376892a-8777-11eb-10fe-ebd9a30329b5
# ╠═03503764-8777-11eb-0793-598490c04eba
# ╟─da97b186-8781-11eb-0591-757d7c34d247
# ╠═da7f0320-8781-11eb-3047-832718ddd4d0
# ╟─da634a7c-8781-11eb-34ce-21239b737208
# ╠═da499c0a-8781-11eb-1ecc-c78a2894425e
# ╠═da312b96-8781-11eb-1a04-e709749917af
# ╟─da19bd4e-8781-11eb-2316-a77ff448d985
# ╠═03fb3db0-8782-11eb-22a3-abdefb466c7b
# ╟─f980ec9a-8784-11eb-3c1f-cbd12e90eb12
# ╟─03e5ccf8-8782-11eb-3b06-3fd79d29ece7
# ╠═03ca941a-8782-11eb-3d2c-992431e1eaf2
# ╠═03b234ea-8782-11eb-259a-0566fa974775
# ╟─70bb451e-8774-11eb-18df-512352f9b934
# ╟─704ad2e8-8774-11eb-2083-812bdb1f4ffe
# ╟─fbd65648-8774-11eb-1556-135f1f66e11a
# ╟─b8924b1c-8774-11eb-12f0-75a0c9247554
# ╠═c2ca6d30-8774-11eb-2f9f-853ae8843263
# ╠═87ccbfea-8782-11eb-28c5-1fd0c3ca976e
# ╠═877b6546-8782-11eb-0787-7f71094de4b0
# ╠═869d4842-8782-11eb-0863-77d366291f93
# ╟─a810df34-8782-11eb-1729-19864c77e6f0
# ╠═c3f6ddc0-8782-11eb-0197-25fdd80a50fd
# ╟─c2a292b0-8774-11eb-11d7-b5a05aa9f616
# ╠═c27462d2-8774-11eb-102e-135655b15bc4
# ╠═d824285c-8782-11eb-09be-7bfa9e2f63dc
# ╠═d809473a-8782-11eb-2331-19eccaff712a
# ╠═c2563758-8774-11eb-0903-dd227f7618d1
# ╠═105a8fe8-8783-11eb-2f0b-8fb550d77372
# ╟─103d8df8-8783-11eb-2c3d-eff502ac833c
# ╠═3de19b84-8783-11eb-139d-9d80c29134d4
# ╠═1025137e-8783-11eb-365c-ab4ec44ed329
# ╟─af1ff426-8783-11eb-239f-1b70da86ce69
# ╠═c5c80b78-8783-11eb-2687-1b29f016d826
# ╠═100e01fc-8783-11eb-1434-bb5bf50dd24d
# ╟─7b44603a-8783-11eb-1c85-37ad050e6e26
# ╠═7b2c41ce-8783-11eb-0e4e-77c7c1f7ca6e
# ╠═7b13911c-8783-11eb-000a-69067847d3ca
# ╠═7aeb0e5e-8783-11eb-1560-0bd4641490ab
# ╠═0ff1eba2-8783-11eb-345b-712a8df94a74
# ╠═09ddfafc-8784-11eb-2bcc-27382c2cd9b9
# ╠═09c5248c-8784-11eb-03b7-09739daa1528
# ╠═09a98e2a-8784-11eb-27ef-2f37c8ce8a05
# ╠═0974d978-8784-11eb-13a2-bf91e30d3b3c
# ╟─2a62df22-8784-11eb-1225-9f45e08c433d
# ╠═2a366eba-8784-11eb-192b-af01d875e165
# ╠═2a193a7a-8784-11eb-110d-4f0a37c0ce40
# ╠═24fc87e0-8784-11eb-018f-87fee51f7013

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
function hyperbolic_cross_grid(node_type::Function,d::S,ind::Array{S,2}) where {S <: Integer}

  multi_index = generate_multi_index(d,ind)

  # Create base nodes to be used in the sparse grid

  base_nodes = node_type(n)
  T = eltype(base_nodes)

  # Determine the unique nodes introduced at each higher level

  unique_base_nodes = Array{Array{T,1},2}(undef,size(multi_index))
  @inbounds for i in eachindex(multi_index)
    if multi_index[i] == 0
      unique_base_nodes[i] = [base_nodes[Int((n-1)/2)+1]]
    else
      unique_base_nodes[i] = [base_nodes[Int((n-1)/2)+1-multi_index[i]], base_nodes[Int((n-1)/2)+1+multi_index[i]]]
    end
  end

  # Construct the hyperbolic cross from the unique nodes

  nodes = Array{T,2}(undef,determine_grid_size(multi_index))
  l = 1
  @inbounds for j in axes(multi_index,1)
    new_nodes = unique_base_nodes[j,1]  # Here new_nodes is a 1d array
    for i = 2:d
      new_nodes = combine_nodes(new_nodes,unique_base_nodes[j,i])  # Here new_nodes becomes a 2d array
    end
    m = size(new_nodes,1)
    nodes[l:l+m-1,:] = new_nodes
    l += m
  end

  return nodes, multi_index

end

function hyperbolic_cross_grid(node_type::Function,d::S,k::S) where {S <: Integer}

  multi_index = generate_multi_index(d,k)

  n = 2*k+1

  # Create base nodes to be used in the sparse grid

  base_nodes = node_type(n)
  T = eltype(base_nodes)

  # Determine the unique nodes introduced at each higher level

  unique_base_nodes = Array{Array{T,1},2}(undef,size(multi_index))
  @inbounds for i in eachindex(multi_index)
    if multi_index[i] == 0
      unique_base_nodes[i] = [base_nodes[Int((n-1)/2)+1]]
    else
      unique_base_nodes[i] = [base_nodes[Int((n-1)/2)+1-multi_index[i]], base_nodes[Int((n-1)/2)+1+multi_index[i]]]
    end
  end

  # Construct the hyperbolic cross from the unique nodes

  nodes = Array{T,2}(undef,determine_grid_size(multi_index))
  l = 1
  @inbounds for j in axes(multi_index,1)
    new_nodes = unique_base_nodes[j,1]  # Here new_nodes is a 1d array
    for i = 2:d
      new_nodes = combine_nodes(new_nodes,unique_base_nodes[j,i])  # Here new_nodes becomes a 2d array
    end
    m = size(new_nodes,1)
    nodes[l:l+m-1,:] = new_nodes
    l += m
  end

  return nodes, multi_index

end

function hyperbolic_cross_grid(node_type::Function,d::S,k::S,n::S) where {S <: Integer}

  multi_index = generate_multi_index(d,k,n)

  # Create base nodes to be used in the sparse grid

  base_nodes = node_type(n)
  T = eltype(base_nodes)

  # Determine the unique nodes introduced at each higher level

  unique_base_nodes = Array{Array{T,1},2}(undef,size(multi_index))
  @inbounds for i in eachindex(multi_index)
    if multi_index[i] == 0
      unique_base_nodes[i] = [base_nodes[Int((n-1)/2)+1]]
    else
      unique_base_nodes[i] = [base_nodes[Int((n-1)/2)+1-multi_index[i]], base_nodes[Int((n-1)/2)+1+multi_index[i]]]
    end
  end

  # Construct the hyperbolic cross from the unique nodes

  nodes = Array{T,2}(undef,determine_grid_size(multi_index))
  l = 1
  @inbounds for j in axes(multi_index,1)
    new_nodes = unique_base_nodes[j,1]  # Here new_nodes is a 1d array
    for i = 2:d
      new_nodes = combine_nodes(new_nodes,unique_base_nodes[j,i])  # Here new_nodes becomes a 2d array
    end
    m = size(new_nodes,1)
    nodes[l:l+m-1,:] = new_nodes
    l += m
  end

  return nodes, multi_index

end

function hyperbolic_cross_grid(node_type::Function,d::S,k::S,n::Array{S,1}) where {S <: Integer}

  multi_index = generate_multi_index(d,k,n)

  # Create base nodes to be used in the sparse grid

  T = Float64
  base_nodes = Array{Array{T,1},1}(undef,d)

  for i = 1:d
    base_nodes[i] = node_type(n[i])
  end

  # Determine the unique nodes introduced at each higher level

  unique_base_nodes = Array{Array{T,1},2}(undef,size(multi_index))
  for i in axes(multi_index,1)
    for j = 1:d
      if multi_index[i,j] == 0
        unique_base_nodes[i,j] = [base_nodes[j][Int((n[j]-1)/2)+1]]
      else
        unique_base_nodes[i,j] = [base_nodes[j][Int((n[j]-1)/2)+1-multi_index[i,j]], base_nodes[j][Int((n[j]-1)/2)+1+multi_index[i,j]]]
      end
    end
  end

  # Construct the hyperbolic cross from the unique nodes

  nodes = Array{T,2}(undef,determine_grid_size(multi_index))
  l = 1
  @inbounds for j in axes(multi_index,1)
    new_nodes = unique_base_nodes[j,1]  # Here new_nodes is a 1d array
    for i = 2:d
      new_nodes = combine_nodes(new_nodes,unique_base_nodes[j,i])  # Here new_nodes becomes a 2d array
    end
    m = size(new_nodes,1)
    nodes[l:l+m-1,:] = new_nodes
    l += m
  end

  return nodes, multi_index

end

function hyperbolic_cross_grid(node_type::Function,ind::Array{S,2},domain::Union{Array{T,1},Array{T,2}}) where {S <: Integer, T <: AbstractFloat}

  if size(domain,2) != d
    error("domain is inconsistent with the number of dimensions")
  end

  (nodes, multi_index) = hyperbolic_cross_grid(node_type,ind)

  # Now scale the nodes to the desired domain

  nodes .= scale_nodes(nodes,domain)

  return nodes, multi_index

end

function hyperbolic_cross_grid(node_type::Function,d::S,k::S,domain::Union{Array{T,1},Array{T,2}}) where {S <: Integer, T <: AbstractFloat}

  if size(domain,2) != d
    error("domain is inconsistent with the number of dimensions")
  end

  (nodes, multi_index) = hyperbolic_cross_grid(node_type,d,k)

  # Now scale the nodes to the desired domain

  nodes .= scale_nodes(nodes,domain)

  return nodes, multi_index

end

function hyperbolic_cross_grid(node_type::Function,d::S,k::S,n::Union{S,Array{S,1}},domain::Union{Array{T,1},Array{T,2}}) where {S <: Integer, T <: AbstractFloat}

  if size(domain,2) != d
    error("domain is inconsistent with the number of dimensions")
  end

  (nodes, multi_index) = hyperbolic_cross_grid(node_type,d,k,n)

  # Now scale the nodes to the desired domain

  nodes .= scale_nodes(nodes,domain)

  return nodes, multi_index

end


text:
#####
##### Sequential Access Kernels
#####

# Conver things to 'Vals' that aren't already Vals
val(x::Val) = x
val(x) = Val(x)

"""
    sequential_read(A::AbstractArray, ::Val{N}, [nontemporal = Val(false), aligned = Val(true)]

Perform a sequential read on the contents of `A` using vectors intrinsics of size `N`.

Optional arguments `nontemporal` and `aligned` should be passes as `Val{True}()` or
`Val{False}()` and control if the vector instructions are nontemporal and aligned,
respectively.
"""
function sequential_read(
        A::AbstractArray{T},
        ::Val{N},
        nontemporal = Val(false),
        aligned = Val(true),
    ) where {T, N}

    kp = KernelParam(;
        kernel = ReadOnly,
        loadtype = nontemporal == Val(true) ? Nontemporal : Standard,
        eltype = T,
        iterator = Sequential,
        vectorsize = N,
        unroll = 4
    )
    return execute!(A, kp)
end

"""
    sequential_write(A::AbstractArray, ::Val{N}, [nontemporal = Val(false), aligned = Val(true)]

Perform a sequential read on the contents of `A` using vectors intrinsics of size `N`.

Optional arguments `nontemporal` and `aligned` should be passes as `Val{True}()` or
`Val{False}()` and control if the vector instructions are nontemporal and aligned,
respectively.
"""
function sequential_write(
        A::AbstractArray{T},
        ::Val{N},
        nontemporal = Val(false),
        aligned = Val(true),
    ) where {T, N}

    kp = KernelParam(;
        kernel = WriteOnly,
        storetype = nontemporal == Val(true) ? Nontemporal : Standard,
        eltype = T,
        iterator = Sequential,
        vectorsize = N,
        unroll = 4
    )
    return execute!(A, kp)
end

"""
    sequential_readwrite(A::AbstractArray, ::Val{N}, [nontemporal = Val(false), aligned = Val(true)]

Perform a sequential read on the contents of `A` using vectors intrinsics of size `N`.

Optional arguments `nontemporal` and `aligned` should be passes as `Val{True}()` or
`Val{False}()` and control if the vector instructions are nontemporal and aligned,
respectively.
"""
function sequential_readwrite(
        A::AbstractArray{T},
        ::Val{N},
        nontemporal = Val(false),
        aligned = Val(true),
    ) where {T, N}

    kp = KernelParam(;
        kernel = ReadWrite,
        loadtype = nontemporal == Val(true) ? Nontemporal : Standard,
        storetype = nontemporal == Val(true) ? Nontemporal : Standard,
        eltype = T,
        iterator = Sequential,
        vectorsize = N,
        unroll = 4,
    )
    return execute!(A, kp)
end

#####
##### Random Access Kernels
#####

"""
    random_read(A::AbstractArray, ::Val{N}, [nontemporal = Val(false), aligned = Val(true)]

Perform a random read on the contents of `A` using vectors intrinsics of size `N`.

Optional arguments `nontemporal` and `aligned` should be passes as `Val{True}()` or
`Val{False}()` and control if the vector instructions are nontemporal and aligned,
respectively.
"""
function random_read(
        A::AbstractArray{T},
        ::Val{N},
        nontemporal = Val(false),
        aligned = Val(true),
        unroll::Val{U} = Val(1),
    ) where {T,N,U}

    kp = KernelParam(;
        kernel = ReadOnly,
        loadtype = nontemporal == Val(true) ? Nontemporal : Standard,
        eltype = T,
        iterator = PseudoRandom,
        vectorsize = N,
        unroll = U,
    )
    return execute!(A, kp)
end

"""
    random_write(A::AbstractArray, ::Val{N}, [nontemporal = Val(false), aligned = Val(true)]

Perform a random read on the contents of `A` using vectors intrinsics of size `N`.

Optional arguments `nontemporal` and `aligned` should be passes as `Val{True}()` or
`Val{False}()` and control if the vector instructions are nontemporal and aligned,
respectively.
"""
function random_write(
        A::AbstractArray{T},
        ::Val{N},
        nontemporal = Val{false}(),
        aligned = Val(true),
        unroll::Val{U} = Val(1),
    ) where {T,N,U}

    kp = KernelParam(;
        kernel = WriteOnly,
        loadtype = nontemporal == Val(true) ? Nontemporal : Standard,
        eltype = T,
        iterator = PseudoRandom,
        vectorsize = N,
        unroll = U,
    )
    return execute!(A, kp)
end

"""
    random_readwrite(A::AbstractArray, ::Val{N}, [nontemporal = Val(false), aligned = Val(true)]

Perform a random read on the contents of `A` using vectors intrinsics of size `N`.

Optional arguments `nontemporal` and `aligned` should be passes as `Val{True}()` or
`Val{False}()` and control if the vector instructions are nontemporal and aligned,
respectively.
"""
function random_readwrite(
        A::AbstractArray{T},
        ::Val{N},
        nontemporal = Val(false),
        aligned = Val(true),
        unroll::Val{U} = Val(1),
    ) where {T,N,U}

    kp = KernelParam(;
        kernel = ReadWrite,
        loadtype = nontemporal == Val(true) ? Nontemporal : Standard,
        eltype = T,
        iterator = PseudoRandom,
        vectorsize = N,
        unroll = U,
    )
    return execute!(A, kp)
end


text:
using Roots
using Statistics
using FileIO
using Parameters
include("random_matrices.jl")
include("network_sims_gpu.jl")

etas = collect(-1.0f0:0.1f0:1.0f0)
ind = parse(Int64, ARGS[1])

N = 3000
eta = etas[ind]
s = 0.0f0

J_seeds = collect(5000:5019)
h0_seed = 6000

h0 = generate_nv(N, Float32, seed=h0_seed)

wait = 200.0f0
t_bound = 2000.0f0
dt = 0.01f0

startstep = round(Int, wait / dt) + 1
lateststep = round(Int, t_bound / dt)

@unpack target_avg_var = load("vanilla_target_avg_var.jld2")

# Difference between the average variance of the network ensemble and target
function f(g)
    vanilla_network_avg_vars = Vector{Float64}(undef, length(J_seeds))
    for (i, J_seed) in enumerate(J_seeds)
        J = generate_J(N, eta, seed=J_seed)
        A = sim_activity_gpu(J, g, h0; s=s, dt=dt, t_bound=t_bound,
                minimum_norm=1e-90)
        endstep = min(lateststep, size(A, 2))
        avg_var = mean(var(view(A, :, startstep:endstep), dims=2))
        vanilla_network_avg_vars[i] = avg_var
    end
    mean(vanilla_network_avg_vars) - target_avg_var
end

lo_bound = 0.0f0
hi_bound = 100f0
@time g = find_root(f, (lo_bound, hi_bound), Bisection(), verbose=true)

out = Dict{String,Any}()
@pack! out = N, g, eta, s, J_seeds, h0_seed, wait, t_bound, dt, target_avg_var
save("Rootfinding_E$(ind)_N=$N,eta=$eta,s=$s.jld2", out)


text:
#!/usr/bin/env julia --project=".."
using RobotOS
@rosimport julia_tutorials.srv.AddTwoInts
rostypegen()
using .julia_tutorials.srv: AddTwoInts, AddTwoIntsResponse, AddTwoIntsRequest
function add_two_ints_client(x, y)
    wait_for_service("add_two_ints")
    try
        add_two_ints = ServiceProxy{AddTwoInts}("add_two_ints")
        respl = add_two_ints(AddTwoIntsRequest(x, y)) # got an AddTwoIntsResponse
        return respl.sum
    catch err
        println("Service call failed: $err")
    end
end
usage() = "$PROGRAM_FILE [x y]"
if !isinteractive()
    if length(ARGS) == 2
        x = parse(Int, ARGS[1])
        y = parse(Int, ARGS[2])
    else
        println(usage())
        exit(1)
    end
    println("Requesting $x + $y")
    println("$x + $y = $(add_two_ints_client(x, y))")
end

max_stars_repo_name: mattborghi/ros-julia-examples

max_issues_repo_name: mattborghi/ros-julia-examples

max_forks_repo_name: mattborghi/ros-julia-examples


text:
#!/usr/bin/julia

#include("./graphtest.jl")

using LightGraphs

using NetCDF

filename = "../netcdf_results/nhept.nc"
println(ncinfo(filename))
x = ncread(filename, "source")

println(" loadit ")

A = [
           0 1 1
           1 0 1
           1 0.4 0
       ]

G = DiGraph(A)

max_stars_repo_name: wolfiex/VisACC

max_issues_repo_name: wolfiex/VisACC

max_forks_repo_name: wolfiex/VisACC


text:
@testset ExtendedTestSet "dice_metric" begin
    @testset ExtendedTestSet "dice_metric" begin
        x = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        y = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        @test dice_metric(x, y) ≈ 1
    end

    @testset ExtendedTestSet "dice_metric" begin
        x = [0 0 0 1; 0 0 0 1; 0 0 0 1; 0 0 0 1]
        y = [1 1 1 0; 1 1 1 0; 1 1 1 0; 1 1 1 0]
        @test dice_metric(x, y) ≈ 0
    end
end


text:
using Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

#
# Setup the following communication graph
#
#   +-----+
#   |     |
#   v     v
#   0<-+  3
#   ^  |  ^
#   |  |  |
#   v  |  v
#   1  +--2
#   ^     |
#   |     |
#   +-----+
#
#

if rank == 0
    dest   = Cint[1,3]
    degree = Cint[length(dest)]
elseif rank == 1
    dest   = Cint[0]
    degree = Cint[length(dest)]
elseif rank == 2
    dest   = Cint[3,0,1]
    degree = Cint[length(dest)]
elseif rank == 3
    dest   = Cint[0,2,1]
    degree = Cint[length(dest)]
end

source = Cint[rank]
graph_comm = MPI.Dist_graph_create(comm, source, degree, dest)

#
# Now send the rank across the edges.
#
# Version 1: use allgather primitive
#

send = [rank]
if rank == 0
    recv = [-1, -1, -1]
elseif rank == 1
    recv = [-1, -1, -1]
elseif rank == 2
    recv = [-1]
elseif rank == 3
    recv = [-1, -1]
end

MPI.Neighbor_allgather!(send, recv, graph_comm);

println("rank = $(rank): $(recv)")

#
# Version 2: use alltoall primitive
#

if rank == 0
    send   = [rank, rank]
    recv   = [-1, -1, -1]
elseif rank == 1
    send = [rank]
    recv   = [-1, -1, -1]
elseif rank == 2
    send = [rank, rank, rank]
    recv   = [-1]
elseif rank == 3
    send = [rank, rank, rank]
    recv   = [-1, -1]
end

MPI.Neighbor_alltoall!(UBuffer(send,1), UBuffer(recv,1), graph_comm);

println("rank = $(rank): $(recv)")

#
# Now send the rank exactly rank times across the edges.
#
# Rank i receives i+1 values from each adjacent process
if rank == 0
    send       = [rank, rank,
                  rank, rank, rank, rank]
    send_count = [2, 4]

    recv   = [-1, -1, -1]
    recv_count = [1, 1, 1]
elseif rank == 1
    send       = [rank]
    send_count = [1]

    recv       = [-1, -1, -1, -1, -1, -1]
    recv_count = [2, 2, 2]
elseif rank == 2
    send       = [rank, rank, rank, rank,
                  rank,
                  rank,rank]
    send_count = [4, 1, 2]

    recv       = [-1, -1, -1]
    recv_count = [3]
elseif rank == 3
    send       = [rank,
                  rank, rank,rank,
                  rank, rank]
    send_count = [1, 3, 2]

    recv       = [-1, -1, -1, -1, -1, -1, -1, -1]
    recv_count = [4, 4]
end

MPI.Neighbor_alltoallv!(VBuffer(send,send_count), VBuffer(recv,recv_count), graph_comm);
println("rank = $(rank): $(recv)")

MPI.Finalize()


text:
@testset "Static graphs" begin
    g = @inferred(CompleteDiGraph(5))
    @test nv(g) == 5 && ne(g) == 20
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = CompleteDiGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = CompleteDiGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = CompleteDiGraph(2)
    @test nv(g) == 2 && ne(g) == 2
    g = @inferred CompleteDiGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 127 * (127 - 1)
    @test eltype(g) == Int8

    g = @inferred(CompleteGraph(5))
    @test nv(g) == 5 && ne(g) == 10
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = CompleteGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = CompleteGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = CompleteGraph(2)
    @test nv(g) == 2 && ne(g) == 1
    g = @inferred CompleteGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 127 * (127 - 1) ÷ 2
    @test eltype(g) == Int8

    g = @inferred(CompleteBipartiteGraph(5, 8))
    @test nv(g) == 13 && ne(g) == 40
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = CompleteBipartiteGraph(0,0)
    @test nv(g) == 0 && ne(g) == 0
    g = CompleteBipartiteGraph(5,0)
    @test nv(g) == 5 && ne(g) == 0
    g = CompleteBipartiteGraph(0,5)
    @test nv(g) == 5 && ne(g) == 0
    g = @inferred CompleteBipartiteGraph(Int8(100), Int8(27))
    @test nv(g) == 127 && ne(g) == 100 * 27
    @test eltype(g) == Int8
    g = CompleteBipartiteGraph(Int8(127), Int8(0))
    @test nv(g) == 127 && ne(g) == 0
    g = CompleteBipartiteGraph(Int8(0), Int8(127))
    @test nv(g) == 127 && ne(g) == 0

    g = @inferred(StarGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = StarGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = StarGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = StarGraph(2)
    @test nv(g) == 2 && ne(g) == 1
    g = @inferred StarGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 127 - 1
    @test eltype(g) == Int8

    g = @inferred(StarDiGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = StarDiGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = StarDiGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = StarDiGraph(2)
    @test nv(g) == 2 && ne(g) == 1
    @test first(edges(g)) == Edge(1, 2) # edges should point outwards from vertex 1
    g = @inferred StarDiGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 127 - 1
    @test eltype(g) == Int8

    g = @inferred(PathDiGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = PathDiGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = PathDiGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = @inferred PathDiGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 126
    @test eltype(g) == Int8

    g = @inferred(PathGraph(5))
    @test nv(g) == 5 && ne(g) == 4
    @test isvalid_simplegraph(g)
    # tests for extreme values
    g = PathGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = PathGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = @inferred PathGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 126
    @test eltype(g) == Int8

    g = @inferred(CycleDiGraph(5))
    @test nv(g) == 5 && ne(g) == 5
    @test isvalid_simplegraph(g)
     # tests for extreme values
    g = CycleDiGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = CycleDiGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = CycleDiGraph(2)
    @test nv(g) == 2 && ne(g) == 2
    g = @inferred CycleDiGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 127
    @test eltype(g) == Int8

    g = @inferred(CycleGraph(5))
    @test nv(g) == 5 && ne(g) == 5
    @test isvalid_simplegraph(g)
     # tests for extreme values
    g = CycleGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = CycleGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = CycleGraph(2)
    @test nv(g) == 2 && ne(g) == 1
    g = @inferred CycleGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 127
    @test eltype(g) == Int8

    g = @inferred(WheelDiGraph(5))
    @test nv(g) == 5 && ne(g) == 8
    @test isvalid_simplegraph(g)
      # tests for extreme values
    g = WheelDiGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = WheelDiGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = WheelDiGraph(2)
    @test nv(g) == 2 && ne(g) == 1
    g = WheelDiGraph(3)
    @test nv(g) == 3 && ne(g) == 4
    g = @inferred WheelDiGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 2 * 126
    @test eltype(g) == Int8

    g = @inferred(WheelGraph(5))
    @test nv(g) == 5 && ne(g) == 8
    @test isvalid_simplegraph(g)
      # tests for extreme values
    g = WheelGraph(0)
    @test nv(g) == 0 && ne(g) == 0
    g = WheelGraph(1)
    @test nv(g) == 1 && ne(g) == 0
    g = WheelGraph(2)
    @test nv(g) == 2 && ne(g) == 1
    g = WheelGraph(3)
    @test nv(g) == 3 && ne(g) == 3
    g = @inferred WheelGraph(Int8(127))
    @test nv(g) == 127 && ne(g) == 2 * 126
    @test eltype(g) == Int8

    g = @inferred(Grid([3, 3, 4]))
    @test nv(g) == 3 * 3 * 4
    @test ne(g) == 75
    @test Δ(g) == 6
    @test δ(g) == 3
    @test isvalid_simplegraph(g)

    g = @inferred(Grid([3, 3, 4], periodic=true))
    @test nv(g) == 3 * 3 * 4
    @test ne(g) == 108
    @test Δ(g) == 6
    @test δ(g) == 6
    @test isvalid_simplegraph(g)

    g = @inferred(CliqueGraph(3, 5))
    @test nv(g) == 15 && ne(g) == 20
    @test g[1:3] == CompleteGraph(3)
    @test isvalid_simplegraph(g)

    g = @inferred(crosspath(3, BinaryTree(2)))
    # f = Vector{Vector{Int}}[[2 3 4];
    # [1 5];
    # [1 6];
    # [1 5 6 7];
    # [2 4 8];
    # [3 4 9];
    # [4 8 9];
    # [5 7];
    # [6 7]
    # ]
    I = [1, 1, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 9]
    J = [2, 3, 4, 1, 5, 1, 6, 1, 5, 6, 7, 2, 4, 8, 3, 4, 9, 4, 8, 9, 5, 7, 6, 7]
    V = ones(Int, length(I))
    Adj = sparse(I, J, V)
    @test Adj == sparse(g)
    @test isvalid_simplegraph(g)
    @test_throws DomainError BinaryTree(Int8(8))

    g = @inferred(DoubleBinaryTree(3))
    # [[3, 2, 8]
    # [4, 1, 5]
    # [1, 6, 7]
    # [2]
    # [2]
    # [3]
    # [3]
    # [10, 9, 1]
    # [11, 8, 12]
    # [8, 13, 14]
    # [9]
    # [9]
    # [10]
    # [10]]
    I = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 5, 6, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 12, 13, 14]
    J = [3, 2, 8, 4, 1, 5, 1, 6, 7, 2, 2, 3, 3, 10, 9, 1, 11, 8, 12, 8, 13, 14, 9, 9, 10, 10]
    V = ones(Int, length(I))
    Adj = sparse(I, J, V)
    @test Adj == sparse(g)
    @test isvalid_simplegraph(g)

    rg3 = @inferred(RoachGraph(3))
    # [3]
    # [4]
    # [1, 5]
    # [2, 6]
    # [3, 7]
    # [4, 8]
    # [9, 8, 5]
    # [10, 7, 6]
    # [11, 10, 7]
    # [8, 9, 12]
    # [9, 12]
    # [10, 11]
    I = [1, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 12, 12]
    J = [3, 4, 1, 5, 2, 6, 3, 7, 4, 8, 9, 8, 5, 10, 7, 6, 11, 10, 7, 8, 9, 12, 9, 12, 10, 11]
    V = ones(Int, length(I))
    Adj = sparse(I, J, V)
    @test Adj == sparse(rg3)
    @test isvalid_simplegraph(g)
end


text:
type Mesh
    nodes::Array{Node}
    local_nodes::Array{Node}
    elements::Array{Element}
    local_elements::Array{Element}
    boundary_info::BoundaryInfo

    # Map Node IDs to the elements they are connected to
    node_to_elem_map::Dict{Int64, Array{Element}}

    Mesh() = new(Array{Node}(0), Array{Node}(0), Array{Element}(0), Array{Element}(0), BoundaryInfo(), Dict{Int64, Array{Element}}())
    Mesh(nodes::Array{Node}, elements::Array{Element}) = new(nodes, Array{Node}(0), elements, Array{Element}(0), BoundaryInfo(), Dict{Int64, Array{Element}}())
end


text:
# Optimization and solution query.

"""
Call Knitro to solve the problem. The return value indicates
the solution status.

"""
function KN_solve(m::Model)
    # Check sanity. If model has Julia callbacks, we need to ensure
    # that Knitro is not multithreaded. Otherwise, the code will segfault
    # as we have trouble calling Julia code from multithreaded C
    # code. See issue #93 on https://github.com/jump-dev/KNITRO.jl.
    if has_callbacks(m)
        KN_set_param(m, "par_numthreads", 1)
        KN_set_param(m, "par_msnumthreads", 1)
    end
    # For KN_solve, we do not return an error if ret is different of 0.
    m.status = @kn_ccall(solve, Cint, (Ptr{Cvoid},), m.env)
    return m.status
end

##################################################
# Solution status and info
##################################################
"Return the solution status, objective, primal and dual variables."
function KN_get_solution(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    nx = KN_get_number_vars(m)
    nc = KN_get_number_cons(m)

    x = zeros(Cdouble, nx)
    lambda = zeros(Cdouble, nx + nc)
    status = Cint[0]
    obj = Cdouble[0.]

    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, x, lambda)

    # Keep solution in cache.
    m.status = status[1]
    m.x = x
    m.mult = lambda
    m.obj_val = obj[1]
    return status[1], obj[1], x, lambda
end

# some wrapper functions for MOI
function get_status(m::Model)
    @assert m.env != C_NULL
    if m.status != 1
        return m.status
    end
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, C_NULL, C_NULL)
    # Keep status in cache.
    m.status = status[1]
    return status[1]
end

function get_objective(m::Model)
    @assert m.env != C_NULL
    if isfinite(m.obj_val)
        return m.obj_val
    end
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, C_NULL, C_NULL)
    # Keep objective value in cache.
    m.obj_val = obj[1]
    return obj[1]
end

function get_solution(m::Model)
    # We first check that the model is well defined to avoid segfault.
    @assert m.env != C_NULL
    if !isempty(m.x)
        return m.x
    end
    nx = KN_get_number_vars(m)
    x = zeros(Cdouble, nx)
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, x, C_NULL)
    # Keep solution in cache.
    m.x = x
    return x
end
get_solution(m::Model, ix::Int) = isempty(m.x) ? get_solution(m)[ix] : m.x[ix]

function get_dual(m::Model)
    # we first check that the model is well defined to avoid segfault
    @assert m.env != C_NULL
    if !isempty(m.mult)
        return m.mult
    end
    nx = KN_get_number_vars(m)
    nc = KN_get_number_cons(m)
    lambda = zeros(Cdouble, nx + nc)
    status = Cint[0]
    obj = Cdouble[0.]
    ret = @kn_ccall(get_solution, Cint,
                    (Ptr{Cvoid}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, status, obj, C_NULL, lambda)
    # Keep multipliers in cache.
    m.mult = lambda
    return lambda
end
get_dual(m::Model, ix::Int) = isempty(m.mult) ? get_dual(m)[ix] : m.mult[ix]

# New getters for Knitro >= 12.0
if KNITRO_VERSION >= v"12.0"
    @define_getters get_var_primal_values
    @define_getters get_var_dual_values
    @define_getters get_con_dual_values
end


text:
### StochasticPolicy ###
# maintained by @etotheipluspi

"""
   StochasticPolicy{D, RNG <: AbstractRNG}

Represents a stochastic policy. Action are sampled from an arbitrary distribution.

Constructor:

    `StochasticPolicy(distribution; rng=Random.GLOBAL_RNG)`

# Fields
- `distribution::D`
- `rng::RNG` a random number generator
"""
mutable struct StochasticPolicy{D, RNG <: AbstractRNG} <: Policy
    distribution::D
    rng::RNG
end
# The constructor below should be used to create the policy so that the action space is initialized correctly
StochasticPolicy(distribution; rng=Random.GLOBAL_RNG) = StochasticPolicy(distribution, rng)

## policy execution ##
function action(policy::StochasticPolicy, s)
    return rand(policy.rng, policy.distribution)
end

## convenience functions ##
updater(policy::StochasticPolicy) = VoidUpdater() # since the stochastic policy does not depend on the belief

# Samples actions uniformly
UniformRandomPolicy(problem, rng=Random.GLOBAL_RNG) = StochasticPolicy(actions(problem), rng)

"""
    CategoricalTabularPolicy

represents a stochastic policy sampling an action from a categorical distribution with weights given by a `ValuePolicy`

constructor:

`CategoricalTabularPolicy(mdp::Union{POMDP,MDP}; rng=Random.GLOBAL_RNG)`

# Fields
- `stochastic::StochasticPolicy`
- `value::ValuePolicy`
"""
mutable struct CategoricalTabularPolicy <: Policy
    stochastic::StochasticPolicy
    value::ValuePolicy
end
function CategoricalTabularPolicy(mdp::Union{POMDP,MDP}; rng=Random.GLOBAL_RNG)
    CategoricalTabularPolicy(StochasticPolicy(Weights(zeros(length(actions((mdp))))), rng), ValuePolicy(mdp))
end

function action(policy::CategoricalTabularPolicy, s)
    policy.stochastic.distribution = Weights(policy.value.value_table[stateindex(policy.value.mdp, s),:])
    return policy.value.act[sample(policy.stochastic.rng, policy.stochastic.distribution)]
end

max_stars_repo_name: JuliaPOMDP/POMDPPolicies

max_issues_repo_name: JuliaPOMDP/POMDPPolicies

max_forks_repo_name: JuliaPOMDP/POMDPPolicies


text:
type GlmResp{V<:FPVector,D<:UnivariateDistribution,L<:Link} <: ModResp       # response in a glm model
    y::V                                       # response
    d::D
    l::L
    devresid::V                                # (squared) deviance residuals
    eta::V                                     # linear predictor
    mu::V                                      # mean response
    mueta::V                                   # derivative of mu w.r.t. eta
    offset::V                                  # offset added to linear predictor (usually 0)
    var::V                                     # (unweighted) variance at current mu
    wts::V                                     # prior weights
    wrkwts::V                                  # working weights
    wrkresid::V                                # working residuals
    function GlmResp(y::V, d::D, l::L, eta::V, mu::V, off::V, wts::V)
        if isa(d, Binomial)
            for yy in y
                0. <= yy <= 1. || error("$yy in y is not in [0,1]")
            end
        else
            for yy in y
                insupport(d, yy) || error("y must be in the support of d")
            end
        end
        n = length(y)
        length(eta) == length(mu) == length(wts) == n || error("mismatched sizes")
        lo = length(off)
        lo == 0 || lo == n || error("offset must have length $n or length 0")
        res = new(y,d,l,similar(y),eta,mu,similar(y),off,similar(y),wts,similar(y),similar(y))
        updatemu!(res, eta)
        res
    end
end

# returns the sum of the squared deviance residuals
deviance(r::GlmResp) = sum(r.devresid)

function updatemu!{T<:FPVector}(r::GlmResp{T}, linPr::T)
    y = r.y
    dist = r.d
    link = r.l
    eta = r.eta
    mu = r.mu
    muetav = r.mueta
    offset = r.offset
    var = r.var
    wts = r.wts
    wrkresid = r.wrkresid
    devresidv = r.devresid

    if isempty(offset)
        copy!(eta, linPr)
    else
        broadcast!(+, eta, linPr, offset)
    end

    @inbounds @simd for i = eachindex(eta,mu,muetav,var,y,wrkresid,devresidv)
        η = eta[i]

        # apply the inverse link function generating the mean vector (μ) from the linear predictor (η)
        μ = mu[i] = linkinv(link, η)

        # evaluate the mueta vector (derivative of μ w.r.t. η) from the linear predictor (eta)
        dμdη = muetav[i] = mueta(link, η)

        var[i] = glmvar(dist, link, μ, η)
        ys = y[i]
        wrkresid[i] = (ys - μ)/dμdη
        devresidv[i] = devresid(dist, ys, μ, wts[i])
    end
    r
end

function wrkresp(r::GlmResp)
    tmp = r.eta + r.wrkresid
    isempty(r.offset) ? tmp : broadcast!(-, tmp, tmp, r.offset)
end

function wrkwt!(r::GlmResp)
    wrkwts = r.wrkwts
    mueta = r.mueta
    var = r.var
    if isempty(r.wts)
        @simd for i = eachindex(var,wrkwts,mueta)
            @inbounds wrkwts[i] = abs2(mueta[i])/var[i]
        end
    else
        wts = r.wts
        @simd for i = eachindex(var,wrkwts,wts,mueta,var)
            @inbounds wrkwts[i] = wts[i] * abs2(mueta[i])/var[i]
        end
    end
    wrkwts
end

abstract AbstractGLM <: LinPredModel

type GeneralizedLinearModel{G<:GlmResp,L<:LinPred} <: AbstractGLM
    rr::G
    pp::L
    fit::Bool
end

function coeftable(mm::AbstractGLM)
    cc = coef(mm)
    se = stderr(mm)
    zz = cc ./ se
    CoefTable(hcat(cc,se,zz,2.0 * ccdf(Normal(), abs.(zz))),
              ["Estimate","Std.Error","z value", "Pr(>|z|)"],
              ["x$i" for i = 1:size(mm.pp.X, 2)], 4)
end

function confint(obj::AbstractGLM, level::Real)
    hcat(coef(obj),coef(obj)) + stderr(obj)*quantile(Normal(),(1. -level)/2.)*[1. -1.]
end
confint(obj::AbstractGLM) = confint(obj, 0.95)

deviance(m::AbstractGLM) = deviance(m.rr)

function loglikelihood(m::AbstractGLM)
    r = m.rr
    wts = r.wts
    y = r.y
    mu = r.mu
    ϕ = deviance(m)/sum(wts)
    d = r.d
    ll = zero(one(eltype(wts)) * one(loglik_obs(d, y[1], mu[1], wts[1], ϕ)))
    @inbounds for i in eachindex(y, mu, wts)
        ll += loglik_obs(d, y[i], mu[i], wts[i], ϕ)
    end
    ll
end

df(x::GeneralizedLinearModel) = dispersion_parameter(x.rr.d) ? length(coef(x)) + 1 : length(coef(x))

function _fit!(m::AbstractGLM, verbose::Bool, maxIter::Integer, minStepFac::Real,
              convTol::Real, start)
    m.fit && return m
    maxIter >= 1 || error("maxIter must be positive")
    0 < minStepFac < 1 || error("minStepFac must be in (0, 1)")

    cvg, p, r = false, m.pp, m.rr
    lp = r.mu
    if start != nothing
        copy!(p.beta0, start)
        fill!(p.delbeta, 0)
        linpred!(lp, p, 0)
        updatemu!(r, lp)
    else
        delbeta!(p, wrkresp(r), wrkwt!(r))
        linpred!(lp, p)
        updatemu!(r, lp)
        installbeta!(p)
    end
    devold = deviance(m)
    for i=1:maxIter
        f = 1.0
        local dev
        try
            delbeta!(p, r.wrkresid, wrkwt!(r))
            linpred!(lp, p)
            updatemu!(r, lp)
            dev = deviance(m)
        catch e
            isa(e, DomainError) ? (dev = Inf) : rethrow(e)
        end
        while dev > devold
            f /= 2.
            f > minStepFac || error("step-halving failed at beta0 = $(p.beta0)")
            try
                updatemu!(r, linpred(p, f))
                dev = deviance(m)
            catch e
                isa(e, DomainError) ? (dev = Inf) : rethrow(e)
            end
        end
        installbeta!(p, f)
        crit = (devold - dev)/dev
        verbose && println("$i: $dev, $crit")
        if crit < convTol
            cvg = true
            break
        end
        devold = dev
    end
    cvg || error("failure to converge in $maxIter iterations")
    m.fit = true
    m
end

StatsBase.fit!(m::AbstractGLM; verbose::Bool=false, maxIter::Integer=30,
              minStepFac::Real=0.001, convTol::Real=1.e-6, start=nothing) =
    _fit!(m, verbose, maxIter, minStepFac, convTol, start)

function initialeta!(dist::UnivariateDistribution, link::Link,
                     eta::AbstractVector, y::AbstractVector, wts::AbstractVector,
                     off::AbstractVector)
    length(eta) == length(y) == length(wts) || throw(DimensionMismatch("argument lengths do not match"))
    @inbounds @simd for i = eachindex(y,eta,wts)
        μ = mustart(dist, y[i], wts[i])
        eta[i] = linkfun(link, μ)
    end
    if !isempty(off)
        @inbounds @simd for i = eachindex(eta,off)
            eta[i] -= off[i]
        end
    end
    eta
end

function StatsBase.fit!(m::AbstractGLM, y; wts=nothing, offset=nothing, dofit::Bool=true,
                        verbose::Bool=false, maxIter::Integer=30, minStepFac::Real=0.001, convTol::Real=1.e-6,
                        start=nothing)
    r = m.rr
    V = typeof(r.y)
    r.y = copy!(r.y, y)
    isa(wts, Void) || copy!(r.wts, wts)
    isa(offset, Void) || copy!(r.offset, offset)
    initialeta!(r.d, r.l, r.eta, r.y, r.wts, r.offset)
    updatemu!(r, r.eta)
    fill!(m.pp.beta0, 0)
    m.fit = false
    if dofit
        _fit!(m, verbose, maxIter, minStepFac, convTol, start)
    else
        m
    end
end

function fit{M<:AbstractGLM,T<:FP,V<:FPVector}(::Type{M},
    X::Union{Matrix{T},SparseMatrixCSC{T}}, y::V,
    d::UnivariateDistribution,
    l::Link = canonicallink(d);
    dofit::Bool = true,
    wts::V = ones(y),
    offset::V = similar(y, 0), fitargs...)

    size(X, 1) == size(y, 1) || throw(DimensionMismatch("number of rows in X and y must match"))
    n = length(y)
    length(wts) == n || throw(DimensionMismatch("length(wts) does not match length(y)"))
    if length(offset) != n && length(offset) != 0
        throw(DimensionMismatch("length(offset) does not match length(y)"))
    end

    wts = T <: Float64 ? copy(wts) : convert(typeof(y), wts)
    off = T <: Float64 ? copy(offset) : convert(Vector{T}, offset)
    eta = initialeta!(d, l, similar(y), y, wts, off)
    rr = GlmResp{typeof(y), typeof(d), typeof(l)}(y, d, l, eta, similar(y), offset, wts)
    res = M(rr, cholpred(X), false)
    dofit ? fit!(res; fitargs...) : res
end

fit{M<:AbstractGLM}(::Type{M},
    X::Union{Matrix,SparseMatrixCSC},
    y::AbstractVector,
    d::UnivariateDistribution,
    l::Link=canonicallink(d); kwargs...) =
    fit(M, float(X), float(y), d, l; kwargs...)

glm(X, y, args...; kwargs...) = fit(GeneralizedLinearModel, X, y, args...; kwargs...)

"""
    dispersion(m::AbstractGLM, sqr::Bool=false)

    Estimated dispersion (or scale) parameter for a model's distribution,
    generally written σ for linear models and ϕ for generalized linear models.
    It is by definition equal to 1 for Binomial and Poisson families.

    If `sqr` is `true`, the squared parameter is returned.
"""
function dispersion(m::AbstractGLM, sqr::Bool=false)
    wrkwts = m.rr.wrkwts
    wrkresid = m.rr.wrkresid

    if isa(m.rr.d, Union{Binomial, Poisson})
        return one(eltype(wrkwts))
    end

    s = zero(eltype(wrkwts))
    @inbounds @simd for i = eachindex(wrkwts,wrkresid)
        s += wrkwts[i]*abs2(wrkresid[i])
    end
    s /= df_residual(m)
    sqr ? s : sqrt(s)
end

## Prediction function for GLMs
function predict(mm::AbstractGLM, newX::AbstractMatrix; offset::FPVector=Array(eltype(newX),0))
    eta = newX * coef(mm)
    if length(mm.rr.offset) > 0
        length(offset) == size(newX, 1) ||
            throw(ArgumentError("fit with offset, so `offset` kw arg must be an offset of length `size(newX, 1)`"))
        broadcast!(+, eta, eta, offset)
    else
        length(offset) > 0 && throw(ArgumentError("fit without offset, so value of `offset` kw arg does not make sense"))
    end
    mu = [linkinv(mm.rr.l, x) for x in eta]
end


text:
@testset "Physical Constants" begin
    @test InteratomicPotentials.ϵ₀ ≈ austrip(8.8541878128e-12u"F / m")
    @test InteratomicPotentials.kₑ ≈ austrip(8.9875517923e9u"N * m^2 / C^2")
end