text:
"""
afmize

simulate afm image from a given structure

example::
pdb = readpbd("1ake.pdb");
decenter!(pdb);
afm_image = afmize(pdb, (150.0, 150.0), (16, 16));
heatmap(afm_image) # using Plots
"""
function afmize_alpha(ta::TrjArray, box, npixel)::Matrix{Float64}
    rpixel = box ./ npixel
    nx = floor.(Int, (ta.x[1, :] .+ 0.5*box[1]) ./ rpixel[1]);
    ny = floor.(Int, (ta.y[1, :] .+ 0.5*box[2]) ./ rpixel[2]);
    #@show nx
    #@show ny
    afm_image = zeros(Float64, npixel)
    z_min = minimum(ta.z[1, :])
    z = ta.z[1, :] .- z_min
    for iatom in 1:ta.natom
        if (1 <= nx[iatom] <= npixel[1]) && (1 <= ny[iatom] <= npixel[2])
            if z[iatom] > afm_image[nx[iatom], ny[iatom]]
                afm_image[nx[iatom], ny[iatom]] = z[iatom]
            end
        end
    end
    return afm_image
end

function defaultParameters()
    return Dict(
     "H" => 1.20,
    "HE" => 1.40,
     "B" => 1.92,
     "C" => 1.70,
     "N" => 1.55,
     "O" => 1.52,
     "F" => 1.47,
    "NE" => 1.54,
    "NA" => 2.27,
    "MG" => 1.73,
    "AL" => 1.84,
    "SI" => 2.10,
     "P" => 1.80,
     "S" => 1.80,
    "CL" => 1.75,
    "AR" => 1.88,
     "K" => 2.75,
    "CA" => 2.31,
    "CYS" => 2.75,
    "PHE" => 3.2,
    "LEU" => 3.1,
    "TRP" => 3.4,
    "VAL" => 2.95,
    "ILE" => 3.1,
    "MET" => 3.1,
    "HIS" => 3.05,
    "HSD" => 3.05,
    "TYR" => 3.25,
    "ALA" => 2.5,
    "GLY" => 2.25,
    "PRO" => 2.8,
    "ASN" => 2.85,
    "THR" => 2.8,
    "SER" => 2.6,
    "ARG" => 3.3,
    "GLN" => 3.0,
    "ASP" => 2.8,
    "LYS" => 3.2,
    "GLU" => 2.95,
    )
end

mutable struct Point2D{T<:Real}
    x::T
    y::T
end

mutable struct Point3D{T<:Real}
    x::T
    y::T
    z::T
end

struct Probe
    x::Float64
    y::Float64
    r::Float64
    angle::Float64
end

mutable struct Sphere
    x::Float64
    y::Float64
    z::Float64
    r::Float64
end

mutable struct AfmizeConfig
    probeAngle::Float64        # Radian
    probeRadius::Float64
    range_min::Point2D
    range_max::Point2D
    resolution::Point2D
    atomRadiusDict::Dict{String, Float64}
end

function AfmizeConfig(r::Float64)
    return AfmizeConfig(10.0 * (pi / 180),
                        r,
                        Point2D(-90, -90),
                        Point2D(90, 90),
                        Point2D(15, 15),
                        defaultParameters())
end

function defaultConfig()
    return AfmizeConfig(10.0 * (pi / 180),
                        20.0,
                        Point2D(-90, -90),
                        Point2D(90, 90),
                        Point2D(15, 15),
                        defaultParameters())
end

function checkConfig(tra::TrjArray ,config::AfmizeConfig)::Union{String, Nothing}
    if (config.range_max.x - config.range_min.x) % config.resolution.x != 0
        return "resolution.x must divide range"
    end
    if (config.range_max.y - config.range_min.y) % config.resolution.y != 0
        return "resolution.y must divide range"
    end
    if config.probeAngle == 0
        return "probeAngle must be positive"
    end

    for name in tra.atomname
        if !haskey(config.atomRadiusDict, name)
            return "config dosen't know atom name $(name)"
        end
    end

    return nothing
end

# 底を0にする
function moveBottom(atoms::Array{Sphere})
    bottom = atoms[1].z - atoms[1].r
    for atom in atoms
        bottom = min(bottom, atom.z - atom.r)
    end
    for atom in atoms
        atom.z -= bottom
    end
end

# 針を球として見たとき、どこで衝突するかの計算
function calcCollisionAsSphere(probe::Probe, atom::Sphere)
    distXY = sqrt((probe.x - atom.x)^2 + (probe.y - atom.y)^2)
    dr = probe.r + atom.r
    if 0 < dr^2 - distXY^2
        return atom.z + sqrt(dr^2 - distXY^2) - probe.r
    else
        return 0.0
    end
end

# 針を円形推台として見たとき、どこで衝突するかの計算
function calcCollisionAsCircularThrusters(probe::Probe, atom::Sphere)
    distXY = sqrt((probe.x - atom.x)^2 + (probe.y - atom.y)^2)
    collisionDist = probe.r + atom.r * cos(probe.angle)
    if collisionDist < distXY
        return atom.z + atom.r * sin(probe.angle) - (distXY - collisionDist) / tan(probe.angle) - probe.r
    elseif probe.r < distXY
        return atom.z + sqrt(atom.r^2 - (distXY - probe.r)^2) - probe.r
    else
        return 0.0
    end
end

# 衝突計算をする範囲を返す
function calRectangle(min_point::Point2D, max_point::Point2D, config::AfmizeConfig)
    if max_point.x <= config.range_min.x || max_point.y <= config.range_min.y
        return nothing
    end
    if config.range_max.x <= min_point.x || config.range_max.y <= min_point.y
        return nothing
    end

    min_point = Point2D(max(min_point.x, config.range_min.x),
                        max(min_point.y, config.range_min.y))
    max_point = Point2D(min(max_point.x, config.range_max.x),
                        min(max_point.y, config.range_max.y))
    width  = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    return (max(1, floor(Int, (min_point.x - config.range_min.x + config.resolution.x - 1) / config.resolution.x)),
            max(1, floor(Int, (min_point.y - config.range_min.y + config.resolution.y - 1) / config.resolution.y)),
            min(width , floor(Int, (max_point.x - config.range_min.x + config.resolution.x - 1) / config.resolution.x)),
            min(height, floor(Int, (max_point.y - config.range_min.y + config.resolution.y - 1) / config.resolution.y)))
end

function calRectangle_sphere(atom::Sphere, config::AfmizeConfig)
    min_point = Point2D(atom.x - atom.r - config.probeRadius,
                        atom.y - atom.r - config.probeRadius)
    max_point = Point2D(atom.x + atom.r + config.probeRadius,
                        atom.y + atom.r + config.probeRadius)

    return calRectangle(min_point, max_point, config)
end

function calRectangle_circularThrusters(atom::Sphere, config::AfmizeConfig)
    dist_xy = (tan(config.probeAngle) * (atom.z + atom.r * sin(config.probeAngle) - config.probeRadius)
                + config.probeRadius + atom.r * cos(config.probeAngle))
    min_point = Point2D(atom.x - dist_xy,
                        atom.y - dist_xy)
    max_point = Point2D(atom.x + dist_xy,
                        atom.y + dist_xy)

    return calRectangle(min_point, max_point, config)
end

