text:
@enum HeadType rounded flat headless tslot
@enum DriveStyle hex phillips slotted torx
@enum Material aluminium brass steel plastic wood

struct Screw
    prodnum::Int
    headtype::HeadType
    drivestyle::DriveStyle
    material::Material
end

#list enum instances
instances(HeadType)
instances(DriveStyle)
instances(Material)

# build randdom inventory
function make_screw(prodnum)
    headtype = rand(instances(HeadType))
    drivestyle = rand(instances(DriveStyle))
    material = rand(instances(Material))

    Screw(prodnum, headtype, drivestyle, material)
end

screws = map(make_screw, 100:150)

# search functions with multiple dispatch and closures
function isround(screw)
    screw.headtype == rounded
end

function isround(screw, material)
    screw.headtype == rounded  &&
    screw.material == material
end

function isnonplastic(screw)
    screw.material != plastic
end

function issteel(screw)
    screw.material == steel
end

function ishex(screw)
    screw.drivestyle == hex
end

function isnonplastic(screw, headtype)
    screw.material != plastic &&
    screw.headtype == headtype
end

# closures - see below how to use inside filter instead of pre-defining
isroundaluminium = screw -> isround(screw, aluminium)
isroundsteel = screw -> isround(screw, steel)
isroundbrass = screw -> isround(screw, brass)
isroundplastic = screw -> isround(screw, plastic)
isroundwood = screw -> isround(screw, wood)

isnonplasticrounded = screw -> isnonplastic(screw, rounded)
isnonplasticflat = screw -> isnonplastic(screw, flat)
isnonplasticheadless = screw -> isnonplastic(screw, headless)
isnonplastictslot = screw -> isnonplastic(screw, tslot)

# queries - with closure either predefined or inside the filter
#filter(isround, screws)
#filter(isroundaluminium, screws)
#filter(screw -> isround(screw, aluminium), screws)
#filter(isroundsteel, screws)
#filter(isroundbrass, screws)
#filter(isroundplastic, screws)
#filter(isroundwood, screws)
#filter(isnonplastic, screws)
#filter(isnonplasticrounded, screws)
#filter(screw -> isnonplastic(screw, rounded), screws)
#filter(isnonplasticflat, screws)
#filter(isnonplasticheadless, screws)
#filter(isnonplastictslot, screws)

# queries with sets
#steel_screws = Set(filter(issteel, screws))
#hex_screws = Set(filter(ishex, screws))

# Dictionary solution
screwdict = Dict(screw.prodnum => screw for screw in screws)
#screwdict[137]

prodnums = keys(screwdict)

function isbrass(prodnum)
    screw = screwdict[prodnum]
    screw.material == brass
end
brass_screws = Set(filter(isbrass, prodnums))

# alternate function to use with closure
function ismaterial(prodnum, material)
    screw = screwdict[prodnum]
    screw.material == material
end

function istorx(prodnum)
    screw = screwdict[prodnum]
    screw.drivestyle == torx
end

function isdrivestyle(prodnum, drivestyle)
    screw = screwdict[prodnum]
    screw.drivestyle == drivestyle
end

brass_screws = Set(filter(isbrass, prodnums))
steel_screws = Set(filter(prodnum ->ismaterial(prodnum,steel), prodnums))
torx_screws = Set(filter(istorx, prodnums))
hex_screws = Set(filter(prodnum ->isdrivestyle(prodnum,hex), prodnums))

# Set queries
#brass_screws ∩ torx_screws
#steel_screws ∩ hex_screws
#brass_screws ∩ hex_screws
#steel_screws ∩ torx_screws
#[screwdict[pn] for pn in steel_screws ∩ torx_screws]

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
@testset "Testing discrete_rv.jl" begin

    # set up
    n = 10
    x = rand(n)
    x ./= sum(x)
    drv = DiscreteRV(x)

    # test Q sums to 1
    @test drv.Q[end] ≈ 1.0

    # test lln
    draws = rand(drv, 100_000)
    c = counter(draws)
    counts = Vector{Float64}(undef, n)
    for i=1:n
        counts[i] = c[i]
    end
    counts ./= sum(counts)

    @test isapprox(maximum(abs, counts - drv.q), 0.0; atol=1e-2)

    draws = Vector{Int}(undef, 100_000)
    rand!(draws, drv)
    c = counter(draws)
    counts = Vector{Float64}(undef, n)
    for i=1:n
        counts[i] = c[i]
    end
    counts ./= sum(counts)

    @test isapprox(maximum(abs, counts - drv.q), 0.0; atol=1e-2)

end  # testset

ext: jl


text:
# ------------------
# # Juliaを使ったParticle Swarm Optimization(PSO)アルゴリズムの実装
#
# - Rosenbrock functionの最適化
# - 多目的最適化にも拡張可能なように解および問題を定義
# ------------------

# Definition of a continuous optimization problem-------------------------
# 連続値最適化問題----------------------------------------------------------
"""
Abstract type of optimization problem.
最適化問題の抽象型
"""
abstract type AbstractProblem{T <: Number} end

"""
Concrete type of continuous optimization problem .
連続値最適化問題の具体型

# Fields
`objfunc`: julia function object of the objective function, 目的関数のオブジェクト
`numobjs`: number of objectives, 目的数
`numvars`: number of variables, 変数の数
`confunc`: julia function object of the constraint function, 制約計算する関数のオブジェクト
`numcons`: number of constraints, 制約の数
`limit`: upper and lower limit of the objective values, 目的関数の上下限値
`bound`: upper and lower bound of the variable values, 変数の上下限値
"""
mutable struct DoubleProblem{T} <: AbstractProblem{T}
    objfunc
    numobjs
    numvars
    confunc
    numcons
    limit::Matrix{Float64}
    bound::Matrix{T}
end
function DoubleProblem(objfunc; numvars=1, numobjs=1, confunc=NaN, numcons=0, limit=fill(NaN, 1, 1), bound=repeat([0.0 1.0], numvars))
    DoubleProblem(objfunc, numobjs, numvars, confunc, numcons, limit, bound)
end

"""
    norm(value, min, max)

Normalize value to 0-1.
`value`の値を0-1の範囲に正規化する
"""
norm(value, min, max) = (value - min) ./ (max - min)

"""
    denorm(value, min, max)

Denormalize value from 0-1 to min-max.
0-1に正規化された`value`の値をmin-maxの範囲に非正規化する
"""
denorm(value, min, max) = value .* (max - min) + min

"""
    evaluate!(problem, solutions)

Evaluate the objective value of each solution.
各解候補を評価し目的関数値を計算する
"""
function evaluate!(problem::DoubleProblem, solutions::Vector{<:AbstractSolution})
    for s in solutions
        vars = denorm(s.variables, problem.bound[:,1], problem.bound[:,2])
        objs = problem.objfunc(vars)
        if (isnan(problem.limit[1,1]))
            s.objectives = objs
        else
            s.objectives = norm(objs, problem.limit[:,1], problem.limit[:,2])
        end
    end
end

# Algorithm definition of the Particle Swarm Optimization---------------------
# Particle Swarm Optimizationの解とアルゴリズム定義-----------------------------
"""
Abstract type of solution.
最適化問題の解の抽象型
"""
abstract type AbstractSolution{T <: Number} end
"""
Concrete type of particle in PSO algorithm.
粒子群最適化における粒子(解)の具体型
"""
mutable struct Particle{T} <: AbstractSolution{T}
    objectives::Vector{Float64}
    variables::Vector{T}
    velocities::Vector{T}
    pbestvars::Vector{T}
    pbestobjs::Vector{Float64}
    constraints::Vector{T}
end

"""
    Particle(numvars[, numobjs, numcons]) -> Particle

Outer constractor of `Particle`.
`Particle`型の外部コンストラクタ

# Arguments 引数
- `numvars`: dimention (number of variables), 次元(変数の数)
- `numobjs = 1`: number of objectives, 目的数
- `numcons = 0`: number of constraints, 制約数
"""
function Particle(numvars, numobjs=1, numcons=0)
    objectives = ones(numobjs) .* Inf
    variables  = rand(numvars)
    velocities = zeros(numvars)
    pbestvars  = copy(variables)
    pbestobjs = ones(numobjs) .* Inf
    constraints = zeros(numcons)
    return Particle(objectives, variables, velocities, pbestvars, pbestobjs, constraints)
end

"""
    updateswarm!(swarm, gbest)

Update the position and velocity of the particles using the global best.
その世代の粒子群のうち最も良い解(gbest)を用いて，粒子の位置と速度を更新する
"""
function updateswarm!(swarm::Vector{Particle{Float64}}, gbest::Particle, c1=NaN, c2=NaN, w=NaN, vecr=false)
    for p in swarm
        if(vecr)
            r1 = rand(length(p.variables))
            r2 = rand(length(p.variables))
        else
            r1 = rand()
            r2 = rand()
        end
        if(isnan(c1)) c1 = rand() * 0.5 + 1.5 end
        if(isnan(c2)) c2 = rand() * 0.5 + 1.5 end
        if(isnan(w))  w = rand() * 0.4 + 0.1 end

        # 位置と速度の更新
        p.velocities = w * p.velocities + c1 * r1 * (p.pbestvars - p.variables) + c2 * r2 * (gbest.variables - p.variables)
        p.variables += p.velocities

        # 境界チェック
        for i = 1:length(p.variables)
            if (p.variables[i] > 1.0)
                p.variables[i] = 1.0
                p.velocities[i] *= -1.0
            elseif (p.variables[i] < 0.0)
                p.variables[i] = 0.0
                p.velocities[i] *= -1.0
            end
        end
    end
end

"""
    perturbation!(swarm)

Apply a mutation to each particles.
粒子位置に突然変異を適用する
"""
function perturbation!(swarm::Vector{Particle{Float64}})

end

"""
    updatepbest!(swarm)

Update the personal best of each particles.
その粒子がこれまで通った最も良い位置(pbest)を更新する
"""
function updatepbest!(swarm::Vector{Particle{Float64}})
    for p in swarm
        if (p.objectives < p.pbestobjs)
            p.pbestvars = copy(p.variables)
            p.pbestobjs = copy(p.objectives)
        end
    end
end

"""
    updategbest!(swarm, gbest)

Update the global best particle.
その世代の粒子群のうち最も良い解(gbest)を更新する
"""
function updategbest!(swarm::Vector{Particle{Float64}}, gbest::Particle{Float64})
    for p in swarm
        if (p.objectives < gbest.objectives)
            gbest.variables = copy(p.variables)
            gbest.objectives = copy(p.objectives)
        end
    end
end

"""
    pso(problem, iter) -> Particle

Algorithm of the Particle Swarm Optimization.
粒子群最適化アルゴリズムを実行しgbestを返す

# Arguments 引数
- `problem`: 問題の型(DoubleProblem型)
- `iter`: iterations, 繰り返し回数
- `nump`: number of particles, 粒子群サイズ．指定しなければ変数の数×10
- `c1`,`c2`,`w`: the parameters of PSO.
- `vecr`: use vectorized `r1`, `r2`. `r1`, `r2`をベクトル化するか否か(t/f)
"""
function pso(problem::DoubleProblem; iter, nump=problem.numvars*10, c1=NaN, c2=NaN, w=NaN, vecr=false)
    # 粒子群の初期化
    swarm = [Particle(problem.numvars) for i = 1:nump]
    gbest = Particle(problem.numvars)
    evaluate!(problem, swarm)
    updategbest!(swarm, gbest)
    updatepbest!(swarm)

    # 最適化の実行
    while (iter > 0)
        updateswarm!(swarm, gbest, c1, c2, w, vecr)
        perturbation!(swarm)
        evaluate!(problem, swarm)
        updategbest!(swarm, gbest)
        updatepbest!(swarm)
        iter -= 1
    end
    return gbest
end

"""
    rosenbrock(x) -> Float64

Function formulation of the Rosenbrock function.
See also: [Rosenbrock function - Wikipedia](https://en.wikipedia.org/wiki/Rosenbrock_function)
ローゼンブロック関数の定義
"""
function rosenbrock(x)
    [sum([ 100(x[i + 1] - x[i]^2)^2 + (x[i] - 1)^2 for i in 1:length(x) - 1])]
end

"""
Main function.
"""
function main()
    # 最適化問題の定義
    numvars = 3
    bound = [[-5.0 5.0];[-5.0 5.0];[-5.0 5.0]]
    problem = DoubleProblem(rosenbrock, numvars=numvars, bound=bound)

    # 最適化の実行
    res = pso(problem, iter=1000, nump=100)

    # 結果の表示
    vars = denorm(res.variables, bound[:,1], bound[:,2])
    objs = res.objectives
    println("objectives: " * string(objs))
    println("variables: " * string(vars))
end

@time main()

max_stars_repo_name: ohtayo/BlogSamples

max_issues_repo_name: ohtayo/BlogSamples

max_forks_repo_name: ohtayo/BlogSamples


text:
## ROCAnalysis.jl  Types and functions for performing ROC analysis
## (c) 2013--2016 David A. van Leeuwen.
##
## Licensed under the MIT software license, see LICENSE.md

module ROCAnalysis

export Roc, TNT, DCF, DataFrame, roc, chllr, eerch, pav, eer, plot, detplot, llrplot, apeplot, nbeplot, qnorm, pnorm, cllr, mincllr, auc, AUC, oeff, peff, plo, ber, minber, dcf, mindcf, setdcf, getdcf

include("modules.jl")

include("roctypes.jl")

include("includes.jl")

end


text:
include("yeast_dyn.jl")
include("shuffle.jl")
include("coop.jl")
include("motif.jl")
srand(1)

scale = 1
n = 11*scale
s = selfact;
for i=2:scale s = [s;selfact] end
selfact = s
n = length(selfact)
ne = 29*scale
np = 15*scale
interc = randomnkp(selfact,ne,np)
ntotal = 100000
nalpha = 16
amin = 0.5
alpha= amin + (1-amin)*[0:nalpha]/nalpha - 0.5*(1-amin)/nalpha # linspace(0,1,n)?
palpha = zeros(Float64,1+nalpha)
nmotif2=3^2
#nmotif3=3^6
nmotif3=138
ninit=2^11
listc = Array(Float64,0)
motifcoop = zeros(Float64,nmotif3,1+nalpha)
#interc = randomnkp(selfact,ne,np)
i=0
while (i<ntotal)
    (intersh,selfsh) = (switchnet(interc),selfact)
#    (intersh,selfsh) = shufflenet(interc,selfact)
    if (isconnected(symadj(intersh)))
#        (a,actratio) = cooppointattr(intersh,selfsh)
        (a,actratio) = coopattr(intersh,selfsh,ninit)
        push!(listc,a)
        if (a >= amin)      # returns -1 if no nontrivial attractor found
            na = 1+ifloor(nalpha*(a-amin)/(1-amin) + 0.5)
            palpha[na] += 1
            motifcoop[:,na] += motifcount3(intersh)
            i += 1
            if i%(ntotal/100)==0 print(".") end
        end
    end
end