function afmize(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atoms = [Sphere(tra.x[i], tra.y[i], tra.z[i],
            config.atomRadiusDict[tra.atomname[i]]) for i in 1:tra.natom]
    moveBottom(atoms)

    stage = zeros(height, width)
    probes = [Probe(config.range_min.x + (w-0.5) * config.resolution.x,
                    config.range_min.y + (h-0.5) * config.resolution.y,
                    config.probeRadius, config.probeAngle)
             for h in 1:height, w in 1:width]

    # 各原子事に計算する価値のある矩形を求めて、その中で探索を行う
    for atom in atoms
        rectangle = calRectangle_sphere(atom, config)

        if isnothing(rectangle) continue end

        for h in rectangle[2]:rectangle[4], w in rectangle[1]:rectangle[3]
            probe = probes[h, w]

            stage[h, w] = max(stage[h, w], calcCollisionAsSphere(probe, atom))
        end
    end

    for atom in atoms
        rectangle = calRectangle_circularThrusters(atom, config)

        if isnothing(rectangle) continue end

        for h in rectangle[2]:rectangle[4], w in rectangle[1]:rectangle[3]
            probe = probes[h, w]

            stage[h, w] = max(stage[h, w], calcCollisionAsCircularThrusters(probe, atom))
        end
    end

    return stage
end

# 高速化前
function afmize_beta(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atoms = [Sphere(tra.x[i], tra.y[i], tra.z[i],
            config.atomRadiusDict[tra.atomname[i]]) for i in 1:tra.natom]
    moveBottom(atoms)

    stage = zeros(height, width)

    for h in 1:height, w in 1:width
        probe = Probe(config.range_min.x + (w-0.5) * config.resolution.x,
                      config.range_min.y + (h-0.5) * config.resolution.y,
                      config.probeRadius, config.probeAngle)
        for atom in atoms
            stage[h, w] = max(stage[h, w], calcCollisionAsSphere(probe, atom))
            stage[h, w] = max(stage[h, w], calcCollisionAsCircularThrusters(probe, atom))
        end
    end

    return stage
end

function afmize_gpu(tra::TrjArray, config::AfmizeConfig)
    message = checkConfig(tra, config)
    if !isnothing(message)
        println(message)
        return zeros(1, 1)
    end

    width = floor(Int, (config.range_max.x - config.range_min.x) / config.resolution.x)
    height = floor(Int, (config.range_max.y - config.range_min.y) / config.resolution.y)
    atom_r = tra.mass
    atom_z = tra.z[:]
    atom_z = atom_z .- (minimum(atom_z) - atom_r[argmin(atom_z)])
    #atom_z .= atom_z .- (minimum(atom_z))
    atom_z_max = maximum(atom_z)
    stage = similar(tra.x, 1, height*width)
    stage .= 0.0

    probe_x = similar(tra.x, 1, width)
    probe_y = similar(tra.x, 1, height)
    probe_x[:] .= range(config.range_min.x, step=config.resolution.x, stop=config.range_min.x + (width-1)*config.resolution.x) .+ (0.5*config.resolution.x)
    probe_y[:] .= range(config.range_min.y, step=config.resolution.y, stop=config.range_min.y + (height-1)*config.resolution.y) .+ (0.5*config.resolution.y)
    probe_r = config.probeRadius
    probe_angle = config.probeAngle

    # collision with sphere
    dist_x2 = (tra.x[:] .- probe_x).^2
    dist_y2 = (tra.y[:] .- probe_y).^2
    index_x_to_xy = [j for j=1:width for i=1:height]
    index_y_to_xy = [i for j=1:width for i=1:height]

    dist_xy2 = dist_x2[:, index_x_to_xy] .+ dist_y2[:, index_y_to_xy]
    #stage_xy = similar(dist_xy2)
    #stage_xy .= 0.0
    dr = probe_r .+ atom_r
    dr2 = dr.^2
    s = dr2 .- dist_xy2
    @show typeof(s)
    index_collide = s .> 0.0
    index_not_collide = .!index_collide
    @show typeof(index_collide)
    if any(index_collide)
        s[Array(index_collide)] .= sqrt.(s[index_collide])
    end
    s .= atom_z .+ s .- probe_r
    if any(index_not_collide)
        s[Array(index_not_collide)] .= 0.0
    end
    s_max = maximum(s, dims=1)
    stage .= max.(stage, s_max)

    # collision with circular thruster
    dist_xy2 .= sqrt.(dist_xy2)
    dist_collision = probe_r .+ atom_r .* cos(probe_angle)
    index_side   = dist_collision .< dist_xy2
    index_corner = (probe_r .< dist_xy2) .& .!index_side
    index_not_collide = .!(index_corner .| index_side)
    s[Array(index_side)] .= (atom_r .* sin.(probe_angle) .- (dist_xy2 .- dist_collision) ./ tan.(probe_angle))[index_side]
    s[Array(index_corner)] .= sqrt.( (atom_r.^2 .- (dist_xy2 .- probe_r).^2)[index_corner]  )
    s .= s .+ atom_z .- probe_r
    s[Array(index_not_collide)] .= 0.0
    s_max = maximum(s, dims=1)
    stage .= max.(stage, s_max)

    stage_reshape = reshape(stage, (height, width))
    return stage_reshape
end

function translateafm(afm, (dx, dy))
    afm_translated = zeros(eltype(afm), size(afm))
    (nx, ny) = size(afm)
    for i in maximum([1, 1-dx]):minimum([nx, nx-dx])
        for j in maximum([1, 1-dy]):minimum([ny, ny-dy])
            afm_translated[i+dx, j+dy] = afm[i, j]
        end
    end
    afm_translated
end

function translateafm_periodic(afm, (dx, dy))
    afm_translated = zeros(eltype(afm), size(afm))
    (nx, ny) = size(afm)
    for x in 1:nx
        for y in 1:ny
            xx = x + dx
            yy = y + dy
            if xx > nx
                xx -= nx
            end
            if yy > ny
                yy -= ny
            end
            afm_translated[xx, yy] = afm[x, y]
        end
    end
    afm_translated
end

function direct_convolution(observed, calculated)
    H, W = size(observed)
    ret = zeros(H, W)
    for y in 1:H, x in 1:W
        for dy in 1:H, dx in 1:W
            ny = y + dy - 1
            nx = x + dx - 1
            if ny > H ny -= H end
            if nx > W nx -= W end
            ret[y, x] += observed[dy, dx] * calculated[ny, nx]
        end
    end
    return ret
end

function fft_convolution(observed, calculated)
    return real.(ifft(fft(observed).*conj.(fft(calculated))))
end

function calcLogProb(observed::AbstractMatrix{T}, calculated::AbstractMatrix{T}, convolution_func) where {T}
    npix = eltype(observed)(size(observed, 1) * size(observed, 2))
    @show npix
    C_o  = sum(observed)
    C_c  = sum(calculated)
    C_cc = sum(calculated.^2)
    C_oo = sum(observed.^2)
    C_oc = convolution_func(observed, calculated)

    log01 = npix .* (C_cc .* C_oo .- C_oc.^2) .+ eltype(observed)(2) .* C_o .* C_oc .* C_c .- C_cc .* C_o.^2 .- C_oo .* C_c.^2
    #log01[Array(log01 .< log(eps(eltype(observed))))] .= eps(eltype(observed))
    log02 = (npix .- eltype(observed)(2)) .* (npix .* C_cc .- C_c.^2)
    #log02 = log02 <= log(eps(eltype(observed))) ? eps(eltype(observed)) : log02
    @show typeof(log01)
    @show typeof(log02)
    logprob = eltype(observed)(0.5) .* (eltype(observed)(3.0) .- npix) .* log.(log01) .+ (eltype(observed)(0.5) .* npix .- eltype(observed)(2.0)) .* log(log02)

    return logprob
end

function getafmposterior(afm::Matrix{Float64}, model_array::TrjArray, q_array::Matrix{Float64}, param_array)
    imax_model = 0
    imax_q = 0
    best_param = param_array[1]
    best_translate = [0, 0]
    best_afm = similar(afm)
    best_posterior = -Inf

    observed = afm
    npix = Float64(size(observed, 1) * size(observed, 2))

    decenter!(model_array)

    Threads.@threads for imodel in 1:size(model_array, 1)
        @show imodel
        model = model_array[imodel, :]
        for iq in 1:size(q_array, 1)
            q = q_array[iq, :]
            model_rotated = MDToolbox.rotate(model, q)
            for iparam in 1:length(param_array)
                param = param_array[iparam]
                calculated = afmize(model_rotated, param)
                C_o  = sum(observed)
                C_c  = sum(calculated)
                #@btime C_oc = sum(observed_translated .* calculated)
                C_oc = sum(observed .* calculated)
                C_cc = sum(calculated.^2)
                C_oo = sum(observed.^2)
                #@btime C_oc_dxdy = real.(ifftshift(ifft(fft(observed_translated).*conj.(fft(calculated)))))
                C_oc_dxdy = real.(ifftshift(ifft(fft(observed).*conj.(fft(calculated)))))
                log01 = npix .* (C_cc .* C_oo .- C_oc_dxdy.^2) .+ 2.0 .* C_o .* C_oc_dxdy .* C_c .- C_cc .* C_o.^2 .- C_oo .* C_c.^2
                log01[log01 .<= 0.0] .= eps(Float64)
                log02 = (npix .- 2.0) .* (npix .* C_cc .- C_c.^2)
                log02 = log02 <= 0 ? eps(Float64) : log02
                logprob = 0.5 .* (3.0 .- npix) .* log.(log01) .+ (0.5 .* npix .- 2.0) .* log.(log02)
                maximum_logprob = maximum(logprob)
                if best_posterior < maximum_logprob
                    best_posterior = maximum_logprob
                    imax_model = imodel
                    imax_q = iq
                    best_param = param
                    x_center = ceil(Int32, (size(observed,1)/2)+1.0)
                    y_center = ceil(Int32, (size(observed,2)/2)+1.0)
                    dx_estimated = argmax(logprob)[1] - x_center
                    dy_estimated = argmax(logprob)[2] - y_center
                    best_translate = (dx_estimated, dy_estimated)
                    best_afm = translateafm(calculated, best_translate)
                end
            end
        end
    end

    return imax_model, imax_q, best_param, best_translate, best_afm, best_posterior
end

mutable struct posteriorResult
    posteriors
    each_best
    afm_results
    posterior_results
    best_translate
    best_posterior
    best_model
    best_quate
    best_model_rotated
    best_param
    best_afm
end

function getafmposteriors_alpha(afm_frames, model_array, quate_array, param_array, opt = "")
    frame_num = size(afm_frames)[1]
    model_num = size(model_array)[1]
    quate_num = size(quate_array)[1]
    param_num = size(param_array)[1]

    convolution_func = fft_convolution
    if opt == "direct"
        convolution_func = direct_convolution
    end

    println("func is $convolution_func")

    results = []
    for i in 1:frame_num
        # 各モデルについて、一つの角度と半径ごとのlogprobの最高値を保持しておく
        posteriors = zeros(model_num, quate_num * param_num)
        # 各モデルについて、最も良い値をだした時のafm画像
        afm_results = [zeros(size(afm_frames[1])) for i in 1:model_num]
        each_best = ones(model_num) .* -Inf

        best_posterior = -Inf
        best_model = model_array[1]
        best_quate = quate_array[1, :]
        best_model_rotated = model_array[1]
        best_param = param_array[1]
        best_afm = zeros(size(afm_frames[1]))
        posterior_results = zeros(model_num)
        push!(results, posteriorResult(posteriors,
                                        each_best,
                                        afm_results,
                                        posterior_results,
                                        (0, 0),
                                        best_posterior,
                                        best_model,
                                        best_quate,
                                        best_model_rotated,
                                        best_param,
                                        best_afm))
    end

    for (model_id, model) in zip(1:model_num, model_array)
        @show model_id
        for quate_id in 1:quate_num
            rotated_model = MDToolbox.rotate(model, quate_array[quate_id, :])
            for param_id in 1:param_num
                cal_frame = afmize(rotated_model, param_array[param_id])
                for frame_id in 1:frame_num
                    prob_mat = calcLogProb(afm_frames[frame_id], cal_frame, convolution_func)
                    max_prob = maximum(prob_mat)
                    id = (quate_id - 1) * param_num + param_id
                    results[frame_id].posteriors[model_id, id] = max_prob

                    if results[frame_id].each_best[model_id] < max_prob
                        results[frame_id].each_best[model_id] = max_prob
                        results[frame_id].afm_results[model_id] = cal_frame
                    end

                    if results[frame_id].best_posterior < max_prob
                        results[frame_id].best_posterior = max_prob
                        results[frame_id].best_model = model
                        results[frame_id].best_quate = quate_array[quate_id, :]
                        results[frame_id].best_model_rotated = rotated_model
                        results[frame_id].best_param = param_array[param_id]
                        results[frame_id].best_translate = Tuple(argmax(prob_mat))
                        results[frame_id].best_afm = translateafm_periodic(cal_frame, results[frame_id].best_translate)
                    end
                end
            end
        end
    end

    for frame_id in 1:frame_num
        results[frame_id].posteriors .-= maximum(results[frame_id].posteriors)
        results[frame_id].posteriors = exp.(results[frame_id].posteriors)
        results[frame_id].posterior_results = sum(results[frame_id].posteriors, dims = 2)
        results[frame_id].posterior_results ./= sum(results[frame_id].posterior_results)
    end

    return results
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# # Clopper-Pearsonの信頼区間
#
# * 黒木玄
# * 2021-12-05
#
# $\newcommand\on{\operatorname}$
#
# 関連ノートブック
#
# * [止め方で結果が変わる？](https://nbviewer.org/github/genkuroki/public/blob/main/0025/%E6%AD%A2%E3%82%81%E6%96%B9%E3%81%A7%E7%B5%90%E6%9E%9C%E3%81%8C%E5%A4%89%E3%82%8F%E3%82%8B%EF%BC%9F.ipynb)
# * [ベイズハッキング](https://nbviewer.org/github/genkuroki/public/blob/main/0025/Bayes%20hacking.ipynb)

# %%
using Distributions
using StatsPlots
using Memoization
using Roots
using StatsFuns: logistic, logit

# %% [markdown]
# ## Clopper-Pearsonの信頼区間のP値函数の様々な表示
#
# ### 二項分布とベータ分布の関係
#
# 右辺の部分積分の繰り返しによって以下の公式を示すことができる(詳しい計算については「[止め方で結果が変わる？](https://nbviewer.org/github/genkuroki/public/blob/main/0025/%E6%AD%A2%E3%82%81%E6%96%B9%E3%81%A7%E7%B5%90%E6%9E%9C%E3%81%8C%E5%A4%89%E3%82%8F%E3%82%8B%EF%BC%9F.ipynb)」を参照):
#
# $$
# \begin{aligned}
# &
# \sum_{j=0}^k \binom{n}{j} \theta^j (1 - \theta)^{n-j} =
# \frac
# {\int_\theta^1 t^{(k+1)-1} (1 - t)^{(n-k)-1}\,dt}
# {B(k+1, n-k)} =
# \on{ccdf}(\on{Beta}(k+1, n-k), \theta),
# \\ &
# \sum_{j=k}^n \binom{n}{j} \theta^j (1 - \theta)^{n-j} =
# \frac
# {\int_0^\theta t^{k-1} (1 - t)^{(n-k+1)-1}\,dt}
# {B(k, n-k+1)} =
# \on{cdf}(\on{Beta}(k, n-k+1), \theta).
# \end{aligned}
# $$
#
# 前者の公式の左辺は二項分布 $\on{Binomial}(n, \theta)$ において成功回数が $k$ 以下になる確率であり, 「パラメータの値は $\theta$ 以上である」という帰無仮説の片側検定の $P$ 値とみなされる. 前者の公式の右辺はベイズ統計での対応する場合にimproper共役事前分布 $\on{Beta}(1, 0)$ から得られる事後分布で測ったパラメータの値が $\theta$ 以上になる確率になっている. 後者の公式は前者から導かれる.
#
# これらの公式は片側検定に関する通常のP値とimproper共役事前分布を用いたベイズ統計の事後分布で測った帰無仮説が成立する確率がぴったり等しいことを意味している.

# %% [markdown]
# ### ベータ分布を使ったClopper-Pearsonの信頼区間のP値函数の表示
#
# 上の結果を使うと二項分布モデルでの片側検定のP値の2倍で定義した両側検定のP値
#
# $$
# \on{pvalue\_dos}(n, k, \theta) =
# \min\left\{
# 1,\;
# 2\sum_{j=0}^k \binom{n}{j} \theta^j (1 - \theta)^{n-j},\;
# 2\sum_{j=k}^n \binom{n}{j} \theta^j (1 - \theta)^{n-j}
# \right\}
# $$
#
# をベータ分布の $\on{cdf}$, $\on{ccdf}$ で表すことができる:
#
# $$
# \begin{aligned}
# \on{pvalue\_dos}(n, k, \theta) &=
# \min\left\{
# 1,\;
# 2\frac
# {\int_\theta^1 t^{(k+1)-1} (1 - t)^{(n-k)-1}\,dt}
# {B(k+1, n-k)},\;
# 2\frac
# {\int_0^\theta t^{k-1} (1 - t)^{(n-k+1)-1}\,dt}
# {B(k, n-k+1)}
# \right\}
# \\ &=
# \min\{
# 1,\;
# 2\on{ccdf}(\on{Beta}(k+1, n-k), \theta),\;
# 2\on{cdf}(\on{Beta}(k, n-k+1), \theta)
# \}.
# \end{aligned}
# $$
#
# ここで $\on{dos}$ は doubuled one-side の略である. この $\on{pvalue\_dos}(n, k, \theta)$ は __Clopper-Pearsonの信頼区間__ に対応するP値函数になっている.

# %% [markdown]
# ### F分布とベータ分布の関係
#
# 自由度 $\nu_1, \nu_2$ のF分布の確率密度函数の $dx$ 倍は以下のように書ける:
#
# $$
# \on{pdf}(\on{FDist}(\nu_1, \nu_2), x)\,dx =
# \frac{1}{B(\nu_1/2, \nu_2/2)}
# \left(\frac{\nu_1 x}{\nu_1 x + \nu_2}\right)^{\nu_1/2}
# \left(1 - \frac{\nu_1 x}{\nu_1 x + \nu_2}\right)^{\nu_2/2}
# x^{-1}
# \,dx \quad (x > 0).
# $$
#
# これを
#
# $$
# t = \frac{\nu_1 x}{\nu_1 x + \nu_2}, \quad
# x = \frac{\nu_2}{\nu_1}\frac{t}{1 - t}, \quad 0 < t < 1
# $$
#
# で変数変換すると,
#
# $$
# \begin{aligned}
# \on{pdf}(\on{FDist}(\nu_1, \nu_2), x)\,dx &=
# \frac{1}{B(\nu_1/2, \nu_2/2)}
# t^{\nu_1/2}
# (1 - t)^{\nu_2/2}
# \,\frac{\nu_1}{\nu_2}\frac{1-t}{t}
# \,\frac{\nu_2}{\nu_1}\frac{dt}{(1-t)^2}
# \\ &=
# \frac{1}{B(\nu_1/2, \nu_2/2)} t^{\nu_1/2 - 1} (1 - t)^{\nu_2/2 - 1} \,dt
# \\ &=
# \on{pdf}(\on{Beta}(\nu_1/2, \nu_2/2), t)\,dt.
# \end{aligned}
# $$
#
# これは上の変数変換によって, 自由度 $\nu_1, \nu_2$ のF分布がパラメータ $\nu_1/2, \nu_2/2$ のベータ分布に変換されることを意味している. ゆえに
#
# $$
# \theta = \frac{\alpha x}{\alpha x + \beta}, \quad
# x = \frac{\beta}{\alpha}\frac{\theta}{1 - \theta}
# $$
#
# とおくと,
#
# $$
# \on{cdf}(\on{Beta}(\alpha, \beta), \theta) = \on{cdf}(\on{FDist}(2\alpha, 2\beta), x), \quad
# \on{ccdf}(\on{Beta}(\alpha, \beta), \theta) = \on{ccdf}(\on{FDist}(2\alpha, 2\beta), x).
# $$

# %% [markdown]
# ### F分布を使ったClopper-Pearsonの信頼区間のP値函数の表示
#
# 以上の結果を使えば __Clopper-Pearsonの信頼区間__ に対応するP値函数をF分布を使って計算することもできる:
#
# $$
# \on{pvalue\_dos}(n, k, \theta) =
# \min\begin{Bmatrix}
# 1 \\
# 2\on{ccdf}\left(\on{FDist}(2(k+1), 2(n-k)), \frac{n-k}{k+1}\frac{\theta}{1-\theta}\right) \\
# 2\on{cdf}\left(\on{FDist}(2k, 2(n-k+1)), \frac{n-k+1}{k}\frac{\theta}{1-\theta}\right)
# \end{Bmatrix}.
# $$

# %%
@memoize function pvalue_dos_naive(n, k, θ)
    bin = Binomial(n, θ)
    min(1, 2sum(pdf(bin, j) for j in 0:k), 2sum(pdf(bin, j) for j in k:n))
end

@memoize function pvalue_dos(n, k, θ)
    k ≤ 0 && return min(1, 2ccdf(Beta(k+1, n-k), θ))
    k ≥ n && return min(1, 2cdf(Beta(k, n-k+1), θ))
    min(1, 2ccdf(Beta(k+1, n-k), θ), 2cdf(Beta(k, n-k+1), θ))
end

@memoize function pvalue_dos_fdist(n, k, θ)
    k ≤ 0 && return min(1, 2ccdf(FDist(2(k+1), 2(n-k)), (n-k)/(k+1)*θ/(1-θ)))
    k ≥ n && return min(1, 2cdf(FDist(2k, 2(n-k+1)), (n-k+1)/k*θ/(1-θ)))
    min(1,
        2ccdf(FDist(2(k+1), 2(n-k)), (n-k)/(k+1)*θ/(1-θ)),
        2cdf(FDist(2k, 2(n-k+1)), (n-k+1)/k*θ/(1-θ)))
end

# %%
# Clopper-Pearsonの信頼区間に対応する3つのP値函数が一致することの確認

n = 10
k = 0:n
θ = 0:0.1:1

p1 = pvalue_dos_naive.(n, k, θ')

# %%
p2 = pvalue_dos.(n, k, θ')

# %%
p1 .≈ p2

# %%
p3 = pvalue_dos_fdist.(n, k, θ')

# %%
p3 .≈ p2

# %% [markdown]
# ### 負の二項分布とベータ分布の関係
#
# 上と同様に右辺の部分積分の繰り返しによって次の公式も示すことができる:
#
# $$
# \sum_{m=n}^\infty \binom{m-1}{k-1} \theta^k (1 - \theta)^{m-k} =
# \frac
# {\int_\theta^1 t^{k-1} (1 - t)^{(n-k)-1}\, dt}
# {B(k, n-k)}
# \quad (n \ge k \ge 1,\ 0 < \theta \le 1).
# $$
#
# 左辺の $\binom{m-1}{k-1} \theta^k (1 - \theta)^{m-k}$ は成功確率を意味するパラメータの値が $\theta$ のBernoulli試行をちょうど $k$ 回成功するまで繰り返したときの試行回数がちょうど $m$ 回になる確率なので, 左辺の和はちょうど $k$ 回の成功するまでの試行回数が $n$ 回以上になる確率になっており, 「パラメータの値は $\theta$ 以上である」という帰無仮説の片側検定の $P$ 値とみなされる.  右辺はベイズ統計での対応する場合にimproper事前分布 $\on{Beta}(0, 0)$ から得られる事後分布で測ったパラメータの値が $\theta$ 以上になる確率になっている.
#
# この公式も片側検定に関する通常のP値とimproper共役事前分布を用いたベイズ統計の事後分布で測った帰無仮説が成立する確率がぴったり等しくなっていることを意味している.

# %% [markdown]
# `Distributions.jl` における負の二項分布の定義は以下の通り:
#
# $$
# \on{pdf}(\on{NegativeBinomial}(r, \theta), k) = \binom{k+r-1}{k} p^r (1 - p)^k \quad (k=0,1,2,\ldots).
# $$
#
# ゆえに
#
# $$
# \on{pdf}(\on{NegativeBinomial}(k, \theta), m-k) =
# \binom{m-1}{m-k} \theta^k (1 - \theta)^{m-k} =
# \binom{m-1}{k-1} \theta^k (1 - \theta)^{m-k} \quad (m\ge k).
# $$

# %%
@memoize pvalue_negbin_naive(n, k, θ) =
    1 - sum(pdf(NegativeBinomial(k, θ), m - k) for m in k:n-1)

@memoize pvalue_negbin(n, k, θ) = ccdf(Beta(k, n-k), θ)

# %%
# 負の二項分布モデルでの片側検定に関する2つのP値函数が一致することの確認

k = 5
n = k+1:15
θ = 0.1:0.1:1

p1 = pvalue_negbin_naive.(n, k, θ')

# %%
p2 = pvalue_negbin.(n, k, θ')

# %%
p1 .≈ p2

# %% [markdown]
# ### もう１つのP値函数
#
# 上の片側確率の2倍版と異なる二項分布の両側検定のP値を次のように定めることもできる:
#
# $$
# \on{pvalue\_exact}(n, k, \theta) =
# \sum_{0\le j\le n,\; P(n,j,\theta)\le P(n,k,\theta)} P(n, k, \theta).
# $$
#
# ここで $
# P(n, k, \theta) = \binom{n}{j} \theta^j (1 - \theta)^{n-j}
# $.

# %%
x ⪅ y = x < y || x ≈ y # 以下の計算で ≤使うと失敗する

@memoize function pvalue_exact(n, k, θ)
    bin = Binomial(n, θ)
    p0 = pdf(bin, k)
    sum(pdf(bin, j) for j in support(bin) if pdf(bin, j) ⪅ p0)
end

# %%
# 別のP値函数が定義されたことの確認

n = 4
k = 0:n
θ = 0:0.1:1

p1 = pvalue_exact.(n, k', θ)
p2 = pvalue_dos.(n, k', θ)
p1 - p2

# %% [markdown]
# ## 信頼区間
#
# ### P値函数を用いた一般的な信頼区間の定義の仕方
#
# これらのＰ値函数 $\on{pvalue\_func} = \on{pvalue\_dos}, \on{pvalue\_exact}$ から信頼度 $1-\alpha$ の信用区間を次のように定義できる:
#
# $$
# \on{confidence\_interval}(\on{pvalue\_func}, n, k, \alpha) =
# \{\, \theta \in [0, 1] \mid \on{pvalue\_func}(n, k, \theta) \ge \alpha \,\}.
# $$
#
# これは $\on{pvalue\_func}(n, k, \theta) = \alpha$ となる $\theta$ を得ることによって数値計算できる.

# %%
@memoize function confidence_interval(pvalue_func, n, k, α = 0.05)
    f(t) = pvalue_func(n, k, logistic(t)) - α
    ci = logistic.(find_zeros(f, -10, 10))
    length(ci) < 2 && return 2k ≤ n ? (0.0, first(ci)) : (first(ci), 1.0)
    (first(ci), last(ci))
end

# %%
[(k, confidence_interval(pvalue_dos, 10, k)) for k in 0:10]

# %%
[(k, confidence_interval(pvalue_exact, 10, k)) for k in 0:10]

# %% [markdown]
# ### Clopper-Pearsonの信頼区間のベータ分布やF分布を用いた表示
#
# 特に $\on{pvalue\_func} = \on{pvalue\_dos}$ の場合には上に述べたことから, 以下が成立することがわかる:
#
# $$
# \on{confidence\_interval}(\on{pvalue\_dos}, n, k, \alpha) = [\theta_L, \theta_U]
# $$
#
# ここで
#
# $$
# \begin{aligned}
# &
# \on{ccdf}(\on{Beta}(k+1, n-k), \theta_U) =
# \frac
# {\int_{\theta_U}^1 t^{(k+1)-1} (1 - t)^{(n-k)-1}\,dt}
# {B(k+1, n-k)} = \frac{\alpha}{2},
# \\ &
# \on{cdf}(\on{Beta}(k, n-k+1), \theta_L) =
# \frac
# {\int_0^{\theta_L} t^{k-1} (1 - t)^{(n-k+1)-1}\,dt}
# {B(k, n-k+1)} = \frac{\alpha}{2}.
# \end{aligned}
# $$
#
# これらを満たす $\theta_L, \theta_U$ はBeta分布に関する $\on{cquantile}$, $\on{quantile}$ 函数で計算できる:
#
# $$
# \theta_U = \on{cquantile}(\on{Beta}(k+1, n-k), \alpha/2), \quad
# \theta_L = \on{quantile}(\on{Beta}(k, n-k+1), \alpha/2).
# $$
#
# $\on{confidence\_interval}(\on{pvalue\_dos}, n, k, \alpha)$ を __Clopper-Pearsonの信頼区間__ と呼ぶ.

# %% [markdown]
# Clopper-Pearsonの信頼区間をF分布を用いて表示することもできる:
#
# $$
# \begin{aligned}
# &
# \theta_U = \frac{(k+1)x_U}{(k+1)x_U + (n-k)}, \quad
# x_U = \on{cquantile}(\on{FDist}(2(k+1), 2(n-k)), \alpha/2),
# \\ &
# \theta_L = \frac{k x_L}{k x_L + (n-k+1)}, \quad
# x_L = \on{quantile}(\on{FDist}(2k, 2(n-k+1)), \alpha/2).
# \end{aligned}
# $$

# %%
@memoize function confidence_interval_dos(n, k, α = 0.05)
    θ_U = k ≥ n ? 1.0 : cquantile(Beta(k+1, n-k), α/2)
    θ_L = k ≤ 0 ? 0.0 :  quantile(Beta(k, n-k+1), α/2)
    (θ_L, θ_U)
end

@memoize function confidence_interval_dos_fdist(n, k, α = 0.05)
    x_U = k ≥ n ? 1.0 : cquantile(FDist(2(k+1), 2(n-k)), α/2)
    x_L = k ≤ 0 ? 0.0 :  quantile(FDist(2k, 2(n-k+1)), α/2)
    θ_U = (k+1)*x_U/((k+1)*x_U + (n-k))
    θ_L = k*x_L/(k*x_L + (n-k+1))
    (θ_L, θ_U)
end

# %%
# confidence_interval(pvalue_dos, n, k) と confidence_interval_dos(n, k) が等しいことの確認

[(k, confidence_interval(pvalue_dos, 10, k) .≈ confidence_interval_dos(10, k)) for k in 0:10]

# %%
# confidence_interval(pvalue_dos, n, k) と confidence_interval_dos_fdist(n, k) が等しいことの確認

[(k, confidence_interval(pvalue_dos, 10, k) .≈ confidence_interval_dos_fdist(10, k)) for k in 0:10]

# %%
# 多重ディスパッチで効率化 (`find_zeros` 函数による計算は効率が悪い)

confidence_interval(::typeof(pvalue_dos), n, k, α = 0.05) = confidence_interval_dos(n, k, α)

# %%
methods(confidence_interval)

# %%
# 多重ディスパッチによる効率化の確認

@which confidence_interval(pvalue_dos, 10, 1)

# %% [markdown]
# ## P値函数と信頼区間のプロット
#
# ### 与えられたデータに対するP値函数と信頼区間のプロット

# %%
"""与えられたデータに対するP値函数と信頼区間のプロット"""
function plot_pvalue_funcs(n, k, α = 0.05)
    m = k/n
    s = √(m*(1 - m)/n)
    a, b = max(0, m - 4s), min(1, m + 4s)
    θ = range(a, b; length=1000)
    p_exact = pvalue_exact.(n, k, θ)
    p_dos   = pvalue_dos.(n, k, θ)

    ci_exact = confidence_interval(pvalue_exact, n, k, α) |> collect
    ci_dos   = confidence_interval(pvalue_dos,   n, k, α) |> collect

    plot(;)
    plot!(θ, p_exact; label="pvalue_exact")
    plot!(θ, p_dos; label="pvalue_dos", ls=:dash)
    plot!(ci_exact, fill(α + 0.005, 2); label="CI for pvalue_exact", c=:blue, lw=2.5)
    plot!(ci_dos,   fill(α - 0.005, 2); label="CI for pvalue_dos", c=:red, lw=2.5, ls=:dash)
    plot!(; ytick=0:0.05:1)
    plot!(; xlabel="θ", ylabel="P-value")
    title!("n = $n, k = $k, α = $α"; titlefontsize=12)
end

# %%
plot_pvalue_funcs(10, 4)

# %%
plot_pvalue_funcs(100, 40)

# %%
plot_pvalue_funcs(1000, 400)

# %% [markdown]
# ### 第一種の過誤が生じる確率

# %%
"""有効サポート"""
function effective_support(bin::Binomial; m = 4)
    μ, σ = mean(bin), std(bin)
    a = max(minimum(bin), round(Int, μ - m*σ))
    b = min(maximum(bin), round(Int, μ + m*σ))
    a:b
end

"""二項分布における期待値"""
@memoize function expectation_value(f, n, θ)
    bin = Binomial(n, θ)
    sum(k -> f(k)*pdf(bin, k), effective_support(bin))
end

"""二項分布における確率"""
probability(f, n, θ) = expectation_value(f, n, θ)

"""第一種の過誤が生じる確率をαを動かしてプロット"""
function plot_prob_typeIerror(n, θ)
    p_exact(α) = probability(k -> pvalue_exact(n, k, θ) ⪅ α, n, θ)
    p_dos(α)   = probability(k -> pvalue_dos(n, k, θ) ⪅ α,   n, θ)

    α = 0:0.001:1
    tick = 0:0.1:1
    P = plot(; legend=:topleft)
    plot!(α, p_exact; label="pvalue_exact")
    plot!(α, p_dos;   label="pvalue_dos", ls=:dash)
    plot!([0, maximum(α)], [0, maximum(α)]; label="", c=:black, ls=:dot)
    plot!(; xtick=tick, ytick=tick)
    plot!(; xlabel="α", ylabel="probabilty of type I error")

    α = 0:0.0001:0.1
    tick = 0:0.01:0.1
    Q = plot(; legend=:topleft)
    plot!(α, p_exact; label="pvalue_exact")
    plot!(α, p_dos;   label="pvalue_dos", ls=:dash)
    plot!([0, maximum(α)], [0, maximum(α)]; label="", c=:black, ls=:dot)
    plot!(; xtick=tick, ytick=tick)
    plot!(; xlabel="α", ylabel="probabilty of type I error")

    plot(P, Q; size=(800, 400), leftmargin=3Plots.mm)
end

# %% [markdown]
# αごとに第一種の過誤が実際に起こる確率のグラフが45度線に近いP値函数の方がより正確であると考えられる.
#
# 以下を見れば分かるように, `pvalue_exact` の方が `pvalue_dos` よりも正確である.

# %%
@time plot_prob_typeIerror(10, 0.4)

# %%
@time plot_prob_typeIerror(100, 0.4)

# %%
@time plot_prob_typeIerror(1000, 0.4)

# %% [markdown]
# ### 信頼区間の被覆確率

# %%
"""区間にθが含まれるかどうかを判定する函数"""
contains(interval, θ) = first(interval) ≤ θ ≤ last(interval)

"""信頼区間に真の値が含まれる確率をプロット"""
function plot_coverage_prob(n, α = 0.05, ytick=0:α/5:1)
    c_exact(θ) = probability(k -> contains(confidence_interval(pvalue_exact, n, k, α), θ), n, θ)
    c_dos(θ)   = probability(k -> contains(confidence_interval(pvalue_dos,   n, k, α), θ), n, θ)

    θ = 0:0.001:1
    y_exact = similar(θ)
    y_dos   = similar(θ)
    Threads.@threads for i in eachindex(θ)
        y_exact[i] = c_exact(θ[i])
        y_dos[i]   = c_dos(θ[i])
    end

    plot(; legend=:bottomright, ylim=(1 - 1.25α, 1 + 0.05α))
    plot!(θ, y_exact; label="pvalue_exact")
    plot!(θ, y_dos;   label="pvalue_dos", ls=:dash)
    hline!([1 - α]; label="", c=:black, ls=:dot)
    plot!(; xtick=0:0.1:1, ytick)
    plot!(; xlabel="θ", ylabel="coverage probability")
    title!("coverage probability for n = $n, α = $α"; titlefontsize=10)
end

# %% [markdown]
# 信頼区間に含まれる確率を被覆確率(coverage probability)と呼ぶ.
#
# 被覆確率が $1-\alpha$ に近い信頼区間の方がより正確であると考えられる.
#
# 以下を見れば分かるように, pvalue_exact の信頼区間の方が pvalue_dos の信頼区間(Clopper-Pearsonの信頼区間)よりも正確である.

# %%
@time plot_coverage_prob(10)

# %%
@time plot_coverage_prob(100)

# %%
@time plot_coverage_prob(1000)

# %%

ext: jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
let
    f(x) = x⋅x
    g(x) = 2x
    x = [1.0, -1.0]
    state = Optini.FirstOrderState(x, f(x), g(x))
    p = -state.∇f

    @testset "Line Search Initial Step Length" begin
        inits = [
            StaticInitial(1f0),
            PreviousDecreaseInitial(1f0),
            QuadraticInitial(1f0)
        ]

        for init in inits
            α = init(state, p)
            @test typeof(α) == Float32 # type stability
            @test α == 1f0 # first iteration
            Optini.update!(init, state, p, 0.1f0)
        end

        next_x = [0.8, -0.8]
        next_state = Optini.FirstOrderState(next_x, f(next_x), g(next_x))
        next_p = -next_state.∇f
        steps = [1f0, 0.15625f0, 0.2840625f0] # expected initial step lengths

        for (init, step) in zip(inits, steps)
            α = init(next_state, next_p)
            @test typeof(α) == Float32 # type stability
            @test α == step # check initial step
        end
    end

    @testset "Line Search Methods" begin
        linesearches = [
            StaticLineSearch(),
            ExactLineSearch(reltol=1e-7),
            BacktrackingLineSearch(),
            InterpolationLineSearch(),
            StrongWolfeLineSearch()
        ]

        for (iter, ls) in enumerate(linesearches)
            α = ls(f, state, p, 1f0)
            @test typeof(α) == Float32 # type stability
            if iter == 1
                @test α == 1f0 # check step length
            else
                @test α == 0.5f0 # check step length
            end
        end
    end
end


text:
using HTSLib
import Downloads

"""
    count_remote_file_records(url::AbstractString)

Iterate over the bam file. Use `curl` or `gsutil` to stream the data.
"""
function count_remote_file_records(url::AbstractString)
    if startswith(url, "http://") || startswith(url, "https://")
        stream = open(pipeline(`curl --silent $url`), "r")
    elseif startswith(url, "gs://")
        stream = open(pipeline(`gsutil cat $url`), "r")
    else
        # assume this is a local file
        stream = open(url, "r")
    end

    ## htsopen may use arbitrary IO object in Julia
    count_remote_file_records(stream)
end

function count_remote_file_records(io::IO)
    nreads::Int = 0
    for record in HTSReadWriter(io, "r")
        nreads += 1
    end
    println("number of reads: $nreads")
end

# sample test files:
#   - http://hgdownload.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeUwRepliSeq/wgEncodeUwRepliSeqBg02esG1bAlnRep1.bam
#   - gs://broad-public-datasets/NA12878/NA12878.cram

if PROGRAM_FILE == relpath(@__FILE__)
    @assert length(ARGS) == 1
    url = ARGS[1]
    count_remote_file_records(url)
end


text:
using Printf
using DelimitedFiles
using Plots

const year_seconds = 31556952

### open filenames and retrun list of files, and list of arrays with strings of data
function open_files(filenames...; begin_index=2, final_index=nothing)
    files = []
    plot_data = []
    for filename in filenames
        file = open(filename, "r")
        push!(files, file)
        temp_data = collect(eachline(file))
        temp_data = temp_data[2:end, :]
        if final_index == nothing
            final_index = size(temp_data)[1]
        end
        data = @view temp_data[begin_index:final_index, :]
        push!(plot_data, data)
    end
    return files, plot_data
end

### parse line of 2D array of strings data into 1D array of Float64s
get_line(data, i) = map(x -> parse(Float64, x), split.(data[i]))

### parse single string into 1D array of FLoat64s
parse_line(line) = map(x -> parse(Float64,x), split.(line))

### take off the first two lines of the slip_file and return them and the new array
function get_y_nn(slip_data)

    nn = parse(Int64, split.(Iterators.take(slip_data, 1))[1][1])
    y = map(x -> parse(Float64, x), split.(slip_data[2]))[3:end]
    slip_data = @view slip_data[3:end, :]

    return y, nn, slip_data

end

### get the number of breaks in a file, and the starting indices of each file
function get_break_indices(slip_data)

    cycle_index = [1]
    count = 0

    for (i, line) in enumerate(slip_data)
        if line == "BREAK"
            count += 1
            push!(cycle_index, i+1)
        end
    end

    return count, cycle_index

end

### make slip contour plot with given slip_data
function plot_slip(slip_data, y, nn, title)

    plt = plot(legend=false, yflip = true, ylabel="Depth(Km)", xlabel="Slip(m)", title=title)

    δ = zeros(nn)
    line = zeros(nn+2)
    break_count = 0
    off_line = get_line(slip_data, 1)
    t = off_line[1]
    δ_off = off_line[3:end]
    @printf "begining time: %f\n" t/year_seconds
    final_index = size(slip_data)[1]
    for index in 1:final_index

        @printf "\r%f" index/final_index
        if slip_data[index] == "BREAK"
            plt = plot!(δ-δ_off, y, linecolor=:black, linewidth = 1)
            break_count += 1
        else
             if break_count % 2 == 0
                 line .= get_line(slip_data, index)
                 δ .= line[3:end]
                 t = line[1]
                 plt = plot!(δ-δ_off, y, linecolor=:blue, linewidth = .2)
             elseif break_count % 2 == 1 && index % 3 == 1
                 δ .= get_line(slip_data, index)[3:end]
                 t = line[1]
                 plt = plot!(δ-δ_off, y, linecolor=:red, linewidth = .2)
            end
        end

    end

    @printf "\nfinal time: %f\n" t/year_seconds

    return plt
end

### Create fault animations, with data gotten from open_files. Assumes all data is of the same length.
function fault_animations(plot_data, y)

    an_length = size(plot_data[1])[1]

    for i in 1:an_length
        plots = []
        for data in plot_data
            line = get_line(data, i)[3:end]
            plt = plot(line, y, yflip=true, legend=false)
            push!(plots, plt)
        end
        plot(plots..., layout = (length(plots), 1))
        gui()
    end
end

### convert begin_cycle, and final_cycle to there corrisponding indices in slip_data
function get_plot_indices(begin_cycle, final_cycle, break_indices)

    index_offset = break_indices[2*begin_cycle - 1]
    if 2*final_cycle + 1 <= length(break_indices)
        final_index = break_indices[2*final_cycle + 1] - 2
    elseif 2*final_cycle + 1 > length(break_indices)
        final_index = break_indices[2*final_cycle] - 2
    end

    return final_index, index_offset

end

### write out new initial conditions
function write_ψδ(ψ, δ, t)

    ψδ = vcat(ψ, δ)

    @printf "Filename: "
    name = chomp(readline())
    filename = string("/home/tharvey2/Basin_Simulations/input_files/initial_cons/", name)
    filename_t = string(filename, "_t")

    fout = open(filename, "w")
    writedlm(fout , ψδ)
    close(fout)

    fout = open(filename_t, "w")
    writedlm(fout, [t])
    close(fout)

    @printf "inital conditions write to %s\n" filename

end

### get coordinates, and data from volume files
function parse_volume_file(file_name, cycle_number, break_indices)

    file = open(file_name, "r")
    iter = eachline(file)

    # get volume coordinates
    x, iter = firstrest(iter)
    x = parse_line(x)
    y, iter = firstrest(iter)
    y = parse_line(y)

    iter = Iterators.drop(iter, 1)

    return file, x, y, iter

end


text:
using StatsBase, JLD,  Hwloc,  DynNets, StatsFuns,Utilities

#initWorkers()
#@everywhere  using DynNets, StatsFuns, Utilities
N_est = 2
NGWpoints = 1
Tvals =  [  600 400 200 ]#[50 250 1250]#
N_Tind = length(Tvals)[1]
scoreRescType = ""#"FISHER-EWMA"#
useStartVal = false# true# false
N_Nmin = 500
Nspacing =  500
Nmax = 10000
Nvals =  Vector(Nmin:Nspacing:Nmax)
Nind = length(Nvals)

est_times = SharedArray{Float64}(zeros(N_Tind,N_Nind,NGWpoints,N_est)) #zeros(N_Tind,N_Nind,NGWpoints,N_est)#
est_conv_flag =  SharedArray{Bool,4}(N_Tind,N_Nind,NGWpoints,N_est) #falses(N_Tind,N_Nind,NGWpoints,N_est) #
est_flag = SharedArray{Bool,4}(N_Tind,N_Nind,NGWpoints,N_est) #falses(N_Tind,N_Nind,NGWpoints,N_est) #
GWvalStore  = Array{Array{Int,1},1}(N_Nind)

@everywhere function loopEval(T,N,simParVec, useStartVal,groupsindTsDgp,scoreRescType,GBAest,GWest,groupsindTsEst)
    Model2Sample = DynNets.SdErgm1(ones(T,N),simParVec,
                                                groupsindTsDgp,scoreRescType)
    Model2Sample,Y_T,Fitness_T  =  DynNets.sampl(Model2Sample,T)
    sampleDegsT = squeeze(sum(Y_T,2),2)

    Model2Est =  DynNets.SdErgm1(sampleDegsT,[zeros(GWest),zeros(GBAest),zeros(GBAest)],
                                                groupsindTsEst,scoreRescType)
    #estimate
    useStartVal ?   start_val = simParVec : start_val = zeros(3,3)
    estTime = @elapsed  hat_pars , conv_flag =
                        DynNets.estimate(Model2Est,start_values = start_val )
    return estTime,conv_flag
end

save_fold = "./data/estimatesTest/compComplexity/"
file_name = "EstTest_varN_Nest_$(N_est)_$(Nmin)_$(Nspacing)_$(Nmax)_varGW_$(NGWpoints)_T_$(Tvals)_corectSpec"* scoreRescType *".jld"
save_path = save_fold*file_name#

#

GW = 1
@progress "N" for indN = 1:N_Nind
#indN = 1
N = Nvals[indN]
#GWvals = unique(Vector{Int}(floor.(linspace(1,N,NGWpoints))))
#GWvalStore[indN] = GWvals

#@progress "GW"  for indGW = NGWpoints:-1:1
indGW = 1
#GW =  GWvals[indGW]
GBAdgp =1; GWdgp = GW ;
GBAest =1; GWest = GW;
#Distribuisci i nodi tra i gruppi in parti uguali
groupsindTsDgp = [[distributeAinVecN(Array{Int}(1:GWdgp),N)];
                    [distributeAinVecN(Array{Int}(1:GBAdgp),GWdgp)]]
groupsindTsEst = [[distributeAinVecN(Array{Int}(1:GWest),N)];
                   [distributeAinVecN(Array{Int}(1:GBAest),GWest)]]

#Initialize the variable to store
realUm = DynNets.linSpacedPar(DynNets.fooSdErgm1, N;NgroupsW = GW )[1]
sim_B = Array(linspace(0.4,0.9,GBAdgp+1)[2:end]) #
sim_A = Array(linspace(0.02,0.2,GBAdgp+1)[2:end]) #
sim_W = realUm.*(1-sim_B[groupsindTsDgp[2]])
simParVec = [sim_W,sim_B,sim_A]

@progress "T" for indT = 1:N_Tind
#  indT=1
T = Tvals[indT]
#@sync @parallel
for n = 1:N_est
#n=1
est_times[indT,indN,indGW,n] , est_conv_flag[indT,indN,indGW,n]  = loopEval(T,N,simParVec, useStartVal,groupsindTsDgp,scoreRescType,GBAest,GWest,groupsindTsEst)
est_flag[indT,indN,indGW,n] = true
println(11)
#end
# save
end
@save(save_path,est_conv_flag,est_times,est_flag,Nvals,GWvalStore ,useStartVal)
println([N,GW,T,sum(est_times[indT,indN,indGW,:]),sum(est_conv_flag[indT,indN,indGW,:])/sum(est_flag[indT,indN,indGW,:])])

end

end

##

 ##


text:
module Strided

import Base: parent, size, strides, tail, setindex
using Base: @propagate_inbounds, RangeIndex, Dims
const SliceIndex = Union{RangeIndex,Colon}

using LinearAlgebra

using TupleTools
using TupleTools: StaticLength

export StridedView, @strided, @unsafe_strided, sreshape, sview

const _NTHREADS = Ref(1)
get_num_threads() = _NTHREADS[]

function set_num_threads(n::Int)
    N = Base.Threads.nthreads()
    if n > N
        n = N
        _set_num_threads_warn(n)
    end
    _NTHREADS[] = n
end
@noinline function _set_num_threads_warn(n)
    @warn "Maximal number of threads limited by number of Julia threads,
            setting number of threads equal to Threads.nthreads() = $n"
end

enable_threads() = set_num_threads(Base.Threads.nthreads())
disable_threads() = set_num_threads(1)

const _use_threaded_mul = Ref(false)
use_threaded_mul() = _use_threaded_mul[]

function disable_threaded_mul()
    _use_threaded_mul[] = false
    return
end

function enable_threaded_mul()
    _use_threaded_mul[] = true
    return
end

function __init__()
    set_num_threads(Base.Threads.nthreads())
end

include("abstractstridedview.jl")
include("stridedview.jl")
include("unsafestridedview.jl")
include("mapreduce.jl")
include("broadcast.jl")
include("macros.jl")

# function _precompile_()
#     ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
#     for T in (Float64, ComplexF64)
#         USMatrix = UnsafeStridedView{T, 2, T, typeof(Base.identity)}
#         SMatrix = StridedView{T, 2, Vector{T}, typeof(Base.identity)}
#
#         precompile(Tuple{typeof(mul!), USMatrix, USMatrix, USMatrix})
#         precompile(Tuple{typeof(mul!), USMatrix, USMatrix, USMatrix, T, T})
#
#         precompile(Tuple{typeof(mul!), SMatrix, SMatrix, SMatrix})
#         precompile(Tuple{typeof(mul!), SMatrix, SMatrix, SMatrix, T, T})
#
#         if T <: Complex
#             CUSMatrix = UnsafeStridedView{T, 2, T, typeof(Base.conj)}
#             CSMatrix = StridedView{T, 2, Vector{T}, typeof(Base.conj)}
#
#             precompile(Tuple{typeof(mul!), USMatrix, USMatrix, CUSMatrix})
#             precompile(Tuple{typeof(mul!), USMatrix, CUSMatrix, USMatrix})
#             precompile(Tuple{typeof(mul!), USMatrix, CUSMatrix, CUSMatrix})
#
#             precompile(Tuple{typeof(mul!), SMatrix, SMatrix, CSMatrix})
#             precompile(Tuple{typeof(mul!), SMatrix, CSMatrix, SMatrix})
#             precompile(Tuple{typeof(mul!), SMatrix, CSMatrix, CSMatrix})
#         end
#
#         for N = 1:4
#             USArray = UnsafeStridedView{T, N, T, typeof(Base.identity)}
#             SArray = StridedView{T, N, Vector{T}, typeof(Base.identity)}
#
#             for ASArray in (USArray, SArray)
#                 precompile(Tuple{typeof(copyto!), ASArray, ASArray})
#                 precompile(Tuple{typeof(mul!), ASArray, T, ASArray})
#                 precompile(Tuple{typeof(mul!), ASArray, ASArray, T})
#                 precompile(Tuple{typeof(rmul!), ASArray, T})
#                 precompile(Tuple{typeof(lmul!), T, ASArray})
#                 precompile(Tuple{typeof(axpy!), T, ASArray, ASArray})
#                 precompile(Tuple{typeof(axpby!), T, ASArray, T, ASArray})
#             end
#
#             if T <: Complex
#                 CUSArray = UnsafeStridedView{T, N, T, typeof(Base.conj)}
#                 CSArray = StridedView{T, N, Vector{T}, typeof(Base.conj)}
#
#                 for (ASArray, CASArray) in ((USArray, CUSArray), (SArray, CSArray))
#                     precompile(Tuple{typeof(copyto!), ASArray, CASArray})
#                     precompile(Tuple{typeof(mul!), ASArray, T, CASArray})
#                     precompile(Tuple{typeof(mul!), ASArray, CASArray, T})
#                     precompile(Tuple{typeof(rmul!), CASArray, T})
#                     precompile(Tuple{typeof(lmul!), T, CASArray})
#                     precompile(Tuple{typeof(axpy!), T, CASArray, ASArray})
#                     precompile(Tuple{typeof(axpby!), T, CASArray, T, ASArray})
#                 end
#             end
#         end
#     end
# end
# _precompile_()

end


text:
using LaserTypes
using Test

ω = 0.057
T₀ = 2π/ω
c = 137.036
λ = c*T₀
w₀ = 75 * λ
a₀ = 2.
LG = setup_laser(LaguerreGaussLaser, :atomic, m = 0, p = 0, λ = λ, a₀ = a₀, w₀ = w₀,
    profile=ConstantProfile())
G = setup_laser(GaussLaser, :atomic, λ = λ, a₀ = a₀, w₀ = w₀,
    profile=ConstantProfile())
for i in 1:100
    x, y, z = w₀ * (0.5 .- rand(3))
    t = T₀ * rand()
    @test Fμν([t*c, x, y, z], G) ≈ Fμν([t*c, x, y, z], LG) rtol = 1e-8
end


text:
# TODO: merge gradcheck and check_grads.
# gradcheck iterates over the elements of the first arg.
# check_grads constructs numerical gradient of all args then compares.
# gradcheck has the ability to sample large arrays, check_grads cannot.
# check_grads can handle Tuples and Dicts, gradcheck cannot.
# gradcheck handles non-scalar functions turning them into scalars.

"""

    gradcheck(f, w, x...; kwargs...)

Numerically check the gradient of `f(w,x...;o...)` with respect to its
first argument `w` and return a boolean result.

The argument `w` can be a Number, Array, Tuple or Dict which in turn
can contain other Arrays etc.  Only the largest 10 entries in each
numerical gradient array are checked by default.  If the output of f
is not a number, gradcheck constructs and checks a scalar function by
taking its dot product with a random vector.

# Keywords

* `gcheck=10`: number of largest entries from each numeric array in
  gradient `dw=(grad(f))(w,x...;o...)` compared to their numerical
  estimates.

* `verbose=false`: print detailed messages if true.

* `kwargs=[]`: keyword arguments to be passed to `f`.

* `delta=atol=rtol=cbrt(eps(w))`: tolerance parameters.  See
  `isapprox` for their meaning.

"""
function gradcheck(f, w, x...; kwargs=[], o...)
    y = f(w, x...; kwargs...)
    if !isa(y,Number); f = gc_scalar(f); end
    g = grad(f)
    d = g(w, x...; kwargs...)
    if isa(w, Number)
        gc_number(d, f, w, x...; kwargs=kwargs, o...)
    elseif isbits(eltype(w))
        gc_array(w, d, f, w, x...; kwargs=kwargs, o...)
    else
        k = gc_indices(w)
        pass = true
        for i in k
            pass &= gc_index(w, d, i, f, w, x...; kwargs=kwargs, o...)
        end
        return pass
    end
end

function gc_number(d, f, w, x...; delta=gc_dx(w),rtol=gc_dx(w),atol=gc_dx(w),verbose=false,kwargs=[])
    (w1, w2) = gc_interval(w, delta)
    (f1, f2) = (f(w1,x...;kwargs...), f(w2,x...;kwargs...))
    nd = (f2-f1) / (w2-w1)
    di = (d===nothing ? zero(nd) : d)
    if !isapprox(di, nd; rtol=rtol, atol=atol)
        if verbose; warn("d=$d nd=$nd"); end
        return false
    else
        if verbose && (d*nd!=0); println("gcheck: d=$d nd=$nd"); end
        return true
    end
end

function gc_index(w, d, i, f, w0, x...; o...)
    di = nothing
    try; di = d[i]; end
    if isa(w[i], Number)
        gc_array(w, d, f, w0, x...; icheck=i, o...)
    elseif isbits(eltype(w[i]))
        gc_array(w[i], di, f, w0, x...; o...)
    else
        k = gc_indices(w[i])
        pass = true
        for j in k
            pass &= gc_index(w[i], di, j, f, w0, x...; o...)
        end
        return pass
    end
end

# TODO: handle Tuples, Dict

function gc_array(w, d, f, worig, x...; gcheck=10, icheck=0, kwargs=[],
                  delta=0, atol=0, rtol=0, verbose=false)
    if icheck > 0
        irange = (icheck:icheck)
    elseif length(w) <= gcheck
        irange = (1:length(w))
    else # if d == nothing
        irange = rand(1:length(w), gcheck)
    #else
    #   irange = sortperm(abs(vec(Array(d))),rev=true)[1:gcheck]
    end
    wi = w[irange[1]]
    if delta == 0; delta = gc_dx(wi); end
    if atol == 0; atol = gc_dx(wi); end
    if rtol == 0; rtol = gc_dx(wi); end
    pass = true
    for i in irange
        w0 = w[i]
        (w1, w2) = gc_interval(w0, delta)
        w[i] = w1
        f1 = f(worig, x...; kwargs...)
        w[i] = w2
        f2 = f(worig, x...; kwargs...)
        w[i] = w0
        nd = (f2-f1) / (w2-w1)
        di = (d===nothing ? zero(nd) : d[i])
        if !isapprox(di, nd; rtol=rtol, atol=atol)
            if verbose; warn("d=$di nd=$nd"); end
            pass = false
        else
            if verbose && (di*nd!=0); println("gcheck: d=$di nd=$nd"); end
        end
    end
    return pass
end

gc_dx(x::Number)=cbrt(eps(x))
gc_dx(x)=cbrt(eps(eltype(x)))
gc_indices(w::Tuple)=(1:length(w))
gc_indices(w)=eachindex(w)

function gc_interval(w,d)
    w1=w-d/2
    w2=w+d/2
    (w1 < 0 < w) && (w1=zero(w))
    (w2 > 0 > w) && (w2=zero(w))
    return (w1,w2)
end

function gc_scalar(f)
    # r = MersenneTwister(0)
    function g(x...; o...)
        try
            y = f(x...; o...)
            # v = getval(y)
            # srand(r,1)
            # a = oftype(v, rand(r, size(v)))
            # return sum(y .* a)
            if isa(getval(y), Associative)
                return sumvalues(y)
            else
                return sum(y)  # TODO: revert this back to y.*a once julia6 compat issues resolved?
            end
        catch e
            Base.warn_once("Cannot convert `$f` to a scalar function: $e")
            return 0
        end
    end
    return g
end

### Testing Utilities:

if !isdefined(:addtest)
let tests=[]
    global addtest,runtests,alltests
    alltests()=tests
    addtest(t...)=push!(tests,t)
    function runtests(a=tests)
        for fx in a
            try
                # tx = fixtest(fx...)
                # check_grads(tx...; fname=fx[1]) || throw(:fail)
                f = eval(AutoGrad,fx[1])
                x = fx[2:end]
                gradcheck(f,x...) || throw(:fail)
            catch e
                warn((fx...,"$e"))
            end
        end
    end
end
end

# gradcheck only checks the first arg, this helper will allow us to check all args

applyN(x,f)=f(x...)
addtestN(f,x...)=addtest(:applyN,collect(x),eval(AutoGrad,f))
gradcheckN(f,x...;o...)=gradcheck(applyN,collect(x),f;o...)

# Generate tests based on given ranges

function addtest1(f,r=(-Inf,Inf))          # unary
    bf = broadcast_func(f)
    addtest(f,randin(r))
    addtest(bf,randin(r,2))
end

function addtest2(f,r1=(-Inf,Inf),r2=r1)   # binary
    bf = broadcast_func(f)
    addtestN(f,randin(r1),randin(r2))
    addtestN(bf,randin(r1),randin(r2,2))
    addtestN(bf,randin(r1,2),randin(r2))
    addtestN(bf,randin(r1,2),randin(r2,2))
end

function randin(range, dims...; eps=0.01)
    if isa(range, UnitRange{Int})
        rand(range, dims...)
    elseif range==(-Inf,Inf)
        r = randn(dims...)
        sign_dot(r)*eps + r
    elseif range==(0,Inf)
        eps-log_dot(rand(dims...))
    elseif range==(1,Inf)
        eps+1-log_dot(rand(dims...))
    elseif range==(-1,Inf)
        eps-1-log_dot(rand(dims...))
    elseif range==(-1,1)
        (1-eps)*(2rand(dims...)-1)
    elseif range==(0,1)
        eps+(1-2eps)*rand(dims...)
    elseif range==(0,2)
        eps+2*(1-eps)*rand(dims...)
    elseif range==(-Inf,-1,1,Inf)
        x = sec_dot(randn(dims...))
        sign_dot(x)*eps + x
    else
        error("Unknown range $range")
    end
end

# Alternative gradient check utility -- deprecated.

# EPS, RTOL, ATOL = 1e-4, 1e-4, 1e-6
EPS, RTOL, ATOL = 1e-4, 1e-2, 1e-4

# Check the computed gradients for fun(args) comparing them with numeric
# approximations.  Deprecated, use `gradcheck` instead.

function check_grads(fun, args...; eps=EPS, rtol=RTOL, atol=ATOL, fname=fun)
    #@dbg 2 (:check_grads,fname,:args,args...)
    isempty(args) && error("No args given")
    exact = ntuple(i->grad(fun,i)(args...), length(args))
    numeric = nd(fun, args...; eps=eps)
    #@dbg 2 (:check_grads,fname,:exact,exact,:numeric,numeric)
    same = isequivalent(exact, numeric; rtol=rtol, atol=atol)
    #same || warn((:check_grads,fname,:args,args,:exact,exact,:numeric,numeric))
    return same
end

function nd(f, args...; eps=EPS)
    #@dbg 2 (:nd,f,args..., :eps, eps)
    unary_f = x->f(x...)
    unary_nd(unary_f, args, eps)
end

unary_nd(f, x::Tuple, eps)         = ntuple(i->unary_nd(indexed_function(f, x, i), x[i], eps), length(x))
unary_nd(f, x::Associative, eps)   = (a=similar(x); for(k,v) in x; a[k] = unary_nd(indexed_function(f, x, k), v, eps); end; a)
unary_nd(f, x::AbstractArray, eps) = reshape(eltype(x)[unary_nd(indexed_function(f, x, i), v, eps) for (i,v) in enumerate(x)], size(x))
unary_nd(f, x::Complex, eps)       = ((f(x + eps/2) - f(x - eps/2)) / eps - im*(f(x + im*eps/2) - f(x - im*eps/2)) / eps)
unary_nd(f, x::Real, eps)          = ((f(x + eps/2) - f(x - eps/2)) / eps)

function indexed_function(fun, arg, index)
    function partial_function(x)
        if isa(arg, Tuple)
            local_arg = (arg[1:index-1]..., x, arg[index+1:end]...)
        else
            local_arg = copy(arg); local_arg[index] = x
        end
        return fun(local_arg)
    end
    return partial_function
end

# isequivalent uses isapprox for Number and AbstractArray{T<:Number}
isequivalent(x::Number,y::Number; o...)=isapprox(x,y;o...)
isequivalent{T<:Number,S<:Number}(x::AbstractArray{T},y::AbstractArray{S}; o...)=(size(x)==size(y) && isapprox(x,y;o...))

# isequivalent extends to Tuple, Associative, and other Arrays, comparing elementwise
isequivalent(x::Tuple, y::Tuple; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))
isequivalent(x::AbstractArray, y::AbstractArray; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))
isequivalent(x::Associative, y::Associative; o...)=all(k->isequivalent(get(x,k,nothing),get(y,k,nothing);o...), unique([keys(x)...,keys(y)...]))

# isequivalent treats `nothing` as equivalent to zero or zero array.
isequivalent(x::Number,z::Void; o...)=isequivalent(z,x;o...)
isequivalent{T<:Number}(x::AbstractArray{T},z::Void; o...)=isequivalent(z,x;o...)
isequivalent(z::Void,x::Number; o...)=isapprox(zero(x),x;o...)
isequivalent{T<:Number}(z::Void,x::AbstractArray{T}; rtol::Real=Base.rtoldefault(T), atol::Real=0, norm::Function=vecnorm) = (norm(x) <= atol/(1-rtol)) # Modified from: linalg/generic.jl:522

function fixtest(f, x...)
    f = eval(f)
    y = f(x...)
    # detect and prevent testing of zero / undefined grads
    plist = Any[]               # define fnew(plist)
    alist = Any[x...]           # to return f(alist)
    fargs = Any[]               # call fnew(fargs...)
    for i=1:length(alist)
        gargs = Any[Grad{i},y,y,x...]
        gargs[i+3] = Rec(gargs[i+3])
        g = nothing
        try
            g = f(gargs...)
        catch e
            if isa(e,MethodError) && e.f === f && e.args[1] === Grad{i}
                continue        # warn("No grad $i for $f: $e")
            else
                error("Error during $f$((gargs...)): $e")
            end
        end
        g === nothing && continue # zero grads
        push!(fargs, alist[i])
        alist[i] = Symbol("x$i")
        push!(plist, alist[i])
    end
    isempty(fargs) && error("$f has no differentiable arguments.")
    f1=f; f = eval(Expr(:->, Expr(:tuple, plist...), Expr(:call, f1, alist...)))
    # if f has non-scalar output, sum it
    isbits(y) || (f2=f; f=(x...)->toscalar(f2(x...)))
    return (f,fargs...)
end


text:
using GRIB

function grib_area(file::String)
    GribFile(file) do reader
        m = Message(reader)
        lons, lats = data(m)
        min_lon = minimum(lons)
        max_lon = maximum(lons)
        if min_lon > 180 || max_lon > 180
            min_lon -= 360
            max_lon -= 360
        end
        if min_lon < -180 || max_lon < -180
            min_lon += 360
            max_lon += 360
        end
        return [maximum(lats), min_lon, minimum(lats), max_lon]
    end
end

function get_key_values(file::String, key::String)
    key_values = Vector()
    GribFile(file) do reader
        for msg in reader
            push!(key_values, string(msg[key]))
        end
    end
    return unique(key_values)
end

function get_keys(file::String)::Vector{String}
    keylist = Vector{String}()
    GribFile(file) do reader
        for msg in reader
            for key in keys(msg)
                push!(keylist, key)
            end
        end
    end
    keylist |> unique
end


text:
# # Convex hull of a set of points

# ## Convex hull in a plane

# This examples shows how to find the convex hull in the context of a plane.
# First we have to create an object representing the vertices. We have multiple ways of doing this

using Polyhedra

v = convexhull([0, 0], [0, 1], [1, 0], [0.1, 0.1]) # list of points
v = vrep([[0, 0], [0, 1], [1, 0], [0.1, 0.1]]) # vector of points
# number of points × dimension matrix
x = [0,0,1,0.1]
y = [0,1,0,0.1]
v = vrep([x y])

# Then we can compute the hull of these points using the planar_hull function

Polyhedra.planar_hull(v)

# ## Convex hull in higher dimension

# In higher dimension, we can do it with a linear programming solver implementing the [MathOptInterface](https://github.com/jump-dev/MathOptInterface.jl), e.g.,

import GLPK

removevredundancy(v, GLPK.Optimizer)

# We can also use any Polyhedral library implementing the interface of this package. If we don't specify any library, it falls back to a default one implementing on this package which will use the `planar_hull` if the dimension is 2 (so it's equivalent to the first approach shown above):

p = polyhedron(v)
removevredundancy!(p)
p

# We can also specify a library. For instance, to use `CDDLib`, write `using CDDLib` and then `p = polyhedron(v, CDDLib.Library())`.

max_stars_repo_path: examples/Convex hull of a set of points.jl

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
"""
    all_libs()
"""
function all_libs()
    libs = Symbol[
        :CUBLAS,
        :GemmKernels,
    ]
    return libs
end

function getbenchfunc(lib::Symbol)
    if lib === :CUBLAS
        return benchmark_gemm_cublas!
    elseif lib === :GemmKernels
        return benchmark_gemm_gemmkernels!
    end
    throw(ArgumentError("Unknown library \"$(lib)\""))
end

function benchmark_gemm_cublas!(C, A, B)
    CUDA.CUBLAS.cublasSetMathMode(CUDA.CUBLAS.handle(), CUDA.CUBLAS.CUBLAS_TENSOR_OP_MATH)
    return BenchmarkTools.@benchmark CUDA.@sync begin
        CUDA.CUBLAS.gemmEx!('N', 'N', 1, $A, $B, 0, $C)
    end
end

function benchmark_gemm_gemmkernels!(C, A, B)
    return BenchmarkTools.@benchmark CUDA.@sync begin
        GemmKernels.BLAS.gemmEx!('N', 'N', 1, $A, $B, 0, $C)
    end
end


text:
###################
# Trace interface #
###################

struct CallRecord{T}
    score::Float64
    retval::T
    args::Tuple
end

"""
    get_call_record(trace)::CallRecord
"""
function get_call_record end

"""
    has_choices(trace)::Bool

If the assignment returned by get_assignment is empty or not.
"""
function has_choices end

"""
get_assignment(trace)

Return a value implementing the assignment interface
"""
function get_assignment end

export CallRecord
export get_call_record
export get_assignment
export has_choices

#############
# Generator #
#############

"""
Generator with return value type T and trace type U
"""
abstract type Generator{T,U} end

get_return_type(::Generator{T,U}) where {T,U} = T
get_trace_type(::Generator{T,U}) where {T,U} = U
get_change_type(::Generator) = :Any

"""
Return a boolean indicating whether a gradient of the output is accepted.
"""
function accepts_output_grad end

"""
Returns a tuple of booleans indicating whether a gradient is available, for each argument
"""
function has_argument_grads end

"""
Return a vector of Types indicating the statically known concrete argument types
"""
function get_concrete_argument_types end

"""
    (new_trace, weight, discard, retchange) = fix_update(gen::Generator, new_args, args_change, trace, constraints)
"""
function fix_update end

"""
Optional Generator method; special case of `fix_update` whree addresses cannot be deleted.

Exists for performance optimization.

    (new_trace::U, weight, retchange) = extend(g::Generator{T,U}, args, args_change, trace::U, constraints)
"""
function extend(gen::Generator{T,U}, new_args, args_change, trace::U, constraints) where {T,U}
    (new_trace, weight, discard, retchange) = fix_update(gen, new_args, args_change, trace, constraints)
    if !isempty(discard)
        error("Some addresses were deleted within extend on generator $gen")
    end
    (new_trace, weight, retchange)
end

# TODO add retchange as return value from predict
"""
    (new_trace::U, retchange) = predict(g::Generator{T,U}, args, args_change, trace::U)
"""
function predict(g::Generator, args, args_change, trace)
    (new_trace, weight, _) = extend(g, args, args_change, trace, EmptyAssignment())
    @assert weight == 0.
    new_trace
end

"""
    (trace::U, weight) = generate(g::Generator{T,U}, args, constraints)
"""
function generate end

# TODO will be removed
"""
    trace = simulate(g::Generator, args)
"""
function simulate(gen::Generator, args)
    (trace, weight) = generate(gen, args, EmptyAssignment())
    if weight != 0.
        error("Got nonzero weight during simulate.")
    end
    trace
end

"""
    (trace, discard) = project(g::Generator, args, constraints)
"""
function project end

"""
choices must contain all random choices, (and no extra random choices, or else
it is an error

    trace = assess(g::Generator, args, choices)
"""
function assess end

"""
change of args.
may not simulate.
addresses may be added or deleted.
constraints may collide with existing random choices.

    (new_trace, weight, discard, retchange) = update(g::Generator, new_args, args_change, trace, constraints)
"""
function update end

"""
    (new_trace, weight, retcahnge) = regenerate(
        gen::Generator, new_args, args_change, trace, selection::AddressSet)
"""
function regenerate end

"""
    weight = ungenerate(g::Generator, trace, constraints)
"""
function ungenerate end

"""
    input_grads::Tuple = backprop_params(gen:Generator, trace, retval_grad)
"""
function backprop_params end

"""
    (input_grads::Tuple, values, gradients) = backprop_trace(gen:Generator, trace, selection::AddressSet, retval_grad)
"""
function backprop_trace end

export Generator
export simulate
export extend
export predict
export generate
export project
export assess
export fix_update
export update
export regenerate
export ungenerate
export backprop_params
export backprop_trace

###########################
# incremental computation #
###########################

struct NoChange end
export NoChange

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
#=
interpolant-basis weighted sum-of-squares polynomials (real or real-valued complex)
parametrized by vector of matrices `Ps` derived from interpolant basis and domain constraints

primal cone
w : ∃ Sₗ ⪰ 0 ∀ l, w = ∑ₗ diag(Pₗ Sₗ Pₗ')

dual cone
w : Pₗ' Diagonal(w) Pₗ ⪰ 0 ∀ l
=#

mutable struct WSOSInterpNonnegative{D <: PrimDual, C <: RealCompF} <: Cache
    oa_s::Vector{AE}
    Ps::Vector{Matrix{C}}
    d::Int
    sep_constr::CR
    WSOSInterpNonnegative{D, C}() where {D <: PrimDual, C <: RealCompF} = new{D, C}()
end

function MOIPajarito.Cones.create_cache(
    oa_s::Vector{AE},
    cone::Hypatia.WSOSInterpNonnegativeCone{RealF, C},
    opt::Optimizer,
) where {C <: RealCompF}
    D = primal_or_dual(cone.use_dual)
    cache = WSOSInterpNonnegative{D, C}()
    cache.oa_s = oa_s
    cache.Ps = cone.Ps
    d = cache.d = MOI.dimension(cone)
    @assert length(oa_s) == d
    @assert size(cache.Ps[1], 1) == d
    if D == Prim
        cache.sep_constr = get_sep_constr(cone, opt)
    end
    return cache
end

# TODO maybe define Hypatia MOI cone equal/hash in Hypatia
function hash_cone(cone::Hypatia.WSOSInterpNonnegativeCone{RealF})
    @assert !cone.use_dual
    return hash(cone.Ps)
end

# primal cone functions

function MOIPajarito.Cones.add_init_cuts(cache::WSOSInterpNonnegative{Prim}, opt::Optimizer)
    # add variable bounds wᵢ ≥ 0
    JuMP.@constraint(opt.oa_model, cache.oa_s .>= 0)
    return
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::WSOSInterpNonnegative{Prim},
    opt::Optimizer,
)
    cut = dot_expr(z, cache.oa_s, opt)
    return [cut]
end

# dual cone functions

function MOIPajarito.Cones.add_init_cuts(cache::WSOSInterpNonnegative{Dual}, opt::Optimizer)
    opt.use_init_fixed_oa || return

    # add cuts to make diagonal of each Pₗ' Diagonal(w) Pₗ nonnegative
    # TODO could add SOC cuts or linearizations for 2x2 principal matrix PSD condition
    w = cache.oa_s
    r = zeros(cache.d)
    for P in cache.Ps, P_i in eachcol(P)
        @. r = abs2(P_i)
        JuMP.@constraint(opt.oa_model, JuMP.dot(r, w) >= 0)
    end
    return
end

function MOIPajarito.Cones.get_subp_cuts(
    z::Vector{RealF},
    cache::WSOSInterpNonnegative{Dual},
    opt::Optimizer,
)
    # TODO to get an extreme ray decomposition, use the technique from PY
    # that uses the primal-dual solution from Hypatia to get values for each Sₗ ⪰ 0
    cut = dot_expr(z, cache.oa_s, opt)
    return [cut]
end

function MOIPajarito.Cones.get_sep_cuts(
    s::Vector{RealF},
    cache::WSOSInterpNonnegative{Dual},
    opt::Optimizer,
)
    # extreme ray decomposition adds eigenvector cuts to make
    # each Pₗ' Diagonal(w) Pₗ PSD
    r = zeros(cache.d)
    cuts = AE[]
    for P in cache.Ps
        Λ = P' * Diagonal(s) * P
        F = eigen!(Hermitian(Λ, :U), -Inf, -opt.tol_feas)
        isempty(F.values) && continue

        PV = P * F.vectors
        for PV_i in eachcol(PV)
            @. r = abs2(PV_i)
            cut = dot_expr(r, cache.oa_s, opt)
            push!(cuts, cut)
        end
    end
    return cuts
end


text:
export Hermite,GaussWeight

#TODO: Add general lines
immutable Hermite{T} <: PolynomialSpace{Line{T}}
    L::T
end
Hermite()=Hermite(1.0)

domain(::Hermite)=Line()
canonicalspace(H::Hermite)=Hermite()
spacescompatible(::Hermite,::Hermite)=true #TODO:L
canonicaldomain(H::Hermite)=Line()
tocanonical(H::Hermite,x)=x

#####
# recα/β/γ are given by
#       x p_{n-1} =γ_n p_{n-2} + α_n p_{n-1} +  p_n β_n
#####

recα(::Type,::Hermite,k)=0;recβ(::Type,::Hermite,k)=0.5;recγ(::Type,::Hermite,k)=k-1
recA(::Type,::Hermite,k)=2;recB(::Type,::Hermite,k)=0;recC(::Type,::Hermite,k)=2k

Derivative(H::Hermite,order)=ConcreteDerivative(H,order)

bandinds{H<:Hermite}(D::ConcreteDerivative{H})=0,D.order
rangespace{H<:Hermite}(D::ConcreteDerivative{H})=domainspace(D)
getindex{H<:Hermite}(D::ConcreteDerivative{H},k::Integer,j::Integer) =
        j==k+D.order?one(eltype(D))*2^m*pochhammer(k,D.order):zero(eltype(D))

function hermitep(r::Range,x::Number)
    n=r[end]+1
    T = typeof(x)
    H = Hermite()
    if n≤2
        v=[1.,2x]
    else
        v=Array(promote_type(Float64,typeof(x)),n)  # x may be complex
        v[1]=1.
        v[2]=2x

        for k=2:n-1
            v[k+1]=((x-recα(T,H,k))*v[k] - recγ(T,H,k)*v[k-1])/recβ(T,H,k)
        end
    end
    v[r+1]
end
hermitep(n::Integer,v::Number)=hermitep(n:n,v)[1]
hermitep(n::Range,v::Vector)=transpose(hcat(map(x->hermitep(n,x),v)...))
hermitep(n::Integer,v::Vector)=map(x->hermitep(n,x),v)

identity_fun(sp::Hermite)=Fun([0.,0.5],sp)

# exp(-Lx^2)
immutable GaussWeight{S,T} <: WeightSpace{S,RealBasis,Line{Float64},1}
    space::S
    L::T
end

GaussWeight(H::Hermite)=GaussWeight(H,H.L)
GaussWeight()=GaussWeight(Hermite())

identity_fun(sp::GaussWeight)=identity_fun(sp.space)

spacescompatible(a::GaussWeight,b::GaussWeight)=spacescompatible(a.space,b.space)&&isapprox(a.L,b.L)

function Derivative(sp::GaussWeight,k::Integer)
   if k==1
        x=Multiplication(Fun(identity,sp.space),sp.space)
        D=Derivative(sp.space)
        D2=D-(2sp.L)x
        DerivativeWrapper(SpaceOperator(D2,sp,GaussWeight(rangespace(D2),sp.L)),1)
    else
        D=Derivative(sp)
        DerivativeWrapper(TimesOperator(Derivative(rangespace(D),k-1).op,D.op),k)
    end
end

weight(H::GaussWeight,x)=exp(-H.L*x.^2)

function Base.sum{H<:Hermite,T}(f::Fun{GaussWeight{H,T}})
    @assert space(f).space.L==space(f).L  # only implemented with matching weight
    f.coefficients[1]*sqrt(π)/sqrt(space(f).L)
end

include("hermitetransform.jl")

function Multiplication{H<:Hermite}(f::Fun{H},S::GaussWeight{H})
    M=Multiplication(f,S.space)
    rs=rangespace(M)
    MultiplicationWrapper(f,SpaceOperator(M,S,GaussWeight(rs,rs.L)))
end

function Multiplication{H<:Hermite,T}(f::Fun{GaussWeight{H,T}},S::Hermite)
    M=Multiplication(Fun(f.coefficients,space(f).space),S)
    rs=rangespace(M)
    MultiplicationWrapper(f,SpaceOperator(M,S,GaussWeight(rs,rs.L)))
end


text:
using PersistenceDiagrams
using ScientificTypes

diagram = PersistenceDiagram([(1, Inf), (2, 3)], dim=0)
@test scitype(diagram) == PersistenceDiagram

diagrams = [diagram, diagram, diagram]
@test scitype(diagrams) == AbstractVector{PersistenceDiagram}

true


text:
module JuliaKara
using Blink
include("JuliaKara_noGUI.jl"); using .JuliaKara_noGUI

export
    World,
    place_kara,
    place_tree,
    place_leaf,
    place_mushroom,
    move,
    turnLeft,
    turnRight,
    putLeaf,
    removeLeaf,
    treeLeft,
    treeFront,
    treeRight,
    mushroomFront,
    onLeaf,
    reloadGUI,
    @World,
    save_world,
    load_world,
    get_kara,
    store!,
    reset!,
    world_state_save

import .JuliaKara_noGUI:World,
    place_kara,
    place_tree,
    place_leaf,
    place_mushroom,
    move,
    turnLeft,
    turnRight,
    removeLeaf,
    putLeaf,
    treeLeft,
    treeRight,
    treeFront,
    mushroomFront,
    onLeaf,
    orientation_to_rad,
    Actor,
    get_kara,
    reset!

import .JuliaKara_noGUI.ActorsWorld:Location,
get_actors_at_location_on_layer,
location_move,
location_fix_ooBound,
LocationFullError,
actor_delete!

"""
    World_GUI(world::World,canvas::GtkCanvas,saved_world::World_State,drawing_delay::Float64)

Creates a new World with a GUI component. Contains the actual `world` the `canvas` used for drawing.
A state `saved_world` the world can be reverted to and a `drawing_delay`.
Is used for every GUI communication.
"""
mutable struct World_GUI
    world::World
    window::Window
    saved_world::World_State
    drawing_delay::Float64
    actor_ids::Dict{Actor,Int}
    id_counter::Int
    World_GUI(world,window,saved_world,drawing_delay) = begin
        new(world,window,saved_world,drawing_delay,Dict{Actor,Int}(),0)
    end
end

function get_next_id(w::World_GUI)
    w.id_counter += 1
    return w.id_counter
end

"""
    World(height::Int,width::Int,name::AbstractString)

Creates a new world of size `width` x `height`. `name` is used as a title for the
GTK window.
"""
function World(height::Int,width::Int,name::AbstractString)
    world = World(height,width)
    world_gui = World_GUI(
        world,
        Window(),
        world_state_save(world),
        0,
    )
    Blink.AtomShell.title(world_gui.window,name)
    loadGUI(world_gui)
    return world_gui
end

function reloadGUI(wo::World_GUI)
    wo.window = Window()
    loadGUI(wo)
    drawActors(wo)
    return wo
end

function drawActors(wo::World_GUI)
    @js_ wo.window world.imgs = []
    for ac in wo.world.actors
        if ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:kara]
            draw_kara(wo,ac)
        elseif ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:mushroom]
            draw_mushroom(wo,ac)
        elseif ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:tree]
            draw_tree(wo,ac)
        elseif ac.actor_definition == JuliaKara_noGUI.ACTOR_DEFINITIONS[:leaf]
            draw_leaf(wo,ac)
        else
            error("Missing actor definition, cant draw shape.")
        end
    end
end

function loadGUI(world_gui::World_GUI)
    window = world_gui.window
    wait_until_defined(window,:Blink)
    # Load Vue
    load!(window,joinpath(@__DIR__,"..","res/js/vue.js"),async=false)
    wait_until_defined(window,:Vue)
    # Load Resources
    joinpath(@__DIR__,"res/icons/bugnorth.png")
    Blink.resource(joinpath(@__DIR__,"..","res/icons/load.png"))
    Blink.resource(joinpath(@__DIR__,"..","res/icons/save.png"))
    Blink.resource(joinpath(@__DIR__,"..","res/icons/bugnorth.png"))
    Blink.resource(joinpath(@__DIR__,"..","res/icons/object_tree.png"))
    Blink.resource(joinpath(@__DIR__,"..","res/icons/object_leaf.png"))
    Blink.resource(joinpath(@__DIR__,"..","res/icons/object_mushroom.png"))
    Blink.resource(joinpath(@__DIR__,"..","res/icons/trashcanbig.png"))
    # Load body
    body!(window,read(joinpath(@__DIR__,"..","res/main.html"),String),async=false)
    wait_until_defined(window,:grid)
    blink_set_Grid!(world_gui)
    @js window grid.draw()
    @js window grid.draw()
    # Register handles
    handle(drop_actor_wrapper(world_gui),window,"drop_actor")
    handle(delete_actor_wrapper(world_gui),window,"delete_actor")
    handle(create_actor_wrapper(world_gui),window,"create_actor")
    handle(change_delay_wrapper(world_gui),window,"change_speed")
    handle(load_world_wrapper(world_gui),window,"load")
    handle(save_world_wrapper(world_gui),window,"save")
end

function change_delay_wrapper(world::World_GUI)
    return function (args)
        speed = parse(Int,args)
        world.drawing_delay = 2.0-(2.0/100.0*speed)
    end
end

function create_actor_wrapper(world::World_GUI)
    return function (args)
        typeof = args["type"]
        x = args["x"]
        y = args["y"]

        create_funs = [
            (x,y) -> place_kara(world,x,y,:NORTH)
            (x,y) -> place_tree(world,x,y)
            (x,y) -> place_mushroom(world,x,y)
            (x,y) -> place_leaf(world,x,y)
        ]
        try
            ac = create_funs[typeof](x,y)
        catch e
            if !isa(e,LocationFullError)
                throw(e)
            end
        end
    end
end

function delete_actor_wrapper(world::World_GUI)
    return function (args)
        id = parse(Int,args["id"])
        for a in world.world.actors
            if world.actor_ids[a] == id
                js_remove_actor(world,a)
                delete!(world.actor_ids,a)
                actor_delete!(world.world,a)
                return
            end
        end
        @error "Actor with $id not found!"
    end
end

function drop_actor_wrapper(world::World_GUI)
    return function (args)
        id = parse(Int,args["id"])
        x = args["x"]
        y = args["y"]
        for a in world.world.actors
            if world.actor_ids[a] == id
                try
                    JuliaKara_noGUI.actor_moveto!(
                        world.world,
                        a,
                        JuliaKara_noGUI.Location(
                            x,y
                        )
                    )
                    js_update_actor(world,a)
                catch e
                    if !isa(e,LocationFullError)
                        throw(e)
                    end
                end
                return
            end
        end
        @error "Actor with $id not found!"
    end
end

function load_world_wrapper(wo::World_GUI)
    return function(args)
       path = js(wo.window.shell, Blink.js"""
       var electron = require('electron');
       electron.dialog.showOpenDialog(
        {
            title: "Open Kara world-file",
            properties: ['openFile'],
            filters: [
                       {name: 'world-file', extensions: ['world',]},
                     ]

        })
       """)
        # path = open_dialog("Pick a World-File", b["win_main"], ("*.world",))
        if path != nothing
            wo.world = JuliaKara_noGUI.load_world(path[1])
            wo.saved_world = world_state_save(wo.world)
            drawActors(wo)
            blink_set_Grid!(wo)
        end
    end
end

function save_world_wrapper(wo::World_GUI)
    return function(args)
       path = js(wo.window.shell, Blink.js"""
       var electron = require('electron');
       electron.dialog.showSaveDialog(
        {
            title: "Save Kara world-file",
            filters: [
                       {name: 'World-File', extensions: ['world',]},
                     ]

        })
       """)
        if path != nothing
            save_world(
                wo,
                path
            )
        end
    end
end

function js_isdefined(d::Dict{String,Any})
    length(d) > 0
end

function js_isdefined(d::AbstractString)
    d != "undefined"
end

function wait_until_defined(w::Window,var::Symbol)
    for i in 1:100
        res = @js w typeof($var)
        if js_isdefined(res)
            return
        else
            sleep(0.1)
        end
    end
    @error "Error after loading $var. It is still not defined."
end

function blink_set_Grid!(w::World_GUI)
    width = w.world.size.width
    height = w.world.size.height
    r = @js_ w.window grid.rows = $height
    h = @js_ w.window grid.cols = $width
end

function to_js_direction(direction::Symbol)
    return findfirst(x->x==direction,JuliaKara_noGUI.ActorsWorld.DIRECTIONS)
end

"""
    place_kara(wo::World_GUI,x::Int,y::Int,direction::Symbol=:NORTH)

Places kara in world `wo` at location `x`, `y` in direction `direction`.
`direction` is either of :NORTH, :EAST, :SOUTH, :WEST.
Returns a reference to the created object.

This function is a wrapper around [`JuliaKara_noGUI.place_kara`](@ref) to support GUI.
"""
function place_kara(wo::World_GUI,x::Int,y::Int,direction::Symbol=JuliaKara_noGUI.ActorsWorld.DIRECTIONS[1])
    ac = place_kara(wo.world,x,y,direction)
    draw_kara(wo,ac)
    return ac
end

function draw_kara(wo::World_GUI,ac::Actor)
    newid = get_next_id(wo)
    wo.actor_ids[ac] = newid
    jsdir = to_js_direction(ac.orientation.value)
    x = ac.location.x
    y = ac.location.y
    @js_ wo.window world.place($newid,$x,$y,$jsdir,"bugnorth.png",$(ac.actor_definition.layer))
end

"""
    place_tree(wo::World_GUI,x::Int,y::Int)

Places a tree in world `wo` at location `x`, `y`.
Returns a referenc to the created object.

This function is a wrapper around [`JuliaKara_noGUI.place_tree`](@ref) to support GUI.
"""
function place_tree(wo::World_GUI,x::Int,y::Int)
    ac = place_tree(wo.world,x,y)
    draw_tree(wo,ac)
    return ac
end

function draw_tree(wo::World_GUI,ac::Actor)
    newid = get_next_id(wo)
    wo.actor_ids[ac] = newid
    x = ac.location.x
    y = ac.location.y
    @js_ wo.window world.place($newid,$x,$y,1,"object_tree.png",$(ac.actor_definition.layer))
end

"""
    place_leaf(wo::World_GUI,x::Int,y::Int)

Places a leaf in world `wo` at location `x`, `y`.
Returns a referenc to the created object.

This function is a wrapper around [`JuliaKara_noGUI.place_leaf`](@ref) to support GUI.
"""
function place_leaf(wo::World_GUI,x::Int,y::Int)
    ac = place_leaf(wo.world,x,y)
    draw_leaf(wo,ac)
    return ac
end

function draw_leaf(wo::World_GUI,ac::Actor)
    newid = get_next_id(wo)
    wo.actor_ids[ac] = newid
    x = ac.location.x
    y = ac.location.y
    @js_ wo.window world.place($newid,$x,$y,1,"object_leaf.png",$(ac.actor_definition.layer))
end

"""
    place_mushroom(wo::World_GUI,x::Int,y::Int)

Places a mushroom in world `wo` at location `x`, `y`.
Returns a referenc to the created object.

This function is a wrapper around [`JuliaKara_noGUI.place_mushroom`](@ref) to support GUI.
"""
function place_mushroom(wo::World_GUI,x::Int,y::Int)
    ac = place_mushroom(wo.world,x,y)
    draw_mushroom(wo,ac)
    return ac
end

function draw_mushroom(wo::World_GUI,ac::Actor)
    newid = get_next_id(wo)
    wo.actor_ids[ac] = newid
    x = ac.location.x
    y = ac.location.y
    @js_ wo.window world.place($newid,$x,$y,1,"object_mushroom.png",$(ac.actor_definition.layer))
end

function js_update_actor(wo::World_GUI,ac::Actor)
    @js_ wo.window world.update_actor($(wo.actor_ids[ac]),$(ac.location.x),$(ac.location.y),$(to_js_direction(ac.orientation.value)))
end

function js_remove_actor(wo::World_GUI,ac::Actor)
    @js_ wo.window world.remove_actor($(wo.actor_ids[ac]))
end

function delay(wo::World_GUI)
    if wo.drawing_delay > 0
        sleep(wo.drawing_delay)
    end
end

"""
    move(wo::World_GUI,ac::Actor)

Moves the actor `ac` a step forward in the world `wo`.

This function is a wrapper around [`JuliaKara_noGUI.move`](@ref) to support GUI.
"""
function move(wo::World_GUI,ac::Actor)
    # In case the actor in front of kara is moveable,
    # additionally two fields need to be repainted.
    front_lo = location_move(ac.location,ac.orientation)
    actors_in_front = get_actors_at_location_on_layer(
        wo.world,
        front_lo,
        ac.actor_definition.layer
    )
    move(wo.world,ac)
    js_update_actor(wo,ac)
    js_update_actor.(Ref(wo),actors_in_front)
    delay(wo)
    return nothing
end

"""
    turnLeft(wo::World_GUI,ac::Actor)

Turns the actor `ac` counter clockwise.

This function is a wrapper around [`JuliaKara_noGUI.turnLeft`](@ref) to support GUI.
"""
function turnLeft(wo::World_GUI,ac::Actor)
    turnLeft(wo.world,ac)
    js_update_actor(wo,ac)
    delay(wo)
end

"""
    turnRight(wo::World_GUI,ac::Actor)

Turns the actor `ac` clockwise.

This function is a wrapper around [`JuliaKara_noGUI.turnRight`](@ref) to support GUI.
"""
function turnRight(wo::World_GUI,ac::Actor)
    turnRight(wo.world,ac)
    js_update_actor(wo,ac)
    delay(wo)
end

"""
    removeLeaf(wo::World_GUI,ac::Actor)

Removes an actor of type leaf from the location `ac` is at.

This function is a wrapper around [`JuliaKara_noGUI.removeLeaf`](@ref) to support GUI.
"""
function removeLeaf(wo::World_GUI,ac::Actor)
    layer = JuliaKara_noGUI.ACTOR_DEFINITIONS[:leaf].layer
    leaf = get_actors_at_location_on_layer(wo.world,ac.location,layer)
    if length(leaf) > 0
        js_remove_actor(wo,leaf[1])
        delete!(wo.actor_ids,leaf[1])
    end
    removeLeaf(wo.world,ac)
    delay(wo)
end

"""
    putLeaf(wo::World_GUI,ac::Actor)

Places an actor of type leaf the location `ac` is at.

This function is a wrapper around [`JuliaKara_noGUI.putLeaf`](@ref) to support GUI.
"""
function putLeaf(wo::World_GUI,ac::Actor)
    place_leaf(wo,ac.location.x,ac.location.y)
    delay(wo)
end

"""
    treeLeft(wo::World_GUI,ac::Actor)

Checks if there is an actor of type tree left of actor `ac`.

This function is a wrapper around [`JuliaKara_noGUI.treeLeft`](@ref) to support GUI.
"""
treeLeft(wo::World_GUI,ac::Actor) = treeLeft(wo.world,ac)

"""
    treeRight(wo::World_GUI,ac::Actor)

Checks if there is an actor of type tree right of actor `ac`.

This function is a wrapper around [`JuliaKara_noGUI.treeRight`](@ref) to support GUI.
"""
treeRight(wo::World_GUI,ac::Actor) = treeRight(wo.world,ac)

"""
    treeFront(wo::World_GUI,ac::Actor)

Checks if there is an actor of type tree in front of actor `ac`.

This function is a wrapper around [`JuliaKara_noGUI.treeFront`](@ref) to support GUI.
"""
treeFront(wo::World_GUI,ac::Actor) = treeFront(wo.world,ac)

"""
    mushroomFront(wo::World_GUI,ac::Actor)

Checks if there is an actor of type mushroom in front of actor `ac`.

This function is a wrapper around [`JuliaKara_noGUI.mushroomFront`](@ref) to support GUI.
"""
mushroomFront(wo::World_GUI,ac::Actor) = mushroomFront(wo.world,ac)

"""
    onLeaf(wo::World_GUI,ac::Actor)

Checks if there is an actor of type leaf below of actor `ac`.

This function is a wrapper around [`JuliaKara_noGUI.onLeaf`](@ref) to support GUI.
"""
onLeaf(wo::World_GUI,ac::Actor) = onLeaf(wo.world,ac)

"""
    @World [name] defintion

`definition` is either a `String` describing the path to a world-file which
should be loaded or a `Tuple{Int,Int}` describing the height and the width
of the world

In case a `name` is provided (Must be a name that can be used as a variable name)
a variable in global scope named `name` and a macro named `@name` are created
which allow access to the world (See Examples).

In case no `name` is provided the world is stored in global scope in a variable
named `world` and kara is placed at location 1,1 and refereced with
a global variable named `kara`. Furthermore all function used for interaction
with kara (`move()`, `turnLeft()`, ...) are extended with methods to allow
calls like `move(kara)`.

# Examples
```julia-repl
julia> @World (10,10)
julia> move(kara) # moves kara in world
julia> @world testw (10,10)
julia> lara = @testw place_kara()
julia> @testw move(lara) # moves lara in testw
```
"""
macro World(definition)
    esc(quote
            if typeof($definition) == String
                world = load_world(
                    $definition,
                    "JuliaKara"
                )
            else
                world = World($definition...,"JuliaKara")
                place_kara(world,1,1)
            end
            kara = get_kara(world)
            import JuliaKara.JuliaKara_noGUI:move,
                turnLeft,
                turnRight,
                putLeaf,
               removeLeaf,
                onLeaf,
                treeFront,
                treeLeft,
                treeRight,
                mushroomFront
            function move(ac::JuliaKara.JuliaKara_noGUI.Actor)
                move(world,ac)
            end
            function turnLeft(ac::JuliaKara.JuliaKara_noGUI.Actor)
                turnLeft(world,ac)
            end
            function turnRight(ac::JuliaKara.JuliaKara_noGUI.Actor)
                turnRight(world,ac)
            end
            function putLeaf(ac::JuliaKara.JuliaKara_noGUI.Actor)
                putLeaf(world,ac)
            end
            function removeLeaf(ac::JuliaKara.JuliaKara_noGUI.Actor)
                removeLeaf(world,ac)
            end
            function onLeaf(ac::JuliaKara.JuliaKara_noGUI.Actor)
                onLeaf(world,ac)
            end
            function treeFront(ac::JuliaKara.JuliaKara_noGUI.Actor)
                treeFront(world,ac)
            end
            function treeLeft(ac::JuliaKara.JuliaKara_noGUI.Actor)
                treeLeft(world,ac)
            end
            function treeRight(ac::JuliaKara.JuliaKara_noGUI.Actor)
                treeRight(world,ac)
            end
            function mushroomFront(ac::JuliaKara.JuliaKara_noGUI.Actor)
                mushroomFront(world,ac)
            end
            nothing
        end
        )
end

macro World(name,definition)
    str_name = string(name)
    esc(quote
        if typeof($definition) == String
            $name = load_world(
                $definition,
                $str_name
            )
        else
            $name = World($definition...,$str_name)
        end
        macro $name(command)
            if command.head == :block
                for ca in command.args
                    if ca.head == :call
                        insert!(ca.args,2,$name)
                    end
                end
            elseif command.head == :call
                insert!(command.args,2,$name)
        end
        return command

        end
        nothing
        end)
end

"""
    load_world(path::AbstractString,name::AbstractString)

Loads a world-file from `path` and names the new window `name`.
Creates a new GTK window.
"""
function load_world(path::AbstractString,name::AbstractString)
    loaded_wo = JuliaKara_noGUI.load_world(path)
    wo = World(loaded_wo.size.width,loaded_wo.size.height,name)
    wo.world = loaded_wo
    wo.saved_world = world_state_save(wo.world)
    blink_set_Grid!(wo)
    drawActors(wo)
    return wo
end

"""
    save_world(wo::World_GUI,path::AbstractString)

Saves a world `wo` into a world-file at `path`.
"""
function save_world(wo::World_GUI,path::AbstractString)
    JuliaKara_noGUI.save_world(wo.world,path)
end

get_kara(wo::World_GUI) = get_kara(wo.world)

"""
    store!(wo::World_GUI)

Stores a state of a world `wo` in `wo.saved_world`.
Can be restored by using `reset!(wo::World)`.

# Exaples
```julia-repl
julia> store!(wo)
julia> # do something in wo
julia> reset!(wo)
```
"""
function store!(wo::World_GUI)
    wo.saved_world = world_state_save(wo.world)
end

function reset!(wo::World_GUI,wst::World_State)
    reset!(wo.world,wst)
    drawActors(wo)
    nothing
end

"""
    reset!(wo::World_GUI)

Resets a world `wo` back to a given state `wo.saved_world`.
Can be stored using `store!(wo)`.
Loading a world from a file stores to state at time of loading.

# Examples
```julia-repl
julia> store!(wo)
julia> # Do something in wo
julia> reset!(wo)
```
"""
function reset!(wo::World_GUI)
    reset!(wo,wo.saved_world)
end

end


text:
###################################################
## projectors/networkprojector.jl
## link geographical and network data
###################################################

"""
    `NetworkProjector`: abstract type that creates a link between NetworkData and GeoData
    is used to project geographical data onto network and to use network results
    must implement attributes:
    - `network::Network`
    - `trips::GeoData` : reference to trip data
    methods:
    - `preloadData!`: takes GeoData and load it
    - `getNetworkTrips`: returns network trips associated with given geo-trips
    - `getTripTiming` : returns timing estimation for a geographical trip (random or known)
"""
abstract type NetworkProjector end

function Base.show(io::IO, p::NetworkProjector)
    typeName = split(string(typeof(p)),".")[end]
    nTrips = length(p.trips)
    println(io,"Network Projector: $(typeName)")
    println(io,"$nTrips trips pre-loaded")
end

"""
    `NetworkData` : convenience function to create NetworkData using Projector
"""
function NetworkData(
        proj::NetworkProjector,
        trips::AbstractArray{Int64,1} = eachindex(proj.trips),
        minTimes::AbstractArray{Float64,2} = maxSpeedTimes(proj.network))
    # Create the "NetworkTrip" array
    return NetworkData(proj.network, getNetworkTrips(proj, trips), minTimes)
end


text:
using ParallelRNGs
using Base.Test

rngstate = 1344328
nprocs() < 2 && addprocs(1)

## Test minimal functionality: do we get the same random variables
## from the default RNG with parallel execution as we would get
## sequentially? Note that if nprocs() > 2 the order they're returned
## may be different than the order that they're produced, so we sort
## them.
srand(rngstate)
a = replicate(rand, 4)
srand(rngstate)

@test sort(a) == sort(rand(4))


text:
# Autogenerated wrapper script for MUMPS_seq_MKL_jll for x86_64-w64-mingw32-libgfortran5-cxx03
export libcmumps, libdmumps, libsmumps, libzmumps

using CompilerSupportLibraries_jll
using METIS_jll
using MKL_jll
JLLWrappers.@generate_wrapper_header("MUMPS_seq_MKL")
JLLWrappers.@declare_library_product(libcmumps, "libcmumps.dll")
JLLWrappers.@declare_library_product(libdmumps, "libdmumps.dll")
JLLWrappers.@declare_library_product(libsmumps, "libsmumps.dll")
JLLWrappers.@declare_library_product(libzmumps, "libzmumps.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, METIS_jll, MKL_jll)
    JLLWrappers.@init_library_product(
        libcmumps,
        "bin\\libcmumps.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdmumps,
        "bin\\libdmumps.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsmumps,
        "bin\\libsmumps.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libzmumps,
        "bin\\libzmumps.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
struct Route
    route_short_name::Union{String, Missing}
    route_long_name::Union{String, Missing}
    route_type::UInt16
end