palpha /= ntotal
writedlm("alphas.$n.$ne.$np.$ntotal.dat",listc,' ')
writedlm("motifcoop.$n.$ne.$np.$ntotal.dat",motifcoop',' ')
writedlm("palpha.$n.$ne.$np.$ntotal.dat",palpha,' ')
motifcoopnorm = similar(motifcoop)
for i=1:nmotif3
    motifcoopnorm[i,:] = motifcoop[i,:] ./ palpha'
end
writedlm("motifcoop_norm.$n.$ne.$np.$ntotal.dat",motifcoopnorm',' ')

## for i=1:nmotif3
##     readline()
##     plot(alpha,motifcoopnorm[i,:])
## end

max_stars_repo_name: akabakcioglu/network-dynamics

max_issues_repo_name: akabakcioglu/network-dynamics

max_forks_repo_name: akabakcioglu/network-dynamics


text:
function grpdelay(b,a=1,nfft=512,whole="",Fs=0)
  if whole != "whole"; nfft = 2*nfft; end

  # this splatting behavior ensures a becomes a vector,
  # even if scalar is passed in
  b = [b...]; # 'splatting', unfortunate but correct as of 0.6
  a = [a...]; # 'splatting', unfortunate but correct as of 0.6
  oa = length(a)-1;         # order of a(z)
  oc = oa + length(b)-1;    # order of c(z)
  c = conv(b,flipdim(a,1)); # c(z) = b(z)*a(1/z)*z^(-oa)
  cr = c .* collect(0:oc);  # derivative of c wrt 1/z

  # Update nfft if necessary, after looking at a and b
  if length(cr)>nfft; nfft = nextpow2(length(cr)); end;
  w = 2*pi*collect(0:nfft-1)/nfft;
  if Fs>0; w = Fs*w/(2*pi); end

  cr = cat(1,cr,zeros(nfft-length(cr))); # zero pad for fft
  c  = cat(1,c ,zeros(nfft-length(c ))); # zero pad for fft
  num = fft(cr,1); # nfft not supported :-\
  den = fft(c ,1); # nfft not supported :-\
  minmag = 10*eps();
  polebins = find(abs(den) .< minmag);
  for b=polebins
    println("*** grpdelay: group delay singular! setting to 0");
    num[b] = 0;
    den[b] = 1;
  end
  gd = real(num ./ den) - oa;

  if whole != "whole"
    ns = Int(nfft/2); # Matlab convention - should be nfft/2 + 1
    gd = gd[1:ns];
    w = w[1:ns];
  end

  w = reshape(vec(w),(length(w),1)); # Matlab returns column vectors
  gd = reshape(vec(gd),(length(gd),1)); # Matlab returns column vectors

  return gd,w
end


text:
# This file can be used to create the validation file for deterministic scc values.
# It will create a validation file containing all possibilities of parameter values
# defined in the specs dictionary below produce the same results.

using MimiDICE2013
using DataFrames
using Query
using CSVFiles
using Test

specs = Dict([
    :year => [2020, 2050],
    :eta => [0, 1.5],
    :prtp => [0.015, 0.03],
    :last_year => [2200, 2305],
])

results = DataFrame(year = [], eta = [], prtp = [], last_year = [], SC = [])

for year in specs[:year]
    for eta in specs[:eta]
        for prtp in specs[:prtp]
            for last_year in specs[:last_year]
                sc = MimiDICE2013.compute_scc(year=Int(year), eta=eta, prtp=prtp, last_year=Int(last_year))
                push!(results, (year, eta, prtp, last_year, sc))
            end
        end
    end
end

path = joinpath(@__DIR__, "deterministic_sc_values_v1-0-2.csv")
save(path, results)


text:
using VAWTools
using Test, Random
using SHA

function sha(fn)
    open(fn) do f
        sha2_256(f)
    end
end

include("smoothing-functions.jl")
include("elevation-bands.jl")

## Gridded
##########
g = Gridded(1.:37, 1:17., rand(37,17))
gc = Gridded(1.:10,1:11., 5*ones(10,11) )

@test g.v==downsample(g,1).v
@test_throws ArgumentError downsample(g,0)

for i=1:9
    @test all(downsample(gc,i).v.==5)
    @test downsample(g,i).v==boxcar(g.v,Int(floor((i)/2)))[1:i:end,1:i:end]
end
@test_throws AssertionError downsample(gc,10)

gc.v[4,5] = NaN
for i=1:9
    ds = downsample(gc,i).v
    @test all((ds.==5) .| (isnan.(ds)))
    @test downsample(g,i).v==boxcar(g.v,Int(floor((i)/2)))[1:i:end,1:i:end]
end
averagemask = trues(size(gc.v))
for i=1:9
    ds = downsample(gc,i,1,true,averagemask).v
    @test all((ds.==5) .| (isnan.(ds)))
    @test downsample(g,i,1,true,trues(size(g.v))).v==boxcar(g.v,Int(floor((i)/2)))[1:i:end,1:i:end]
end
averagemask = trues(size(g.v))
averagemask[2,1] = false
for i=1:9
    ds = downsample(g,i,1,true,averagemask).v
    bc = boxcar(g.v,Int(floor((i)/2)),averagemask)[1:i:end,1:i:end]
    @test bc==ds
end

## File IO
#################
const todelete = []
function tempfn(ext="")
    out = tempname()*ext
    push!(todelete, out)
    out
end
function delfls()
    for f in todelete
        if isfile(f)
            rm(f)
        end
    end
end

## .agr readers:
g1 = VAWTools._read_agr("testfiles/wiki.agr")
g1b = VAWTools._read_agr("testfiles/wiki.bin")
@test g1.v == [-9999 -9999 5 2
                -9999 20 100 36
                3 8 35 10
                32 42 50 6
                88 75 27 9
                13 5 1 -9999]
@test g1.nc==4
@test g1.nr==6
@test g1.xll==0.0
@test g1.yll==-40.0
@test g1.dx==50.0
@test g1.NA==-9999
# write and read back:
tfn1 = tempfn(".agr")
write_agr(g1, tfn1)
g2 = VAWTools._read_agr(tfn1, Float64)
@test g1==g2
tfn2 = tempfn(".agr")
write_agr(g2, tfn2)
@test sha(tfn1)==sha(tfn2)

# reading from io thingies
io = open("testfiles/wiki.agr", "r")
g2 = VAWTools.read_agr(io)
close(io)
g1 = VAWTools.read_agr("testfiles/wiki.agr")
@test g1.x==g2.x
@test g1.y==g2.y
@test g1.midpoint==g2.midpoint

# Matthias abuse of NODATA_value field as UTM zone:
VAWTools.read_agr("testfiles/wiki.agr")

## RasterIO agr reading
g1 = VAWTools.read_agr("testfiles/wiki.agr")
# println("Expected output: ERROR 1: ERROR - failed to load SRS ...")
# g2,proj4 = VAWTools.read_rasterio("testfiles/wiki.agr")
# @test isequal(g1.v,g2.v)
# @test g1.x==g2.x
# @test g1.y==g2.y
# @test g1.midpoint==g2.midpoint
# @test proj4==""

# write and read back binary:
tfn3 = tempfn(".bin")
write_agr(g1, tfn3, NA_agr=-9999.0)
g2 = VAWTools._read_agr(tfn3)
g1_ = VAWTools._read_agr("testfiles/wiki.agr", Float32)
@test g1_==g2
@test sha(tfn3)==sha("testfiles/wiki.bin")

# larger file:
g3 = VAWTools._read_agr("testfiles/t2.bin", Float32)
@test g3.nc==121
@test g3.nr==81
@test g3.xll==678750.0f0
@test g3.yll==140250.0f0
@test g3.dx==25.0f0
@test g3.NA==-9999.0f0

tfn4 = tempfn(".bin")
write_agr(g3, tfn4)
@test sha("testfiles/t2.bin")==sha(tfn4)

# Gridded
g1 = VAWTools._read_agr("testfiles/wiki.agr")
gg = Gridded(g1)
@test isequal(VAWTools.AGR(gg, NA_agr=-9999).v, g1.v)

# geotiff
gt = VAWTools.read_geotiff("testfiles/small_world.tif")
gt_ag = VAWTools.read_geotiff("testfiles/small_world.asc")
@test gt.proj=="+proj=longlat +datum=WGS84 +no_defs"
@test gt.x==gt_ag.x
@test gt.y==gt_ag.y
@test isa(gt.v, Array{Float32,2})
@test gt.v==gt_ag.v

# Trajectory
@test_throws AssertionError Traj(1:5, 6:11, [0.0:5.0;])

tr = Traj(1.0:5, 6.0:10)
@test tr.x==collect(1.0:5)
@test tr.y==collect(6.0:10)
@test tr.splits==[1:5]
@test !VAWTools.haserror(tr)
@test !VAWTools.hasvalues(tr)

tr = Traj(1.0:5, 6.0:10, Nothing[], Nothing[], [1:3,4:5])
@test tr.x==collect(1.0:5)
@test tr.y==collect(6.0:10)
@test tr.splits==[1:3,4:5]
@test !VAWTools.haserror(tr)
@test !VAWTools.hasvalues(tr)

tr = Traj(1.0:5, 6.0:10, 1.0:5.0)
@test tr.x==collect(1.0:5)
@test tr.y==collect(6.0:10)
@test tr.v==collect(1.0:5)
@test tr.splits==[1:5]
@test !VAWTools.haserror(tr)
@test VAWTools.hasvalues(tr)

tr = Traj(1.0:5, 6.0:10, 1.0:5.0, Float64[], [1:3,4:5])
@test tr.x==collect(1.0:5)
@test tr.y==collect(6.0:10)
@test tr.v==collect(1.0:5)
@test tr.splits==[1:3,4:5]
@test !VAWTools.haserror(tr)
@test VAWTools.hasvalues(tr)

tr = Traj(1.0:5, 6.0:10, 1.0:5.0, 1.0:5.0)
@test tr.x==collect(1.0:5)
@test tr.y==collect(6.0:10)
@test tr.v==collect(1.0:5)
@test tr.err==collect(1.0:5)
@test tr.splits==[1:5]
@test VAWTools.haserror(tr)
@test VAWTools.hasvalues(tr)

tr = Traj(1.0:5, 6.0:10, 1.0:5.0, 1.0:5.0, [1:3,4:5])
@test tr.x==collect(1.0:5)
@test tr.y==collect(6.0:10)
@test tr.v==collect(1.0:5)
@test tr.err==collect(1.0:5)
@test tr.splits==[1:3,4:5]
@test VAWTools.haserror(tr)
@test VAWTools.hasvalues(tr)
@test all(VAWTools.distances(tr).==sqrt(2))

# split it
tr = Traj([1,2,4,5], 1:4)
@test tr.splits==[1:4]
@test !VAWTools.haserror(tr)
@test !VAWTools.hasvalues(tr)
VAWTools.split_traj!(tr,2)
@test tr.splits==[1:2,3:4]
@test !VAWTools.haserror(tr)
@test !VAWTools.hasvalues(tr)

## .xyn reader
po = read_xyn("testfiles/poly.xyn")
@test length(po)==1
@test size(po[1])==(2,24)
cpo = VAWTools.concat_poly(po)
@test po[1]==cpo[1]
@test po==VAWTools.split_poly(cpo...)

po = read_xyn("testfiles/multipoly.xyzn", hasz=true)
@test length(po)==4
@test size(po[1])==(3,4)
@test size(po[2])==(3,5)
@test size(po[3])==(3,4)
@test size(po[4])==(3,3)
cpo = VAWTools.concat_poly(po)
@test po==VAWTools.split_poly(cpo...)

# tidy up temp-files
delfls()

## IN poly
@test VAWTools.leftorright(0.5,0.5, 1,0,1,1)==-1
@test VAWTools.leftorright(1.5,.5, 1,0,1,1)==1
@test VAWTools.leftorright(1,0.5, 1,0,1,1)==0

poc = VAWTools.concat_poly(po)
@test !inpoly([0.0,0], poc[1][1:2,:])
@test inpoly([6.28793e5, 90129.7+10], poc[1][1:2,:])

# counterclockwise poly
poly = Float64[0  0  1  1  0
               0  1  1  0  0]
# plot with:
# using Plots; plot(poly[1,:], poly[2,:])

p1 = [0.5, 0.5]
p2 = [0.5, 0.99]
p22 = [0.5, 1] # on top edge (per def outside)
p22_ = [1, 0.5] # on right edge (per def outside)
p23 = [0.5, 0] # on bottom edge (per def inside)
p23_ = [0, 0.5] # on left edge (per def inside)
p24 = [0, 0]   # on corner
p25 = [0, .4]   # on edge
p3 = [0.5, 1.1]

@test inpoly(p1, poly)
@test inpoly(p2, poly)
@test !inpoly(p22, poly)
@test !inpoly(p22_, poly)
@test inpoly(p23, poly)
@test inpoly(p23_, poly)
@test inpoly(p24, poly)
@test inpoly(p25, poly)
@test !inpoly(p3, poly)

# clockwise poly
poly = Float64[0  1  1  0  0
               0  0  1  1  0]

@test inpoly(p1, poly)
@test inpoly(p2, poly)
@test !inpoly(p22, poly)
@test !inpoly(p22_, poly)
@test inpoly(p23, poly)
@test inpoly(p23_, poly)
@test inpoly(p24, poly)
@test inpoly(p25, poly)
@test !inpoly(p3, poly)

# cross-over poly
poly = Float64[ 0  1  0  1  0
                0  0  1  1  0]

eval(:(@test_broken inpoly(p1, poly) )) # should be true
@test inpoly(p2, poly)
@test !inpoly(p22, poly)
@test !inpoly(p22_, poly)
@test inpoly(p23, poly)
@test !inpoly(p23_, poly)
@test inpoly(p24, poly)
@test !inpoly(p25, poly) # different
@test !inpoly(p3, poly)

# with interior region
poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside interior poly: i.e. labeled as outside
@test !inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               # in-interior
               0.4 0.4
               0.4 0.2
               0.2 0.2
               0.2 0.4
               0.4 0.4
               # interior
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside in-interior poly
@test inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               # in-interior
               0.4 0.4
               0.2 0.4
               0.2 0.2
               0.4 0.2
               0.4 0.4
               # interior
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside in-interior poly
@test inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior #1
               0.1 0.1
               0.1 0.6
               0.4 0.6
               0.4 0.6
               0.4 0.1
               0.1 0.1
               0 0
               # interior #2
               0.6 0.4
               0.6 0.6
               0.8 0.6
               0.8 0.4
               0.6 0.4
               0 0
               # exterior
               0 1
               1 1
               1 0
               0 0]'
@test !inpoly([0.2,0.4], poly)
@test !inpoly([0.3,0.15], poly)
@test inpoly([0.5,0.4], poly)
@test inpoly([0.5,0.2], poly)
@test !inpoly([0.7,0.5], poly)

# Test points which are on same y-level as a endpoint of poly edge
poly = Float64[0.1 1 2 1  0.1
               0.1 -1 0 1 0.1];
@test inpoly([1,0.0], poly)
@test inpoly([1,-0.0001], poly)
@test inpoly([1,0.0001], poly)
poly = Float64[0.0 1 2 1  0.0
               0.0 -1 0 1 0.0];
@test inpoly([1,0.0], poly)
@test inpoly([1,-0.0001], poly)
@test inpoly([1,0.0001], poly)

#####
# boxcar
#####
## 1D
orig = [1:10.0;]
window = 3
@inferred VAWTools.boxcar(orig, window)
filt1 = VAWTools.boxcar(orig, window)
@test size(filt1)==size(orig)
@test filt1==[2.5, 3.0, 3.5, 4.0, 5.0, 6.0, 7.0, 7.5, 8.0, 8.5]

window = 3.3
@inferred VAWTools.boxcar(orig, window)
filt1 = VAWTools.boxcar(orig, window)
@test size(filt1)==size(orig)
@test filt1≈[2.67442, 3.16981, 3.66667, 4.16438, 5.0, 6.0, 6.83562, 7.33333, 7.83019, 8.32558] atol=1e-4

windows = (2.4,3.3)
@inferred VAWTools.boxcar(orig, windows)
filt1 = VAWTools.boxcar(orig, windows)
@test size(filt1)==size(orig)
@test filt1≈[2.67442, 3.16981, 3.66667, 4.44776, 5.44776, 6.44776, 7.28125, 7.77778, 8.27273, 8.76471] atol=1e-4

## 2D
#nr,nc = 4,5
T = Float32
nr,nc = 20,31
#nr,nc = 600,500
window = 3
# orig = (1.0:nr)''*(1:nc)'
# weights = ones(Int,nr,nc)
Random.seed!(1)
orig = rand(T,nr,nc)
weights = rand(T,nr,nc)
weightsb = bitrand(nr,nc)
weightsbb = convert(Matrix{Bool}, weightsb)

dropmask = falses(nr,nc)
dropmask[5,6] = true
# with weights
filt1 = VAWTools.boxcar(orig, window, weights)
@inferred VAWTools.boxcar(orig, window, weights)
M = VAWTools.boxcar_matrix(T, window, weights)
@inferred VAWTools.boxcar_matrix(T, window, weights)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = VAWTools.apply_boxcar_matrix(M, orig)
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test filt1[i] ≈ filt2[i]
end
# with weightsb
filt1 = VAWTools.boxcar(orig, window, weightsb)
@inferred VAWTools.boxcar(orig, window, weightsb)
M = VAWTools.boxcar_matrix(T, window, weightsb)
@inferred VAWTools.boxcar_matrix(T, window, weightsb)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = VAWTools.apply_boxcar_matrix(M, orig)
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test filt1[i] ≈ filt2[i]
end
# with weightsbb
filt1 = VAWTools.boxcar(orig, window, weightsbb)
@inferred VAWTools.boxcar(orig, window, weightsbb)
M = VAWTools.boxcar_matrix(T, window, weightsbb)
@inferred VAWTools.boxcar_matrix(T, window, weightsbb)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = VAWTools.apply_boxcar_matrix(M, orig)
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test filt1[i] ≈ filt2[i]
end

# with weights and dropmask
filt1 = VAWTools.boxcar(orig, window, weightsbb, dropmask)
@inferred VAWTools.boxcar(orig, window, weightsbb, dropmask)
M = VAWTools.boxcar_matrix(T, window, weightsbb, dropmask)
@inferred VAWTools.boxcar_matrix(T, window, weightsbb, dropmask)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = VAWTools.apply_boxcar_matrix(M, orig)
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test filt1[i] ≈ filt2[i]
end

# with NaN poisoning
orig[1,1] = NaN
weightsbb[1,1] = 1
filt1 = VAWTools.boxcar(orig, window, weightsbb, dropmask)
M = VAWTools.boxcar_matrix(T, window, weightsbb, dropmask)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = VAWTools.apply_boxcar_matrix(M, orig)
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    if isnan(filt1[i])
        @test isequal(filt1[i], filt2[i])
    else
        @test filt1[i] ≈ filt2[i]
    end
end

# NaN is masked: no poisoning
orig[1,1] = NaN
weightsbb[1,1] = 0
filt1 = VAWTools.boxcar(orig, window, weightsbb, dropmask)
M = VAWTools.boxcar_matrix(T, window, weightsbb, dropmask)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = VAWTools.apply_boxcar_matrix(M, orig)
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    if isnan(filt1[i])
        @test isequal(filt1[i], filt2[i])
    else
        @test filt1[i] ≈ filt2[i]
    end
end

## different left and right windows
windows1 = (3,4)
windows2 = ((2,3), (1,2))

Random.seed!(1)
orig = rand(T,nr,nc)

for w in [windows1, windows2]
    @inferred VAWTools.boxcar(orig, w)
    @inferred VAWTools.boxcar(orig, w, weights)
    @inferred VAWTools.boxcar(orig, w, weights, dropmask)
end

@inferred VAWTools.boxcar(orig, (3.3,4.2))

####
# piecewiselinear
####
f = VAWTools.piecewiselinear([78.0], [-12.0])
@test f(-1)==-12
@test f(100)==-12

f = VAWTools.piecewiselinear([0, 1, 2], [-1, 1, -2])
@test f(0)==-1
@test f(0.5)==0
@test f(1)==1
@test f(1.5)==-0.5
@test f(2)==-2

#####
# projections
####
import Proj4
wgs84 = Proj4.Projection("+proj=longlat +datum=WGS84 +no_defs")
utm56 = Proj4.Projection("+proj=utm +zone=56 +south +datum=WGS84 +units=m +no_defs")

pt  = [150.0 -27 110]
@test isapprox(VAWTools.transform_proj(VAWTools.transform_proj(pt, wgs84, utm56), utm56, wgs84), pt)


text:
# =========================================================================== #

# =========================================================================== #

#Implementing an x0 viable solution huristic

function Glouton(cost, matrix, n, m)

    #Creating a set of lines that will be avaluated
    actif= zeros(Bool, m) #if actif[j]=0 the line will be avaluated
    stop= ones(Bool,m)

    #creating my utility array for my m variables
    util = zeros(Float64, n)

    #initialisation de la Solution
    SOL = zeros(Float64, n)

    println( "\n Matrix: ", matrix, "\n Cost: ", cost,"\n n= ", n, "   m= ", m,"\n actif: ", actif, "\n stop: ", stop, "\n util: ", util, "\n SOL: ", SOL, "\n \n")

    while actif!=stop
        #Foction d'utilité
        util = Utilite(cost, matrix, actif, n, m)

        println("\n **********************************\n utilité finiiiiiii! :D \n *************************************\n")

        #Choix parmis les candidats de util
        PosCandidat = PosMax(util,n)

        println("\n +++++++++++++ \n PosMax est fini! =) \n Le candidat à la position: ", PosCandidat, "\n Type de de candidat: ", typeof(PosCandidat), "\n +++++++++++++++++ \n")

        #On ajoute le candidat à la solution
        SOL[PosCandidat] = 1

        #Desactive! le candidat selectionné
        Desactive!(PosCandidat, matrix, actif, m)
        println("\n Desactive! est fini avec actif: ", actif,"\n Et SOL: ", SOL,  "\n /////////////////////////////////////////////////")

    end

    return(SOL)

end

# =========================================================================== #

#Detremine an utility  based on an active (actif) set of matrix's lines
function Utilite(cost, matrix, actif, n, m)

    util = zeros(Float64, n)
    i=1
    #On each column of matrix
    for j=1:n
        println("\n \n ============================= \n j= ", j, "  i= ", i, "\n util: ", util, "\n actif[i]: ", actif[i])

        #For each line
        for i=1:m
            println(" actif=", actif, "\n util: ", util)

            #checking if the line is active
            if actif[i]==0

                println("\n i=", i,"  j= ", j, "\n util: ", util)

                K=matrix[i,j]

                println("\n k= ", K)

                util[j]=util[j]+K

                println("\n util: ", util)

            end
        end

        # dividing the number if iterations of the variable j by its cost (not a zero)

        println("\n ***** conversion: ***** \n ", " Util[j]= ", util[j], "  cost[j] = ", cost[j])

        if util[j]!=false

            util[j] =  cost[j] / util[j]

        println(" util= ", util)
        end
    end
    return util
end

# =========================================================================== #

function PosMax(util,n) #On choisie un candidat parmie les différents coûts, par ordre lexicographique

    Pos=1
    ValCan=0.0

    println("\n +++++++++++++++++++\n Debut PosMax\n ", " Pos: ", Pos, " ValCan: ", ValCan, "\n util: ", util, "\n n: ", n, "\n" )

    for i=1:n

        println("\n boucle: ", i)
        if util[i] > ValCan

            ValCan = util[i]
            Pos = i
        end
    end

    println("\n PosMax est fini avec: \n Max: ", ValCan, " Position: ", Pos, "\n")
    return Pos
end

# =========================================================================== #

function Desactive!(PosCandidat, matrix, actif, m) #On désactive les lignes où est le candidat

    println("\n \n ////////////////////////////////////////// \n Debut de Desactive! \n Var check: \n position candidat: ", PosCandidat, "\n matrice : ", matrix, "\n actif: ", actif, "\n n: ", n, "\n")
    for i=1:m
        if matrix[i,PosCandidat] == 1
            actif[i] = 1

        end
    end
    println("Fin Desactive! ")
end


text:
# Beta distribution

export Beta

@parameterized Beta(α,β) ≃ Lebesgue(𝕀)

@kwstruct Beta(α, β)

@kwalias Beta [
    a     => α
    alpha => α
    b     => β
    beta  => β
]

TV.as(::Beta) = as𝕀

function logdensity(d::Beta{(:α, :β)}, x)
    return (d.α - 1) * log(x) + (d.β - 1) * log(1 - x) - logbeta(d.α, d.β)
end

Base.rand(rng::AbstractRNG, T::Type, μ::Beta) = rand(rng, Dists.Beta(μ.α, μ.β))

distproxy(d::Beta{(:α, :β)}) = Dists.Beta(d.α, d.β)

asparams(::Type{<:Beta}, ::Val{:α}) = asℝ₊
asparams(::Type{<:Beta}, ::Val{:β}) = asℝ₊


text:
using .Makie: Makie

Makie.plottype(::AOI) = Makie.Image

function Makie.convert_arguments(P::Type{<:Makie.Image}, aoi::AOI)
  Makie.convert_arguments(P,
    range(aoi.topleft.lon, aoi.bottomright.lon; length=size(aoi.mapimg,2)),
    range(aoi.topleft.lat, aoi.bottomright.lat; length=size(aoi.mapimg,1)),
    rotr90(aoi.mapimg))
end

Makie.convert_single_argument(x::LLA) = Makie.Point2f(x.lon, x.lat)
Makie.convert_single_argument(x::AbstractVector{<:LLA}) = [Makie.Point2f(v.lon, v.lat) for v ∈ x]


text:
# returns the interpolated value of f(x), where f is defined on the
# grid R
# input parameters
# R ... grid
# n ... length of the grid
# f ... function defined on the grid R
# x ... point at which to interpolate
# i ... the integer number of the radial point at which "x" lies
# output parameters
# V ... the interpolated function value
function radial_interp(r, f, x::Float64, i::Int64)
    #IMPLICIT NONE
    #INTEGER, INTENT(in) :: n
    #REAL(8), INTENT(in) :: f(n), x, R(n)
    #REAL(8), INTENT(out) :: V
    #INTEGER, INTENT(in) :: i

    N = size(r,1)
    if N < 4
        error("radial_interp: n >= 4 required")
    end

    j1 = i - 1
    j2 = j1 + 1

    n1 = j1 - 1
    n2 = j2 + 1

    if n1 < 1
        n2 = n2 - n1 + 1
        n1 = 1
    end

    if n2 > N
        n1 = n1 - (n2 - N)
        n2 = N
    end

    return @views lagrange_interp( r[n1:n2], f[n1:n2], x )
end


text:
const GWPos = SVector{2,Int}

"""
    SimpleGridWorld(;kwargs...)

Create a simple grid world MDP. Options are specified with keyword arguments.

# States and Actions
The states are represented by 2-element static vectors of integers. Typically any Julia `AbstractVector` e.g. `[x,y]` can also be used for arguments. Actions are the symbols `:up`, `:left`, `:down`, and `:right`.

# Keyword Arguments
- `size::Tuple{Int, Int}`: Number of cells in the x and y direction [default: `(10,10)`]
- `rewards::Dict`: Dictionary mapping cells to the reward in that cell, e.g. `Dict([1,2]=>10.0)`. Default reward for unlisted cells is 0.0
- `terminate_from::Set`: Set of cells from which the problem will terminate. Note that these states are not themselves terminal, but from these states, the next transition will be to a terminal state. [default: `Set(keys(rewards))`]
- `tprob::Float64`: Probability of a successful transition in the direction specified by the action. The remaining probability is divided between the other neighbors. [default: `0.7`]
- `discount::Float64`: Discount factor [default: `0.95`]
"""
@with_kw struct SimpleGridWorld <: MDP{GWPos, Symbol}
    size::Tuple{Int, Int}           = (10,10)
    rewards::Dict{GWPos, Float64}   = Dict(GWPos(4,3)=>-10.0, GWPos(4,6)=>-5.0, GWPos(9,3)=>10.0, GWPos(8,8)=>3.0)
    terminate_from::Set{GWPos}      = Set(keys(rewards))
    tprob::Float64                  = 0.7
    discount::Float64               = 0.95
end

# States

function POMDPs.states(mdp::SimpleGridWorld)
    ss = vec(GWPos[GWPos(x, y) for x in 1:mdp.size[1], y in 1:mdp.size[2]])
    push!(ss, GWPos(-1,-1))
    return ss
end

function POMDPs.stateindex(mdp::SimpleGridWorld, s::AbstractVector{Int})
    if all(s.>0)
        return LinearIndices(mdp.size)[s...]
    else
        return prod(mdp.size) + 1 # TODO: Change
    end
end

struct GWUniform
    size::Tuple{Int, Int}
end
Base.rand(rng::AbstractRNG, d::GWUniform) = GWPos(rand(rng, 1:d.size[1]), rand(rng, 1:d.size[2]))
function POMDPs.pdf(d::GWUniform, s::GWPos)
    if all(1 .<= s[1] .<= d.size)
        return 1/prod(d.size)
    else
        return 0.0
    end
end
POMDPs.support(d::GWUniform) = (GWPos(x, y) for x in 1:d.size[1], y in 1:d.size[2])

POMDPs.initialstate_distribution(mdp::SimpleGridWorld) = GWUniform(mdp.size)

# Actions

POMDPs.actions(mdp::SimpleGridWorld) = (:up, :down, :left, :right)
Base.rand(rng::AbstractRNG, t::NTuple{L,Symbol}) where L = t[rand(rng, 1:length(t))] # don't know why this doesn't work out of the box

const dir = Dict(:up=>GWPos(0,1), :down=>GWPos(0,-1), :left=>GWPos(-1,0), :right=>GWPos(1,0))
const aind = Dict(:up=>1, :down=>2, :left=>3, :right=>4)

POMDPs.actionindex(mdp::SimpleGridWorld, a::Symbol) = aind[a]

# Transitions

POMDPs.isterminal(m::SimpleGridWorld, s::AbstractVector{Int}) = any(s.<0)

function POMDPs.transition(mdp::SimpleGridWorld, s::AbstractVector{Int}, a::Symbol)
    if s in mdp.terminate_from || isterminal(mdp, s)
        return Deterministic(GWPos(-1,-1))
    end

    destinations = MVector{length(actions(mdp))+1, GWPos}(undef)
    destinations[1] = s

    probs = @MVector(zeros(length(actions(mdp))+1))
    for (i, act) in enumerate(actions(mdp))
        if act == a
            prob = mdp.tprob # probability of transitioning to the desired cell
        else
            prob = (1.0 - mdp.tprob)/(length(actions(mdp)) - 1) # probability of transitioning to another cell
        end

        dest = s + dir[act]
        destinations[i+1] = dest

        if !inbounds(mdp, dest) # hit an edge and come back
            probs[1] += prob
        else
            probs[i+1] += prob
        end
    end

    return SparseCat(destinations, probs)
end

function inbounds(m::SimpleGridWorld, s::AbstractVector{Int})
    return 1 <= s[1] <= m.size[1] && 1 <= s[2] <= m.size[2]
end

# Rewards

POMDPs.reward(mdp::SimpleGridWorld, s::AbstractVector{Int}) = get(mdp.rewards, s, 0.0)
POMDPs.reward(mdp::SimpleGridWorld, s::AbstractVector{Int}, a::Symbol) = reward(mdp, s)

# discount

POMDPs.discount(mdp::SimpleGridWorld) = mdp.discount

# Conversion
function POMDPs.convert_a(::Type{V}, a::Symbol, m::SimpleGridWorld) where {V<:AbstractArray}
    convert(V, [aind[a]])
end
function POMDPs.convert_a(::Type{Symbol}, vec::V, m::SimpleGridWorld) where {V<:AbstractArray}
    actions(m)[convert(Int, first(vec))]
end

# Fallback Render
function POMDPModelTools.render(m::SimpleGridWorld, step; kwargs...)
    return "$step\n\nPlease import Compose to enable SimpleGridWorld visualization."
end


text:
const FrontMatter = Dict{String,Any}

"""
    frontmatter(nb::Notebook; raise::Bool=false)::Dict{String,Any}

Extract frontmatter from a notebook, which is extra meta-information that the author attaches to the notebook, often including *title*, *description*, *tags*, *author*, and more. Search for *frontmatter* online to learn more.

If `raise` is true, then parsing errors will be rethrown. If `false`, this function will always return a `Dict`.

Currently, you can give frontmatter to a notebook by defining a global variable `frontmatter` anywhere in the notebook, which should be a named tuple. Like so:

```julia
frontmatter = (
	title = "⚡ JavaScript for widgets",
	description = "A simple notebook demonstrating the Julia syntax essentials",
	layout = "layout.jlhtml",
	tags = ["docs", "advanced", "widgets"],
);
```

This definition will be picked up (using syntax analysis) and evaluated. Note:
- You are not allowed to use variables defined in other cells. (But you can *use* `frontmatter` in your notebook.)
- You can use `;` to hide it from being displayed, only the definition matters.

!!! warning
    This will change in the future, because we might make this a GUI thing! Consider this function in the experimental phase.
"""
function frontmatter(nb::Notebook; raise::Bool=false)
    top = updated_topology(nb.topology, nb, nb.cells)

	cs = where_assigned(top, Set([:frontmatter]))
	if isempty(cs)
		FrontMatter()
	else
		try
			c = only(cs)
			m = Module()
			Core.eval(m, Meta.parse(c.code))
			result = Core.eval(m, :frontmatter)
			if result isa FrontMatter
				result
			else
				FrontMatter(String(k)=>v for (k,v) in pairs(result))
			end
		catch e
            if raise
                rethrow(e)
            else
                @error "Error reading frontmatter. Make sure that `frontmatter` is a `NamedTuple` or a `Dict{String,Any}`, and that it does not use global variables." nb.path exception=(e,catch_backtrace())
                FrontMatter()
            end
		end
	end
end

function frontmatter(abs_path::String)
	# this will load the notebook to analyse, it won't run it
    frontmatter(load_notebook_nobackup(abs_path))
end


text:
using HighlyAdaptiveLasso
using Documenter

DocMeta.setdocmeta!(HighlyAdaptiveLasso, :DocTestSetup, :(using HighlyAdaptiveLasso); recursive=true)

makedocs(;
    modules=[HighlyAdaptiveLasso],
    authors="Olivier Labayle",
    repo="https://github.com/olivierlabayle/HighlyAdaptiveLasso.jl/blob/{commit}{path}#{line}",
    sitename="HighlyAdaptiveLasso.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://olivierlabayle.github.io/HighlyAdaptiveLasso.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/olivierlabayle/HighlyAdaptiveLasso.jl",
)


text:
function MALA(U_fun::Function, grad_U::Function, h::Real, dim::Integer, n_iter::Integer, Xi_init = zeros(dim))

  Xi = zeros(dim, n_iter + 1)
  Xi[:,1] = Xi_init
  Z = randn(dim, n_iter)
  U = rand(n_iter)
  q(x, y) = exp(-1/(2*h) * sumabs2(y - x + (h/2)*grad_U(x)))
  acceptance_count = 0
  for i=1:n_iter
    Xi_current = Xi[:,i]
    Xi_proposed = Xi_current - (h/2) * grad_U(Xi_current) + sqrt(h) * Z[:,i];
    accept_prob = min(1, exp(U_fun(Xi_current)-U_fun(Xi_proposed)) * q(Xi_proposed, Xi_current) /q(Xi_current,Xi_proposed))
    if U[i] <= accept_prob
      Xi[:,i+1] = Xi_proposed
      acceptance_count += 1
    else
      Xi[:,i+1] = Xi_current
    end
  end
  println("MALA: Fraction of accepted proposals: ", acceptance_count / n_iter)
  return Xi
end

function consensusMALA(U_fun::Vector{Function}, grad_U::Vector{Function}, h::Real, dim::Integer, n_samples::Integer, Xi_init::Vector{Float64} = zeros(dim))
  n_batches = length(U_fun)
  samples = zeros(dim, n_samples, n_batches)
  W = Vector{Array{Float64,2}}(n_batches)  # weights
  cumulative_samples = zeros(dim, n_samples)
  cumulative_weights = zeros(dim,dim)
  for i=1:n_batches
    samples[:,:,i] = MALA(U_fun[i], grad_U[i], h, dim, n_samples, Xi_init)[:,2:end]
    # println(size(samples[:,:,i]))
    # print(cov(samples[:,:,i],2)); println()
    W[i] = inv(cov(samples[:,:,i],2))
    cumulative_samples += W[i] * samples[:,:,i]
    cumulative_weights += W[i]
  end
  println(cumulative_weights)
  consensus_samples = cumulative_weights \ cumulative_samples
end

max_stars_repo_name: jbierkens/zigzag-experiments

max_issues_repo_name: jbierkens/zigzag-experiments

max_forks_repo_name: jbierkens/zigzag-experiments


text:
function game(decisionMatrix::DataFrame; verbose=false)::GameResult

    newDecisionMatrix = similar(decisionMatrix)

    nrow, ncol = size(decisionMatrix)

    minmat = minimum(decisionMatrix)
    modified = false

    if minmat < 0
        modified = true
        newDecisionMatrix = decisionMatrix .- minmat
    end

    dm = convert(Matrix, newDecisionMatrix)

    model = Model(Cbc.Optimizer);
    MOI.set(model, MOI.Silent(), !verbose)

    @variable(model, x[1:nrow])
    @objective(model, Min, sum(x))

    for i in 1:nrow
        @constraint(model, sum(x[1:nrow] .* dm[:, i]) >= 1)
    end

    for i in 1:nrow
        @constraint(model, x[i] >= 0)
    end

    optimize!(model);

    values = JuMP.value.(x)

    gamevalue = 1 / objective_value(model)
    p = gamevalue * values

    if modified
        gamevalue -= abs(minmat)
    end

    result = GameResult(
        p,
        gamevalue
    )

    return result
end

max_stars_repo_name: ersagunkuruca/JMcDM

max_issues_repo_name: ersagunkuruca/JMcDM

max_forks_repo_name: ersagunkuruca/JMcDM


text:
using Pkg
Pkg.activate(@__DIR__)
using RCall
using SigmaRidgeRegression
using StatsBase
using Random
using MLJ
using Tables
using LaTeXTabulars
using LaTeXStrings
using BenchmarkTools

# Workaround until the following is resolved
# https://github.com/alan-turing-institute/MLJBase.jl/issues/428
# Else the session will hang upon printing the X matrix
# and wrapping it as Tables.MatrixTable
function Base.show(io::IO, tbl::Tables.MatrixTable)
    return Base.show(io, typeof(tbl))
end

# Load CLL dataset
cll_data = CLLData.load();
foldid = cll_data.foldid;
X = cll_data.X;
y = cll_data.y;
ngr = cll_data.ngr

resample_ids = [(findall(foldid .!= k), findall(foldid .== k)) for k in 1:10]

groups = GroupedFeatures(ngr)

Random.seed!(1)
rand_idx = sample(1:121, 121; replace=false)
Xdrug_resample = X[rand_idx, 1:305]
Xgaussian = randn(121, 100)
Xnoise = [X Xdrug_resample Xgaussian]

groups_noise = GroupedFeatures([ngr; 305; 100])

X_table = MLJ.table(X)
X_plus_noise = MLJ.table(Xnoise)

function single_table_line(
    X, y, resampling, _model, model_name; tuning_name=nothing, sigdigits=3
)
    ridge_machine = machine(_model, X, y)
    fit!(ridge_machine)
    best_param = ridge_machine.report.best_param
    if isnothing(tuning_name)
        tuning_string = ""
    else
        best_param = round(best_param; sigdigits=sigdigits)
        tuning_string = L"%$(tuning_name) = %$(best_param)"
    end

    λs = round.(deepcopy(ridge_machine.report.best_λs), sigdigits=sigdigits)

    ridge_benchmark = @benchmark fit!(machine($(_model), $(X), $(y)))

    time_loo = round(mean(ridge_benchmark.times) / (1e9); sigdigits=sigdigits)

    eval_ridge = evaluate!(ridge_machine; resampling=resampling, measure=rms)
    @show eval_ridge.measurement
    _rmse = round(eval_ridge.measurement[1]; sigdigits=sigdigits)

    return [model_name, tuning_string, λs..., time_loo, _rmse]
end

loo_single_ridge = LooRidgeRegressor(;
    ridge=SingleGroupRidgeRegressor(; groups=groups, center=true, scale=true)
)
loo_single_ridge_noise = LooRidgeRegressor(;
    ridge=SingleGroupRidgeRegressor(; groups=groups_noise, center=true, scale=true)
)

loo_sigmaridge = LooSigmaRidgeRegressor(; groups=groups, center=true, scale=true)
loo_sigmaridge_noise = LooSigmaRidgeRegressor(;
    groups=groups_noise, center=true, scale=true
)

loo_multi_ridge = LooRidgeRegressor(;
    ridge=MultiGroupRidgeRegressor(; groups=groups, center=true, scale=true),
    rng=MersenneTwister(1),
)

loo_multi_ridge_noise = LooRidgeRegressor(;
    ridge=MultiGroupRidgeRegressor(; groups=groups_noise, center=true, scale=true),
    rng=MersenneTwister(1),
)

cv_glasso = TunedSeagull(;groups=groups, center=true, scale=true)

cv_glasso_noise = TunedSeagull(;groups=groups_noise, center=true, scale=true)

line_single_ridge = single_table_line(
    X_table, y, resample_ids, loo_single_ridge, "\textbf{Single Ridge}"
)
line_single_ridge_noise = single_table_line(
    X_plus_noise, y, resample_ids, loo_single_ridge_noise, "\textbf{Single Ridge}"
)

line_sigma_ridge = single_table_line(
    X_table,
    y,
    resample_ids,
    loo_sigmaridge,
    L"$\sigmacv$\textbf{-Ridge}";
    tuning_name=L"\sigmacv",
)
line_sigma_ridge_noise = single_table_line(
    X_plus_noise,
    y,
    resample_ids,
    loo_sigmaridge_noise,
    L"$\sigmacv$\textbf{-Ridge}";
    tuning_name=L"\sigmacv",
)

line_multi_ridge = single_table_line(
    X_table, y, resample_ids, loo_multi_ridge, "\textbf{Multi Ridge}"
)

line_multi_ridge_noise = single_table_line(
    X_plus_noise, y, resample_ids, loo_multi_ridge_noise, "\textbf{Multi Ridge}"
)

line_glasso = single_table_line(
    X_table,
    y,
    resample_ids,
    cv_glasso,
    L"\textbf{Group Lasso}";
    tuning_name=L"$\widehat{\lambda}^{gl}",
)

line_glasso_noise = single_table_line(
    X_plus_noise,
    y,
    resample_ids,
    cv_glasso_noise,
    L"\textbf{Group Lasso}";
    tuning_name=L"$\widehat{\lambda}^{gl}",
)

tbl_spec = Tabular("lllllll")
line1 = [
    "",
    "Tuning",
    L"$\widehat{\lambda}_{\text{Drugs}}$",
    L"\widehat{\lambda}_{\text{Methyl}}",
    L"\widehat{\lambda}_{\text{RNA}}",
    "Time (s)",
    "10-fold RMSE",
]
lines = [line1, Rule(), line_sigma_ridge, line_single_ridge, line_multi_ridge, line_glasso]
latex_tabular("cll_analysis.tex", tbl_spec, lines)

tbl_spec2 = Tabular("lllllllll")
line1 = [
    "",
    "Tuning",
    L"$\widehat{\lambda}_{\text{Drugs}}$",
    L"\widehat{\lambda}_{\text{Methyl}}",
    L"\widehat{\lambda}_{\text{RNA}}",
    L"\widehat{\lambda}_{\text{Noise1}}",
    L"\widehat{\lambda}_{\text{Noise2}}",
    "Time (s)",
    "10-fold RMSE",
]
lines = [
    line1,
    Rule(),
    line_sigma_ridge_noise,
    line_single_ridge_noise,
    line_multi_ridge_noise,
    line_glasso_noise,
]
latex_tabular("cll_analysis_noise.tex", tbl_spec2, lines)


text:
module Choose

using ..Mattermost

function choose(req)
    options = split(req.text)
    # Weg mit dem !choose
    popfirst!(options)

    if length(options) == 0
        reply = string("@", req.user_name, ", du musst mir schon sagen was du mit !choose aussuchen möchtest. Das Kristallkugel-Modul bekomme ich erst in Version 2 :(")
    elseif length(options) == 1
        # Ja/Nein
        choice = rand(0:1)
        reply = string("@", req.user_name, ", ich sage: ", choice == 0 ? "Nein" : "Ja")
    else
        choice = 0
        legacy_joke = rand()
        if legacy_joke > 0.98
            push!(options, "Ja")
            choice = length(options)
        else
            choice = rand(1:length(options))
        end
        reply = string("@", req.user_name, ", ich sage: ", options[choice])
    end

    return OutgoingWebhookResponse(reply)
end

end


text:
using Test
using AirfoilMeshers
using LinearAlgebra
using Statistics
using ParametricAirfoils
using StaticArrays

mkpath(joinpath(@__DIR__, "output"))

include("test_tridiag.jl")
include("test_surface.jl")
include("test_hyperbolic.jl")


text:
@testset "adjust parameters to support sampler" begin

	cc = dummyChoiceContext()

	@testset "constructor, getparams, numparams, paramranges, setparams" begin

		subsA = DataGenerators.UniformSampler([7.1,9.2])
		s = DataGenerators.AdjustParametersToSupportSampler(subsA)

	    @testset "numparams and paramranges" begin
			@test DataGenerators.numparams(s) == 0
			prs = DataGenerators.paramranges(s)
	        @test typeof(prs) <: Vector{Tuple{Float64,Float64}}
			@test prs == Tuple{Float64,Float64}[]
	    end

	    @testset "default params" begin
			@test DataGenerators.getparams(s) == Float64[]
	    end

		@testset "setparams"  begin
			DataGenerators.setparams!(s, Float64[])
			@test_throws ErrorException DataGenerators.setparams!(s, [0.0])
		end

		@testset "validates subsampler to supported types" begin
			@test_throws ErrorException DataGenerators.AdjustParametersToSupportSampler(DataGenerators.GeometricSampler())
		end

	end

	@testset "sampling" begin

		subsA = DataGenerators.DiscreteUniformSampler([-1000.0, 214455.0])
		s = DataGenerators.AdjustParametersToSupportSampler(subsA)

		@mtestset "support is applied" reps=Main.REPS alpha=Main.ALPHA for cpsupport in [(188, 199), (-4, -4),]
			# note: includes case where lower bound = upper bound
			x, trace = DataGenerators.sample(s, cpsupport, cc)
			@test cpsupport[1] <= x <= cpsupport[2]
			@mtest_distributed_as DiscreteUniform(cpsupport[1],cpsupport[2]) x
		end

	end

	@testset "estimateparams" begin

		subsparams = [9.3, 15.8]
		subs1A = DataGenerators.UniformSampler(subsparams)
		s1 =DataGenerators.AdjustParametersToSupportSampler(subs1A)

		othersubsparams = [-55.1, -4.2]
		subs2A = DataGenerators.UniformSampler(othersubsparams)
		s2 =DataGenerators.AdjustParametersToSupportSampler(subs2A)

		traces = map(1:100) do i
			x, trace = DataGenerators.sample(s1, (0,1), cc)
			trace
		end

		estimateparams!(s2, traces)

		@test DataGenerators.getparams(s2) == Float64[]
		@test DataGenerators.getparams(subs2A) == othersubsparams
		# because no parameters of its own to estimate

	end

end


text:
"""
Calculate the best peak and centroid and their uncertainties using the median of the
distributions
"""
function bestcent(inpcent)
    perclim = 84.1344746
    centau = percentile(inpcent, 50)
    centau_uperr = (percentile(inpcent, perclim))-centau
    centau_loerr = centau-(percentile(inpcent, (100.0-perclim)))
    #print('Centroid, error: %10.3f  (+%10.3f -%10.3f)'%(centau, centau_loerr, centau_uperr))
    return centau,centau_loerr,centau_uperr
end

function centerror(centdistr)
    centau,centau_loerr,centau_uperr = bestcent(centdistr)
    return centau,centau_loerr,centau_uperr
end


text:
struct MapLevelFactor{S<:Unsigned, T<:Unsigned} <: AbstractFactor{T}
    name::String
    levels::AbstractVector{<:AbstractString}
    basefactor::AbstractFactor{S}
    transform::Function
    newindex::Vector{T}
end

Base.length(var::MapLevelFactor) = length(var.basefactor)

function MapLevelFactor(name::String, basefactor::AbstractFactor{S}, transform::Function) where {S<:Unsigned}
    baselevels = getlevels(basefactor)
    levelmap = Dict{String, S}()
    newlevels = Vector{String}()
    newindex = Vector{S}(length(baselevels) + 1)
    for (index, level) in enumerate(baselevels)
        newlevel = transform(level)
        if newlevel == MISSINGLEVEL
            newindex[index + 1] = zero(S)
        else
            levelcount = length(levelmap)
            levelindex = get(levelmap, newlevel, levelcount + 1)
            if levelindex > levelcount
                levelmap[newlevel] = levelindex
                newindex[index + 1] = levelindex
                push!(newlevels, newlevel)
            else
                newindex[index + 1] = levelindex
            end
        end
    end
    newmissing = transform(MISSINGLEVEL)
    if newmissing == MISSINGLEVEL
        newindex[1] = 0
    else
        levelcount = length(levelmap)
        levelindex = get(levelmap, newmissing, levelcount + 1)
        if levelindex > levelcount
            newindex[1] = levelindex
            push!(newlevels, newmissing)
        end
    end
    MapLevelFactor{S, S}(name, newlevels, basefactor, transform, newindex)
end

function slice(factor::MapLevelFactor{S, T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned} where {S<:Unsigned}
    newindex = factor.newindex
    f = (i::S -> newindex[i + 1])
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slices = slice(factor.basefactor, fromobs, toobs, slicelength)
    mapslice(f, slices, slicelength, T)
end

function Base.map(factor::MapLevelFactor{S, T}, dataframe::AbstractDataFrame; permute::Bool = false) where {S<:Unsigned} where {T<:Unsigned}
    MapLevelFactor(factor.name, map(factor.basefactor, dataframe; permute = permute), factor.transform)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text: #specific stuff for sp4 theory