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

# %% [markdown]
# There are a few bugs in the code of https://discourse.julialang.org/t/solving-difference-equation-part-2/67057.
#
# * The code `for k in 2:length(x)` in the function `seqper_new(x; tol=1e-3)` must be corrected as `for k in 2:(length(x) ÷ 2 + 1)`.
# * `nxblock * (ii - 1) + jj` must be `nyblock * (ii - 1) + jj`.
# * `(i - 1) * xpts + j` must be `(i - 1) * ypts + j`
# * The `Threads.@threads for` loop of the function `isoperiodic_test` is not thread-safe, so that the result changes with each execution.  Threads are conflicting at the line `set_parameter!(ds, [par1range[i], par2range[j]])`, where the same `ds` will be modified by multiple threads.
#
# Memory allocations are also unnecessarily too much.
#
# * The size of the array `sol_last` in memory is (NIter + 1) * xpts * ypts * 8 bytes = 2001 * 100 * 100 * 8 bytes = 153 MB, which is quite large, but it is discarded at the end of the function `isoperiodic_test`.
# * Each execution of the `DynamicalSystems.trajectory` function results in a wasted memory allocation, which is repeated tens of thousands of times.
#
# Other points that need improvement
#
# * The calculation results are saved to a file, and the results cannot be reused without going through the file.
# * The code is unnecessarily complicated.
# * You have created the complicated function `isoperiodic_test` to do most of all the work you want to do.
#
# Furthermore, it is not necessary to use `DynamicalSystems.jl` in this case.
#
# __My test implementation:__

# %%
using Base.Threads
using Plots
pyplot() # https://github.com/JuliaPlots/Plots.jl/issues/3560
using CSV
using DataFrames

function f(x, p)
   r, k = p
   a, mu, d = 5.0, 0.5, 0.2
   xnew1 = x[1]*exp(r*(1-x[1]/k)-x[2]/(a+x[1]^2))
   xnew2 = x[2]*exp(mu*x[1]/(a+x[1]^2)-d)
   xnew1, xnew2
end

function traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
    x = x0
    for _ in 1:nwarmup
        x = f(x, p)
    end
    @inbounds tr1[1], tr2[1] = x
    for i in 2:niters+1
        x = f(x, p)
        @inbounds tr1[i], tr2[i] = x
    end
end

function findperiod(x; tol=1e-3)
    n = length(x)
    @inbounds for k in 2:(n ÷ 2 + 1)
        if abs(x[k] - x[1]) ≤ tol
            all(j -> abs(x[j] - x[j-k+1]) ≤ tol, k:n) && return k - 1
        end
    end
    return n
end

function calcperiod2d(f = f, niters = 2000, nwarmup = 50000, x0 = (0.4, 0.5),
        r = (1:0.05:5)[1:end-1], k = (2:0.05:5)[1:end-1], tol=1e-3
    )
    tmp1 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    tmp2 = [Vector{Float64}(undef, niters+1) for _ in 1:nthreads()]
    period = Matrix{Int}(undef, length(r), length(k))
    @threads for j in eachindex(k)
        tid = threadid()
        tr1, tr2 = tmp1[tid], tmp2[tid]
        for i in eachindex(r)
            p = (r[i], k[j])
            traj2d!(tr1, tr2, f, p, x0, niters, nwarmup)
            @inbounds period[i, j] = findperiod(tr1; tol)
        end
    end
    period, r, k
end

xm1ox(x) = (x - 1)/x # x minus 1 over x : {x ≥　1} -> {0 ≤ y < 1}
inv1my(y) = 1/(1 - y) # inverse of 1 minus 1 : {0 ≤ y < 1} -> {x ≥ 1}

function plot_period(r, k, period; kwargs...)
    cbtick = [1, 2, 3, 4, 5, 6, 10, 30]
    colorbar_ticks = (xm1ox.(cbtick), string.(cbtick))
    heatmap(r, k, xm1ox.(period)'; xlabel="r", ylabel="k", title="period",
        colorbar_ticks, kwargs...)
end

function save_period(r, k, period; fn = "result.csv")
    ii, jj = reim(complex.(axes(r, 1), axes(k, 1)'))
    rr, kk = reim(complex.(r, k'))
    df = DataFrame(i = vec(ii), j = vec(jj), r = vec(rr), k = vec(kk), period = vec(period))
    CSV.write(fn, df)
end

function load_period(fn = "result.csv")
    df = CSV.read(fn, DataFrame)
    imax = maximum(df.i)
    jmax = maximum(df.j)
    r = reshape(df.r, imax, jmax)[1:imax]
    k = reshape(df.k, imax, jmax)[1:imax:end]
    period = reshape(df.period, imax, jmax)
    period, r, k
end

# %%
using Logging
disable_logging(Logging.Warn)

# %%
period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
period, r, k = @time calcperiod2d()
plot_period(r, k, period)

# %%
period, r, k = calcperiod2d()
fn = save_period(r, k, period; fn = "result.csv")
df = CSV.read(fn, DataFrame)

# %%
period, r, k = load_period(fn)
plot_period(r, k, period)

# %%
niters = 2000
nwarmup = 50000
x0 = (0.4, 0.5)
r = range(1, 5, length=100+1)[1:end-1]
k = range(2, 5, length=100+1)[1:end-1]
period, = @time calcperiod2d(f, niters, nwarmup, x0, r, k)
period, = @time calcperiod2d(f, niters, nwarmup, x0, r, k)
period, = @time calcperiod2d(f, niters, nwarmup, x0, r, k)
plot_period(r, k, period)

# %%
niters = 2000
nwarmup = 50000
x0 = (0.4, 0.5)
r = range(1, 5, length=200+1)[1:end-1]
k = range(2, 5, length=200+1)[1:end-1]
period, = @time calcperiod2d(f, niters, nwarmup, x0, r, k)
period, = @time calcperiod2d(f, niters, nwarmup, x0, r, k)
period, = @time calcperiod2d(f, niters, nwarmup, x0, r, k)
plot_period(r, k, period)

# %%

max_stars_repo_path: 0018/DynamicalSystem/On Solving difference equation Part 2.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# Test backward compatibility

### Model
nsrc = 2
model, model0, dm = setup_model(tti, viscoacoustic, nlayer)
q, srcGeometry, recGeometry = setup_geom(model; nsrc=nsrc)
dt = srcGeometry.dt[1]
nt = srcGeometry.nt[1]
nrec = length(recGeometry.xloc[1])

@testset "Backward compatibility" begin
    @timeit TIMEROUTPUT "Backward compatibility" begin
        info =  Info(prod(model.n), nt, nsrc)
        @test_logs (:warn, "Info is deprecated and will be removed in future versions") Info(prod(model.n), nt, nsrc)
        @test_logs (:warn, "judiModeling(info::Info, ar...; kw...) is deprecated, use judiModeling(ar...; kw...)") judiModeling(info, model)
        @test_logs (:warn, "judiModeling(info::Info, ar...; kw...) is deprecated, use judiModeling(ar...; kw...)") judiModeling(info, model; options=Options())
        @test_logs (:warn, "judiProjection(info::Info, ar...; kw...) is deprecated, use judiProjection(ar...; kw...)") judiProjection(info, recGeometry)
        @test_logs (:warn, "judiWavefield(info::Info, ar...; kw...) is deprecated, use judiWavefield(ar...; kw...)") judiWavefield(info, dt, nsrc, randn(nt, model.n...))

        @test_throws ArgumentError judiRHS(info, recGeometry, randn(Float32, nt, nrec))
        @test_throws ArgumentError judiLRWF(info, nsrc, randn(nt))

    end
end

ext: jl


text:
module TestOperations

using Test
using MLJBase
using ..Models

@testset "error for operations on nodes" begin
    X = rand(4)
    m = machine(UnivariateStandardizer(), X) |> fit!
    @test_throws ArgumentError inverse_transform(m)
#    @test_deprecated transform(m)
    X = source(rand(4))
    m = machine(UnivariateStandardizer(), X) |> fit!
    @test_throws ArgumentError inverse_transform(m)
#    @test_deprecated transform(m)
end

end

true


text:
"""
    group(iter)
    group(by::Union{Function, Type}, iter)
    group(by::Union{Function, Type}, f::Union{Function, Type}, iter...)

Return a dictionary of groups of elements `x` of the iterable `iter` into groups labeled by
`by(x)`, transforming each element to `f(x)`, where `by` and `f` default to the `identity`
function. If multiple collections (of the same length) are provided, the transformations
`by` and `f` occur elementwise.

# Example

```jldoctest
julia> group(iseven, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Dictionary{Bool,Array{Int64,1}} with 2 entries:
 false │ [1, 3, 5, 7, 9]
 true  │ [2, 4, 6, 8, 10]

julia> group(iseven, x -> x ÷ 2, [1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
Dictionary{Bool,Array{Int64,1}} with 2 entries:
 false │ [0, 1, 2, 3, 4]
 true  │ [1, 2, 3, 4, 5]
```
"""
group(iter) = group(iter, iter)
group(by::Callable, iter) = group(mapview(by, iter), iter)
group(by::Callable, f::Callable, iter) = group(mapview(by, iter), mapview(f, iter))

# TODO multi-input mapview
#group(by::Callable, f::Callable, iter1, iter2) = group(mapview(by, iter1, iter2), mapview(f, iter1, iter2))
#group(by::Callable, f::Callable, iter1, iter2, iters...) = group(mapview(by, iter1, iter2, iters...), mapview(f, iter1, iter2, iters...))
group(by::Callable, f::Callable, iter1, iter2, iters...) = group((x -> by(x...)), (x -> f(x...)), zip(iter1, iter2, iters...))

"""
    group(groups, values)

Return a dictionary of the elements of `values` grouped by the label inidcated by by the
matching element from `groups`.

# Example

```julia
julia> group([true,false,true,false,true], [1,2,3,4,5])
2-element Dictionary{Bool,Array{Int64,1}}
  true │ [1, 3, 5]
 false │ [2, 4]
```
"""
function group(groups, values)
    I = eltype(groups) # TODO EltypeUnknown
    T = eltype(values) # TODO EltypeUnknown

    out = Dictionary{I, Vector{T}}()
    @inbounds for (group, value) in zip(groups, values)
        push!(get!(Vector{T}, out, group), value)
    end

    return out
end

function group(groups::AbstractVector, values::AbstractVector)
    I = eltype(groups)
    T = eltype(values)

    if keys(groups) != keys(values)
        throw(DimensionMismatch("dimensions must match"))
    end

    VT = Core.Compiler.return_type(Base.emptymutable, Tuple{typeof(values)})
    out = Dictionary{I, VT}()
    @inbounds for i in keys(groups)
        group = groups[i]
        value = values[i]
        push!(get!(() -> Base.emptymutable(values), out, group), value)
    end

    return out
end

function group(groups::AbstractDictionary, values::AbstractDictionary)
    I = eltype(groups)
    T = eltype(values)

    out = Dictionary{I, Vector{T}}()

    if sharetokens(groups, values)
        @inbounds for token in tokens(groups)
            group = gettokenvalue(groups, token)
            value = gettokenvalue(values, token)
            push!(get!(Vector{T}, out, group), value)
        end
    else
        if length(groups) != length(values)
            throw(KeyError("Indices must match"))
        end

        for (i, group) in pairs(groups)
            value = values[i]
            push!(get!(Vector{T}, out, group), value)
        end
    end

    return out
end

"""
    groupreduce(by, f, op, iter...; [init])

Applies a `mapreduce`-like operation on the groupings labeled by passing the elements of
`iter` through `by`. Mostly equivalent to `map(g -> reduce(op, g), group(by, f, iter))`,
but designed to be more efficient. If multiple collections (of the same length) are
provided, the transformations `by` and `f` occur elementwise.
"""
groupreduce(by::Callable, op::Callable, iter) = groupreduce(by, identity, op, iter)

groupreduce(by::Callable, f::Callable, op::Callable, iter1, iter2, iters...; kw...) = groupreduce((x -> by(x...)), (x -> f(x...)), op, zip(iter1, iter2, iters...); kw...)

function groupreduce(by::Callable, f::Callable, op::Callable, iter; kw...)
    groupreduce(op, mapview(by, iter), mapview(f, iter); kw...)
end

function groupreduce(op::Callable, groups, values; kw...)
    I = eltype(groups)
    nt = Base.values(kw)
    if nt isa NamedTuple{()}
        T = eltype(values)
    elseif nt isa NamedTuple{(:init,)}
        T = Core.Compiler.return_type(op, Tuple{typeof(nt.init), eltype(values)})
    else
        throw(ArgumentError("groupreduce doesn't support these keyword arguments: $(setdiff(keys(nt), (:init,)))"))
    end

    out = Dictionary{I, T}()
    @inbounds for (group, value) in zip(groups, values)
        (hadtoken, token) = gettoken!(out, group)
        if hadtoken
            settokenvalue!(out, token, op(gettokenvalue(out, token), value))
        else
            if nt isa NamedTuple{()}
                settokenvalue!(out, token, value)
            else
                settokenvalue!(out, token, op(nt.init, value))
            end
        end
    end

    return out
end

"""
    groupsum(by, [f], iter)

Sum the elements of `iter` belonging to different groups, optionally mapping by `f`.
"""
groupsum(iter) = groupsum(identity, iter)
groupsum(by::Callable, iter) = groupsum(by, identity, iter)
groupsum(by::Callable, f::Callable, iter) = groupreduce(by, f, +, iter)

"""
    groupprod(by, [f], iter)

Multiply the elements of `iter` belonging to different groups, optionally mapping by `f`.
"""
groupprod(iter) = groupprod(identity, iter)
groupprod(by::Callable, iter) = groupprod(by, identity, iter)
groupprod(by::Callable, f::Callable, iter) = groupreduce(by, f, *, iter)

"""
    groupcount([by], iter)

Determine the number of elements of `iter` belonging to each group.
"""
groupcount(by::Callable, iter) = groupreduce(by, x->1, +, iter)
groupcount(iter) = groupcount(identity, iter) # A handy extension of `unique`

"""
    groupfirst([by], [f], iter)

Find the first element of each grouping.
"""
groupfirst(iter) = groupreduce(identity, identity, former, iter)
groupfirst(by::Callable, iter) = groupreduce(by, identity, former, iter)
groupfirst(by::Callable, f::Callable, iter) = groupreduce(by, f, former, iter)
groupfirst(groups, values) = groupreduce(former, groups, values)

"""
    grouplast([by], [f], iter)

Find the last element of each grouping.
"""
grouplast(iter) = groupreduce(identity, identity, latter, iter)
grouplast(by::Callable, iter) = groupreduce(by, identity, latter, iter)
grouplast(by::Callable, f::Callable, iter) = groupreduce(by, f, latter, iter)
grouplast(groups, values) = groupreduce(latter, groups, values)

"""
    grouponly([by], [f], iter)

Return a dictionary mapping the unique elements `x` of iter grouped by `by(x)` with value
`f(x)`, where `by` and `f` default to the identity function.

This is an optimized equivalent of `only.(group(by, f, iter))` and is similar to
`Dictionary(by.(iter), f.(iter))`
"""
grouponly(iter) = grouponly(identity, iter)

grouponly(by::Callable, iter) = grouponly(by, identity, iter)

grouponly(by::Callable, f::Callable, iter) = grouponly(mapview(by, iter), mapview(f, iter))

"""
    grouponly(groups, values)

This is an optimized equivalent of `only.(group(groups, values))`.
"""
grouponly(groups, values) = Dictionary(groups, values)
# should make this more generic... if `groups` is an `AbstractIndices` we can use `similar`
# to generate the output...

"""
    groupunique([by], [f], iter)

Return a dictionary mapping the sets of elements `x` of iter grouped by `by(x)` with values
`f(x)`, where `by` and `f` default to the identity function.

This is similar to `unique.(group(by, f, iter))` but each subgroup is an `AbstractIndices`
of distinct values (rather than an `AbstractArray` with possibly repeated values).

# Example

```julia
julia> groupunique(iseven, [1,2,1,2,3])
2-element Dictionary{Bool,Indices{Int64}}
 false │ {2}
  true │ {1, 3}
```
"""
groupunique(iter) = groupunique(identity, iter)
groupunique(by::Callable, iter) = groupunique(by, identity, iter)
groupunique(by::Callable, f::Callable, iter) = groupunique(mapview(by, iter), mapview(f, iter))

"""
    groupunique(groups, values)

Collect the unique elements of `values` by the corresponding element of `groups`, returning
a dictinary of sets.

# Example

```julia
julia> groupunique([true,false,true,false,true], [1,2,1,2,3])
2-element Dictionary{Bool,Indices{Int64}}
  true │ {1, 3}
 false │ {2}
```
"""
function groupunique(groups, values)
    I = eltype(groups)
    T = eltype(values)

    out = Dictionary{I, Indices{T}}()
    for (group, value) in zip(groups, values)
        grouping = get!(Indices{T}, out, group)
        set!(grouping, value)
    end

    return out
end

@deprecate groupinds groupfind

"""
    groupfind([by], container)

Seperate the indices `i` of `container` into groups labeled by `by(container[i])`, where the
default value of `by` is `identity`.

# Example

```jldoctest
julia> groupfind(iseven, [3,4,2,6,5,8])
2-element Dictionary{Bool,Array{Int64,1}}
 false │ [1, 5]
  true │ [2, 3, 4, 6]
```
"""
groupfind(by::Callable, container) = groupfind(mapview(by, container))

function groupfind(container)
    I = eltype(container)
    T = keytype(container)

    out = Dictionary{I, Vector{T}}()
    @inbounds for i in keys(container)
        tmp = get!(Vector{T}, out, container[i])
        push!(tmp, i)
    end
    return out
end

function groupfind(inds::AbstractIndices)
    I = eltype(inds)
    T = keytype(inds)

    out = Dictionary{I, Vector{T}}()
    for i in inds
        tmp = get!(Vector{T}, out, i)
        push!(tmp, i)
    end
    return out
end

function groupfind(dict::AbstractDictionary)
    I = eltype(dict)
    T = keytype(dict)

    inds = keys(dict)
    out = Dictionary{I, Vector{T}}()
    @inbounds if istokenizable(dict)
        for t in tokens(dict)
            tmp = get!(Vector{T}, out, gettokenvalue(dict, t))
            push!(tmp, gettokenvalue(inds, t))
        end
    else
        for i in inds
            tmp = get!(Vector{T}, out, dict[i])
            push!(tmp, i)
        end
    end
    return out
end

# Semi-lazy grouping container (TODO docstring)
struct GroupDictionary{I, T, Parent, Inds <: AbstractDictionary{I}} <: AbstractDictionary{I, T}
    parent::Parent
    inds::Inds
end

Base.parent(g::GroupDictionary) = getfield(g, :parent)
_inds(g::GroupDictionary) = getfield(g, :inds)
Base.keys(g::GroupDictionary) = keys(_inds(g))

Base.isassigned(g::GroupDictionary{I}, k::I) where {I} = isassigned(_inds(g), k)
Base.@propagate_inbounds function Base.getindex(g::GroupDictionary{I}, i::I) where {I}
    inds = _inds(g)[i]
    return view(parent(g), inds)
end

Dictionaries.tokentype(g::GroupDictionary) = tokentype(_inds(g))
@propagate_inbounds Dictionaries.gettoken(g::GroupDictionary{I}, i::I) where {I} = gettoken(_inds(g), i)
Dictionaries.istokenassigned(g::GroupDictionary, t) = istokenassigned(_inds(g), t)
@propagate_inbounds function Dictionaries.gettokenvalue(g::GroupDictionary, t)
    inds = gettokenvalue(_inds(g), t)
    return view(parent(g), inds)
end

"""
    groupview([by], container)

Like `group`, but each grouping is a `view` of the indexable input `container`.

# Example

```jldoctest
julia> v = [3,4,2,6,5,8]
6-element Array{Int64,1}:
 3
 4
 2
 6
 5
 8

julia> groups = groupview(iseven, v)
2-element GroupDictionary{Bool,SubArray{Int64,1,Array{Int64,1},Tuple{Array{Int64,1}},false},Array{Int64,1},Dictionary{Bool,Array{Int64,1}}}
 false │ [3, 5]
  true │ [4, 2, 6, 8]

julia> groups[false] .= 99  # set all the odd values to 99
2-element view(::Array{Int64,1}, [1, 5]) with eltype Int64:
 99
 99

julia> v
6-element Array{Int64,1}:
 99
  4
  2
  6
 99
  8
```
"""
groupview(container) = groupview(identity, container)
function groupview(by::Callable, container)
    inds = groupfind(by, container)
    T = typeof(container)
    V = Core.Compiler.return_type(view, Tuple{T, eltype(inds)})
    return GroupDictionary{keytype(inds), V, T, typeof(inds)}(container, inds)
end


text:
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
# Copyright (c) 2021, Joseph M. Rutherford
# All rights reserved.
#
# Code provided under the license contained in the LICENSE file.
#
# This file depends upon sudoku_common.jl

"""
    as_text(puzzle)

Return the contents of the puzzle as a grid of strings.

Known values are integers greater than 0.
Unknown values are 0.
"""
function as_text_grid(puzzle::SolvablePuzzle)
    source = as_text(puzzle)
    rank = get_rank(puzzle)
    rank_squared = rank*rank
    # rows are each value in a block row separated by a spacer,
    # and an extra spacer for blocks
    # e.g., "||1|2||3|4||" or "||1|2|3||4|5|6||7|8|9||"
    # for rank*rank numbers, rank*rank-1 spacers,rank doubled spacers, and 2 leading spacers =
    # 2*rank*rank + rank + 2 spaces
    column_count =2*rank*rank + rank + 2
    fill_columns = Array{Integer}(undef,(rank_squared,))
    next_column = 1
    for i = 1:rank # column block
        for j = 1:rank # block column
            fill_columns[next_column] = (i-1)*(2*rank+1)+3+(j-1)*2
            next_column += 1
        end
    end

    row_count = rank_squared+rank+1
    spacer = "|"
    horizontal_space = "-"
    interim = fill(spacer,(row_count,column_count))
    next_puzzle_row = 1
    # Now walk across all rows of output
    for i = 1:row_count
        if 1 + (i - 1) % (rank+1) == 1
            for j = 1:column_count
                interim[i,j] = horizontal_space
            end
        else
            puzzle_row = next_puzzle_row
            next_puzzle_row += 1
            # Now walk across all columns of output
            for j = 1:rank_squared # column block
                interim[i,fill_columns[j]] = source[puzzle_row,j]
            end
        end
    end

    # Convert interim set of strings to lines of text
    all_rows = Array{String}(undef,row_count)
    for i = 1:row_count
        # Concatenate columns of strings
        line = ""
        for j = 1:column_count
            line = line*interim[i,j]
        end
        all_rows[i] = line
    end

    # Accumulate lines with newline
    result = ""
    for i = 1:row_count
        result = result*all_rows[i]*"\n"
    end
    return result
end

max_stars_repo_name: joseph-m-rutherford/sudoku-julia

max_issues_repo_name: joseph-m-rutherford/sudoku-julia

max_forks_repo_name: joseph-m-rutherford/sudoku-julia


text:
# Poisson distribution

export Poisson
import Base
using SpecialFunctions: logfactorial

@parameterized Poisson(λ) ≪ CountingMeasure(ℤ[0:∞])

function logdensity(d::Poisson{(:λ,)}, y)
    λ = d.λ
    return y * log(λ) - λ - logfactorial(y)
end

# function logdensity(d::Poisson{(:log_λ,)}, y)
#     log_λ = d.log_λ
#     return y * log(log_λ) + (1 - y) * log(1 - log_λ)
# end

# function density(d::Poisson{(:log_λ,)}, y)
#     log_λ = d.log_λ
#     return 2*log_λ*y - log_λ - y + 1
# end

sampletype(::Poisson) = Int

Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:λ,)}) = rand(rng, T, Dists.Poisson(d.λ))

≪(::Poisson, ::IntegerRange{lo,hi}) where {lo, hi} = lo ≤ 0 && isinf(hi)


text:
module TestStringIndexing

using Test, DataFrames

@testset "iteration" begin
    df = DataFrame(a=1:2, b=3:4)
    er = eachrow(df)
    ec = eachcol(df)
    @test er.a == er."a" == ec.a == ec."a" == df.a == df."a"
    @test_throws ArgumentError er.c
    @test_throws ArgumentError er."c"
    @test_throws ArgumentError ec.c
    @test_throws ArgumentError ec."c"
    @test_throws ArgumentError df.c
    @test_throws ArgumentError df."c"
    @test hasproperty(er, :a) == hasproperty(er, "a") ==
          hasproperty(ec, :a) == hasproperty(ec, "a") ==
          hasproperty(df, :a) == hasproperty(df, "a") == true
    @test hasproperty(er, :c) == hasproperty(er, "c") ==
          hasproperty(ec, :c) == hasproperty(ec, "c") ==
          hasproperty(df, :c) == hasproperty(df, "c") == false

    @test keys(er) == 1:2
    @test propertynames(er) == propertynames(ec) == propertynames(df) ==
          keys(ec) == [:a, :b]
    @test_throws MethodError keys(df)
end

@testset "joins" begin
    df1 = DataFrame(a=1, b=2)
    df2 = DataFrame(a=1, c=2)
    df3 = DataFrame(a=1, d=2)

    # only check if the output is the same in all cases
    for f in (innerjoin, leftjoin, rightjoin, outerjoin, antijoin, semijoin)
        @test f(df1, df2, on=:a) == f(df1, df2, on="a") ==
              f(df1, df2, on=[:a]) == f(df1, df2, on=["a"]) ==
              f(df1, df2, on=:a => :a) == f(df1, df2, on="a" => "a") ==
              f(df1, df2, on=[:a => :a]) == f(df1, df2, on=["a" => "a"])
        @test_throws TypeError f(df1, df2, on = :a => "a")
        @test_throws ArgumentError f(df1, df2, on = [:a => "a"])

        if f === innerjoin || f === outerjoin
            @test f(df1, df2, df3, on=:a) == f(df1, df2, df3, on="a") ==
                  f(df1, df2, df3, on=[:a]) == f(df1, df2, df3, on=["a"]) ==
                  f(df1, df2, df3, on=:a => :a) == f(df1, df2, df3, on="a" => "a") ==
                  f(df1, df2, df3, on=[:a => :a]) == f(df1, df2, df3, on=["a" => "a"])
            @test_throws TypeError f(df1, df2, df3, on = :a => "a")
            @test_throws ArgumentError f(df1, df2, df3, on = [:a => "a"])
        end
    end
end

@testset "reshape" begin
    df = DataFrame(a=repeat([1:3;], inner=[4]),
                   b=repeat([1:4;], inner=[3]),
                   c=1:12, d=1.0:12.0,
                   e=map(string, 'a':'l'))

    # only check if the output is the same in all cases
    for v in (true, false)
        @test stack(df, [:c, :d], [:a], variable_name=:varn, value_name=:valn, view=v) ==
              stack(df, ["c", "d"], ["a"], variable_name="varn", value_name="valn", view=v)
    end

    wide = DataFrame(id=1:12,
                     a=repeat([1:3;], inner=[4]),
                     b=repeat([1:4;], inner=[3]),
                     c=randn(12),
                     d=randn(12))

    long = stack(wide)
    @test unstack(long, :variable, :value) == unstack(long, "variable", "value")
    @test unstack(long, :id, :variable, :value) ==
          unstack(long, "id", "variable", "value")
    @test unstack(long, [:id, :a], :variable, :value) ==
          unstack(long, ["id", "a"], "variable", "value")
    @test unstack(long, :id, :variable, :value, renamecols=x->Symbol(:_, x)) ==
          unstack(long, "id", "variable", "value", renamecols=x->"_"*x)
end

@testset "selection" begin
    df = DataFrame(a=1:2, b=3:4)
    # only check if the output is the same in all cases
    @test select(df, :a, :b => :d, :b => identity => :d2, :b => identity,
                 [:a, :a] => (+), [:a, :a] => (+) => :e, AsTable(:a) => ByRow(first),
                 nrow => :xxx) ==
          select(df, "a", "b" => "d", "b" => identity => "d2", "b" => identity,
                 ["a", "a"] => (+), ["a", "a"] => (+) => "e", AsTable("a") => ByRow(first),
                 nrow => "xxx")
    @test transform(df, :a, :b => :d, :b => identity => :d2, :b => identity,
                    [:a, :a] => (+), [:a, :a] => (+) => :e, AsTable(:a) => ByRow(first),
                    nrow => :xxx) ==
          transform(df, "a", "b" => "d", "b" => identity => "d2", "b" => identity,
                    ["a", "a"] => (+), ["a", "a"] => (+) => "e", AsTable("a") => ByRow(first),
                    nrow => "xxx")
    @test select(df, [:a]) == select(df, ["a"]) ==
          select(df, :a) == select(df, "a")
    @test transform(df, [:a]) == transform(df, ["a"]) ==
          transform(df, :a) == transform(df, "a")

    df2 = copy(df)
    @test select!(df, :a, :b => :d, :b => identity => :d2, :b => identity,
                  [:a, :a] => (+), [:a, :a] => (+) => :e, AsTable(:a) => ByRow(first),
                  nrow => :xxx) ==
          select!(df2, "a", "b" => "d", "b" => identity => "d2", "b" => identity,
                  ["a", "a"] => (+), ["a", "a"] => (+) => "e", AsTable("a") => ByRow(first),
                  nrow => "xxx")

    df = DataFrame(a=1:2, b=3:4)
    df2 = copy(df)
    @test select!(df, [:a]) == select!(df2, ["a"])

    df = DataFrame(a=1:2, b=3:4)
    df2 = copy(df)
    @test select!(df, :a) == select!(df2, "a")

    df = DataFrame(a=1:2, b=3:4)
    df2 = copy(df)
    @test transform!(df, :a, :b => :d, :b => identity => :d2, :b => identity,
                    [:a, :a] => (+), [:a, :a] => (+) => :e, AsTable(:a) => ByRow(first),
                    nrow => :xxx) ==
          transform!(df2, "a", "b" => "d", "b" => identity => "d2", "b" => identity,
                    ["a", "a"] => (+), ["a", "a"] => (+) => "e", AsTable("a") => ByRow(first),
                    nrow => "xxx")

    df = DataFrame(a=1:2, b=3:4)
    df2 = copy(df)
    @test transform!(df, [:a]) == transform!(df2, ["a"])

    df = DataFrame(a=1:2, b=3:4)
    df2 = copy(df)
    @test transform!(df, :a) == transform!(df2, "a")

    df = DataFrame(a=1:2, b=3:4)
    @test_throws ArgumentError select(df, [:a, "b"])
    @test_throws ArgumentError select(df, ["a", :b])
    @test_throws ArgumentError select(df, ["a", "b", "a"])
end

@testset "tables" begin
    df = DataFrame(a=1:2, b=3:4)

    @test columnindex(df, :a) == columnindex(df, "a")
    @test columnindex(df, :c) == columnindex(df, "c")

    @test Tables.schema(df) == Tables.schema(Tables.columntable(df))
    @test [:a, :b] == Tables.columnnames(df) == Tables.columnnames(df[1, :]) ==
          Tables.columnnames(eachrow(df)) == Tables.columnnames(eachcol(df))
end

@testset "split-apply-combine" begin
    df = DataFrame(g=[1, 1, 1, 2, 2], a=1:5)

    # only check if the output is the same in all cases
    gdf = groupby(df, :g)
    @test gdf == groupby(df, "g")
    @test groupby(df, [:g, :a]) == groupby(df, ["g", "a"])

    @test names(gdf) == ["g", "a"]

    k = keys(gdf)
    @test names(k[1]) == ["g"]
    @test propertynames(k[1]) == keys(k[1]) == [:g]
    @test haskey(k[1], :g) == haskey(k[1], "g") == true
    @test haskey(k[1], :a) == haskey(k[1], "a") == false
    @test  k[1].g == k[1]."g" == k[1][:g] == k[1]["g"]

    @test combine(gdf, :a) == combine(gdf, "a") ==
          combine(gdf, [:a]) == combine(gdf, ["a"])

    @test combine(gdf, "a" => identity, ungroup=false) ==
          combine(gdf, :a => identity, ungroup=false)
    @test combine(gdf, ["a"] => identity, ungroup=false) ==
          combine(gdf, [:a] => identity, ungroup=false)
    @test combine(gdf, nrow => :n, ungroup=false) ==
          combine(gdf, nrow => "n", ungroup=false)

    @test combine(gdf, "a" => identity) == combine(gdf, :a => identity)
    @test combine(gdf, ["a"] => identity) == combine(gdf, [:a] => identity)
    @test combine(gdf, nrow => :n) == combine(gdf, nrow => "n")
end

@testset "DataFrameRow" begin
    dfr = DataFrame(a=1:2, b=3:4, c=5:6)[2, ["c", "a"]]
    @test names(dfr) == ["c", "a"]
    @test names(dfr, "a") == names(dfr, :a) == names(dfr, 2) == names(dfr, Not("c")) ==
          names(dfr, ["a"]) == names(dfr, Cols("a")) == names(dfr, Between("a", "a")) == ["a"]
    @test keys(dfr) == propertynames(dfr) == [:c, :a]
    @test haskey(dfr, :a) == haskey(dfr, "a") == true
    @test haskey(dfr, :z) == haskey(dfr, "z") == false
    @test hasproperty(dfr, :a) == hasproperty(dfr, "a") == true
    @test hasproperty(dfr, :z) == hasproperty(dfr, "z") == false
    @test dfr["a"] == dfr[:a] == dfr[2]
    dfr["a"] = 1000
    @test dfr."a" == 1000
    dfr[1:2] = Dict("a" => 100, "c" => 500)
    @test dfr.a == 100
    @test dfr.c == 500
    @test_throws ArgumentError dfr[1:2] = Dict("a" => 100, "f" => 500)
    @test_throws ArgumentError dfr[1:2] = Dict("a" => 100, :c => 500)
end

@testset "names, propertynames and hasproperty" begin
    df = DataFrame(a=1, x1=2, x2=3, x3=4)
    sdf = view(df, 1:1, 1:3)
    dfr = view(df, 1:1, 1:3)
    er = eachrow(df)
    ec = eachcol(df)
    gdf = groupby(df, :a)

    for v in [df, er, ec]
        @test names(v) == ["a", "x1", "x2", "x3"]
        @test propertynames(v) == [:a, :x1, :x2, :x3]
        @test hasproperty(v, :a)
        @test !hasproperty(v, :x)
        @test hasproperty(v, "a")
        @test !hasproperty(v, "x")
        @test names(v, 1) == names(v, :a) == names(v, "a") == ["a"]
        @test names(v, Not(:a)) == names(v, Not("a")) == names(v, Not(1)) ==
              names(v, r"x") == ["x1", "x2", "x3"]
    end

    @test names(gdf) == ["a", "x1", "x2", "x3"]
    @test names(gdf, 1) == names(gdf, :a) == names(gdf, "a") == ["a"]
    @test names(gdf, Not(:a)) == names(gdf, Not("a")) == names(gdf, Not(1)) ==
          names(gdf, r"x") == ["x1", "x2", "x3"]

    for v in [sdf, dfr]
        @test names(v) == ["a", "x1", "x2"]
        @test propertynames(v) == [:a, :x1, :x2]
        @test names(v, 1) == names(v, :a) == names(v, "a") == ["a"]
        @test names(v, Not(:a)) == names(v, Not("a")) == names(v, Not(1)) ==
              names(v, r"x") == ["x1", "x2"]
        @test hasproperty(v, :a)
        @test !hasproperty(v, :x)
        @test hasproperty(v, "a")
        @test !hasproperty(v, "x")
    end
end

@testset "append! and push!" begin
    df = DataFrame(a=1, b=2)
    append!(df, Dict("a" => 2, "b" => 3))
    @test df == DataFrame(a=1:2, b=2:3)

    df = DataFrame()
    append!(df, Dict("a" => 2, "b" => 3))
    @test df == DataFrame(a=2, b=3)

    df = DataFrame(a=1, b=2)
    push!(df, Dict("a" => 2, "b" => 3))
    @test df == DataFrame(a=1:2, b=2:3)
    @test_throws ArgumentError push!(df, Dict("a" => 2, "b" => 3), cols=:orderequal)
    push!(df, Dict("b" => 4, "c" => 0), cols=:union)
    @test isequal(df, DataFrame(a=[1:2; missing], b=2:4, c=[missing, missing, 0]))
end

end # module


text:
# Functions
# =========
#
# Function definitions that are derived from C macros or static inline
# functions.  These are considered as part of htslib's public API but not
# exported from the dynamic library.
#

# utility function
_bitfield(ptr::Ptr, offset::Int) = (unsafe_load(ptr) >> offset) & 1 != 0

# Defined in sam.h
# ----------------

bam_cigar_op(c)     = c & BAM_CIGAR_MASK
bam_cigar_oplen(c)  = c >> BAM_CIGAR_SHIFT
bam_cigar_opchar(c) = _BAM_CIGAR_STR_PADDED[bam_cigar_op(c)+1]
bam_cigar_gen(l, o) = l<<BAM_CIGAR_SHIFT | o
bam_cigar_type(o)   = BAM_CIGAR_TYPE>>(o<<1) & 0x3

bam_is_rev(b::bam1_t)    = (b.core.flag & BAM_FREVERSE) != 0
bam_is_mrev(b::bam1_t)   = (b.core.flag & BAM_FMREVERSE) != 0
bam_get_qname(b::bam1_t) = Cstring(b.data)
bam_get_cigar(b::bam1_t) = Ptr{UInt32}(b.data + b.core.l_qname)
bam_get_seq(b::bam1_t)   = b.data + (b.core.n_cigar<<2) + b.core.l_qname
bam_get_qual(b::bam1_t)  = b.data + (b.core.n_cigar<<2) + b.core.l_qname + (b.core.l_qseq+1)>>1
bam_get_aux(b::bam1_t)   = b.data + (b.core.n_cigar<<2) + b.core.l_qname + (b.core.l_qseq+1)>>1 + b.core.l_qseq
bam_get_l_aux(b::bam1_t) = b.data - (b.core.n_cigar<<2) - b.core.l_qname - b.core.l_qseq - (b.core.l_qseq+1)>>1

bam_seqi(s::Ptr{UInt8}, i::Integer) = unsafe_load(s, i>>1 + 1) >> ((~i&1)<<2) & 0xf

function bam_set_mempolicy(b::Ptr{bam1_t}, policy::Integer)
    dst = Ptr{UInt8}(b) + fieldoffset(bam1_t, 6)
    unsafe_store!(Ptr{UInt32}(dst), convert(UInt32, policy))
    return
end
bam_get_mempolicy(b::bam1_t) = b.mempolicy

# Derive methods for pointer argument.
for func in [:bam_is_rev, :bam_is_mrev, :bam_get_qname, :bam_get_cigar,
             :bam_get_seq, :bam_get_qual, :bam_get_aux, :bam_get_l_aux,
             :bam_get_mempolicy]
    @eval $(func)(b::Ptr{bam1_t}) = $(func)(unsafe_load(b))
end

function sam_itr_next(htsfp::Ptr{htsFile}, itr::Ptr{hts_itr_t}, r::Ptr{bam1_t})
    is_bgzf = _bitfield(Ptr{UInt32}(htsfp), 4)
    is_cram = _bitfield(Ptr{UInt32}(htsfp), 3)
    if !is_bgzf && !is_cram
        # TODO: hts_log_error
        return Cint(-2)
    end
    if itr == C_NULL
        # TODO: hts_log_error
        return Cint(-2)
    end

    multi = _bitfield(Ptr{UInt32}(itr), 4)
    if multi
        return hts_itr_multi_next(htsfp, itr, r)
    else
        return hts_itr_next(is_bgzf ? unsafe_load(htsfp).fp : C_NULL, itr, r, htsfp)
    end
end


text:
#==

sussesive Van der Waals aproximations

idea: helmolhtz equations are costly to evaluate
what if, instead of using a root finder, we aproximate via interpolation via van der waals?

VdW: P(1+a/v2)(v-b) = RT

P = RT/(v-b) - a/v2

procedure: obtaining v₁

relax factor = 0.9 = λ

if aproximation is from liquid:

    using: a/v2 ≈ 0
    b = RT/P₀ - v₀ #estimated covolume, do not test volumes lesser than this
    v₁ = λv₀ + (1-λ)b
    P₁ = P(v₁)
    if P₁ < Pspec:
        P₀,v₀ = P₁,v₁, REPEAT
    else:
        return  P₀,v₀,P₁,v₁

if aproximation is from gas:
    at high volumes:
        P ≈ RT/v - a/v2
        Pv2 = RTv - a
        a = RTv - Pv2 = v(RT-Pv)
        v2 - RT/P v + a/P = 0
        v = 0.5*(RT/P + sqrt((RT/P)^2 - 4a/P))
    v₁ = v₀/λ
    P₁ = P(v₁)
    if P₁ > Pspec:
        P₀,v₀ = P₁,v₁, REPEAT
    else:
        return  P₀,v₀,P₁,v₁

with  P₀,v₀,P₁,v₁, we can start aproximating P

the basis is the following:

1.  P₀ = RT/(v₀-b) - a/v₀^2  # * v₀^2
    P₁ = RT/(v₁-b) - a/v₁^2  #  *v₁^2

2.  P₀v₀^2 = RTv₀^2 /(v₀-b) - a
    P₁v₁^2 = RTv₁^2 /(v₁-b) - a

joining:
3.  P₁v₁^2 - P₀v₀^2   = RTv₁^2 /(v₁-b) - RTv₀^2 /(v₀-b) # * (v₁-b)*(v₀-b)

4. (P₁v₁^2 - P₀v₀^2)*(v₁-b)*(v₀-b) = RTv₁^2(v₀-b) - RTv₀^2(v₁-b)

5. (P₁v₁^2 - P₀v₀^2)*(v₁-b)*(v₀-b) = RTv₁^2(v₀-b) - RTv₀^2(v₁-b)

# M = (P₁v₁^2 - P₀v₀^2)

6. M*(v₁-b)*(v₀-b) = -(RTv₁^2 - RTv₀^2)*b + (RTv₀v₁^2 - RTv₁v₀^2)

# N = RT(v₁^2 - v₀^2), L = RT(v₀v₁^2 - v₁v₀^2)

7. = (v₁-b)*(v₀-b) + (N/M)b - L/M == 0

#expanding:
b^2 (- N/M-v₁-v₀) b + v₁v₀ - L/M

using cuadratic solving: b is obtained:

#

then a is obtaining as follows:

(RT/(v₁-b)-P₁)=  a #using ₁ as is the most recent aproximation

with b and , now we have a local VdW aproximation, we can obtain the new P as the following:

A = a*Pspec/(RT)^2
B = b*Pspec/RT

#generalized cubic compressibility factor polynomial
Z3 + ((c1+c2-1)B -1)Z2 + (c1c2B^2 - (c1+c2)*(B2+B)+A)Z - B(c1c2(B2+B)+A)

#in VdW: c1 = c2 = 0
Z3 + ((-1)B -1)Z2 + (+A)Z - B(+A)
Z3 + -(B +1)Z2 + AZ - AB
s
using cubic solving, a trio of z (zsols) is obtained

if liquid:
    z = minimum(reals(zsols))
if gas
    z = maximum(reals(zsols))

vnew = zRT/Pspec
Pnew = P(vnew)

check convergence, if not:
P₀,v₀,P₁,v₁ = P₁,v₁,Pnew,vnew,  REPEAT

==#
"""
uses van der waals aproximations to calculate the function p(v) = p0

"""
function suva_fzero(p,Pspec,v0,t;
    liquid=true,
    relax=0.9*one(float(Pspec)),
    maxevals=100,
    atol = zero(float(Pspec)),
    rtol = 8*eps(one(float(Pspec))))
    pspec = Pspec

    _Pspec,_v0,_t,_atol,_rtol,_relax = promote(pspec,v0,t,atol,rtol,relax)
    return _suva_fzero(p,_Pspec,_v0,t,liquid,_relax,maxevals,_atol,_rtol)
end

function _suva_fzero(P,pspec::TT,v0::TT,T::TT,liquid=true,relax = TT(0.9),maxevals=100, atol=zero(TT), rtol=8*eps(TT)) where TT
    Pspec = pspec
    R = RGAS #constant defined outside, 8.314...
    v₀ = v0
    P₀= P(v0)

    λ = relax
    _1 = one(TT)
    _0 = zero(TT)
    P₁ = _0
    v₁ = _0
    nan = _0/_0
    RT = R*T
    RT_P = RT/pspec
    invRT = _1/RT
    preA = Pspec*invRT*invRT
    preB = Pspec*invRT
    count = 0

    if liquid
        while count < 110
            # at liquid volumes, a/v2 ≈ 0
            # p ≈ RT/(v-b)
            b = RT/P₀ - v₀ #estimated covolume, do not test volumes lesser than this
            #@show b
            if b < 0 #patological v0, changing to z aprox
                P₀= P(v0)
                b = RT/P₀ - v₀
                #@show v₀,P₀
                #@show b
            end

            v₁ = 1.5*b
            P₁ = P(v₁)
            #@show v₁,P₁

            if (P₁ < Pspec < P₀) | (P₀ < Pspec < P₁) #bounded P, excelent
                break
            elseif (P₀ < P₁< Pspec) #step in the right direction, P1 increases, but is lower than Pspec
               # @show P₁
                P₀,v₀ = P₁,v₁
                count +=1
            elseif (P₁ < P₀ < Pspec) #step in the wrong direction P1 decreased, danger zone
                v₁ = v₀*0.9
                P₁ = P(v₁)
               # @show P₁
            elseif min(Pspec,P₀,P₁) == Pspec #we have increasing pressures, good
                P₀,v₀ = P₁,v₁
                b = RT/P₀ - v₀
                v₁ = RT_P - b
                P₁ = P(v₁)
               # @show P₁
                break
            elseif count == 100
                return error("cannot find initial value for gas")
            else
                P₀,v₀ = P₁,v₁
               # @show P₁
                count +=1
            end
        end
    else
        while count < 110
            #at high volumes, P ≈ RT/v - a/v2
            a = v₀*(RT-P₀*v₀)
           # @show a
            v₁ = 0.5*(RT_P + sqrt(RT_P^2 - 4*a/Pspec))
            P₁ = P(v₁)
            if (P₁ < Pspec < P₀) | (P₀ < Pspec < P₁) #bounded P, excelent
                break
            elseif (Pspec < P₁ < P₀) #step in the right direction, P₁ < P₀, but Pspec < P₁
                P₀,v₀ = P₁,v₁
                count +=1
            elseif (Pspec < P₀ < P₁) #step in the wrong direction, reduce volume directly
                v₁ = v₀*1.1
                P₁ = P(v₁)
            elseif max(Pspec,P₀,P₁) == Pspec #very good, and extra iteration to refine, as (P₀,v₀) could be very bad
                P₀,v₀ = P₁,v₁
                a = v₀*(RT-P₀*v₀)
               # @show a
                v₁ = 0.5*(RT_P + sqrt(RT_P^2 - 4*a/Pspec))
                P₁ = P(v₁)
                break
            elseif count == 100
                return error("cannot find initial value for gas")
            else
                P₀,v₀ = P₁,v₁
                count +=1
            end
        end
    end

    evalscount = 0
    uatol = atol / oneunit(atol) * oneunit(real(v₁))
    adjustunit = oneunit(real(P₀))/oneunit(real(v₀))
    vnew = v₁
    Pnew = P₁

    while evalscount < maxevals

        v₁2 = v₁*v₁
        v₀2 = v₀*v₀
        v₀3 = v₀2*v₀
        v₁3 = v₁2*v₁

        M = (P₁*v₁2 - P₀*v₀2)
        _a = M
        _b = M*(v₀+v₁)
        _c =  RT*(v₀ - v₁) - M*v₁*v₀
        sqrtΔ = sqrt(_b*_b - 4*_a*_c)
        adiv = 0.5/_a
        b1 = (-_b + sqrtΔ)*adiv#covolume
        b2 = (-_b - sqrtΔ)*adiv
        b = b1
        a = (RT/(v₁-b) - P₁)*v₁2 #temperature coefficient

        A = a*preA
        B = b*preB
        #@show -A*B
        #Z3 + ((-1)B -1)Z2 + (+A)Z - B(+A)
        #-AB + AZ +(-B-1)Z2 + Z3
        poly = (-A*B, A, -B-_1, _1)
        zvals = cardan(poly)
        #realvals = cardan2(poly)
        realvals = cardan_reals(zvals)
        #@show zvals
        #@show realvals
        #realvals = filter(>(_0),_realvals)
        z = maximum(realvals)

        vnew = z*RT_P #stops negative volumes
        Pred = RT/(vnew-b) - a/(vnew*vnew)
        #@show Pred
        #@show (a,b)
        if vnew < 0
            vnew= 0.5*(v₀+v₁)
        end
        Pnew = P(vnew)
        ΔP = Pnew-Pspec
        ΔP₀ = P₀-Pspec
        ΔP₁ = P₁-Pspec
        #@show ΔP,ΔP₁,ΔP₀

        if (abs(ΔP) > abs(ΔP₀)> abs(ΔP₁))
            correction = "corrected spec,P₀ is better"
         #   @show correction
            vnew = v₀ - (v₀-v₁)*ΔP₀/(P₀ - P₁)
            #@show vnew
            if vnew < 0 #correction again for negative volumes
                vnew= 0.5*(v₀+v₁)
            end
            Pnew = P(vnew)
            ΔP = Pnew-Pspec
            #@show ΔP,ΔP₁,ΔP₀
        elseif (abs(ΔP) > abs(ΔP₁))
            correction ="corrected spec P₁ is better"
          #  @show correction
            vnew = v₁ - (v₀-v₁)*ΔP₁/(P₀ - P₁)
            if vnew < 0 #correction again for negative volumes
                vnew= 0.5*(v₀+v₁)
            end
            Pnew = P(vnew)
            ΔP = Pnew-Pspec
            #@show ΔP,ΔP₁,ΔP₀
        else
            correction ="no correction"
           # @show correction
        end

        vcor = (vnew-b)
        pcor = (Pnew+a/(vnew*vnew))
        #@show  vcor,pcor
        #@show  vcor*pcor
        #@show  RT
        if ΔP == eps(Pnew)
            return vnew
        end
        #@show ΔP,ΔP₁,ΔP₀
      #  @show Pnew,P₁,P₀
        #progress = abs(ΔP)/abs(ΔP₁)
        #@show progress
        abs(ΔP) <= adjustunit * max(uatol, abs(vnew) * rtol) && return vnew
        iszero(ΔP) && return vnew
        isnan(ΔP) && return nan
        if (ΔP == ΔP₁)
            sdp = sign(ΔP)
            v_sprev = sign(P(prevfloat(vnew))-Pspec)
            v_snext = sign(P(prevfloat(vnew))-Pspec)
            p_sprev = sign(prevfloat(ΔP))
            p_snext = sign(nextfloat(ΔP))
            sdp * v_sprev <= 0 && return vnew
            sdp * v_snext <= 0 && return vnew
            sdp * p_sprev <= 0 && return vnew
            sdp * p_snext <= 0 && return vnew
            return nan
        end

            P₀,v₀,P₁,v₁ = P₁,v₁,Pnew,vnew
        evalscount += 1
        #@show evalscount
    end
    return vnew
end

# P = RT/v-b - a/v2
#  a= v2(RT/v-b - P)
#(p+a/v2)(v-b) > 0

# P (1- aV/RT) = b

#==
v₀2 = v₀^2
v₁2 = v₁^2
v₀3 = v₀2*v₀
v₁3 = v₁2*v₁
M = (P₁*v₁2 - P₀*v₀2)
p₀*v₀2()
_b = (v₁+v₀)
b = -b +1- (1/2M)*sqrt(_b^2 - 4M(k(v₀ - v₁) + P₀*v₀3*v₁ - P₁*v₀*v₁3)) + P₀ v₀3 + P₀ v₀^2 v₁ - p₁ v₀ v₁^2 - p₁ v₁3)/(2 (P₀ v₀^2 - p₁ v₁^2))
==#

# Z = PV/RT


text:
module Patchwork

using FunctionalCollections
using Compat

import Base:
       convert,
       promote_rule,
       isequal,
       ==,
       >>,
       &,
       show,
       <<

export Node,
       Elem,
       properties,
       children,
       haschildren,
       hasproperties,
       withchild,
       withlastchild,
       TextNode,
       text,
       NodeVector,
       Props,
       props,
       attrs,
       EmptyNode,
       MaybeKey,
       tohtml

const MaybeKey = @compat Union{(@compat Void), Symbol}

# A Patchwork node
@compat abstract type Node end

immutable TextNode <: Node
    text::AbstractString
end
text(xs...) =
    TextNode(string(xs...))

convert(::Type{Node}, s::AbstractString) = text(s)
promote_rule(::Type{Node}, ::Type{AbstractString}) = Node

# Abstract out the word "Persistent"
const NodeVector = PersistentVector{Node}
const Props = Dict{Any, Any}

const EmptyNode = NodeVector([])

convert(::Type{NodeVector}, x) =
    NodeVector([convert(Node, y) for y in x])

convert(::Type{NodeVector}, x::Node) =
    NodeVector([x])

convert(::Type{NodeVector}, x::NodeVector) =
    x

convert(::Type{NodeVector}, x::AbstractString) =
    NodeVector([text(x)])

convert(::Type{Props}, x::AbstractArray) = Props(x)

# A DOM Element
immutable Elem{ns, tag} <: Node
    count::Int
    children::NodeVector
    properties::Props
end

function (::Type{Elem{ns, tag}}){ns, tag}(properties, children)
    childvec = convert(NodeVector, children)
    if isempty(properties)
        Elem{ns,tag}(count(childvec), childvec, Props())
    else
        Elem{ns,tag}(count(childvec), childvec, properties)
    end
end

(::Type{Elem{ns, tag}}){ns, tag}() = Elem{ns,tag}(0, EmptyNode)

hasproperties(el::Elem) = isdefined(el, :properties)
haschildren(el::Elem) = !isempty(el.children)
properties(el::Elem) = isdefined(el, :properties) ? el.properties : Props()
children(el::Elem) = el.children

_count(t::TextNode) = 1
_count(el::Elem) = el.count + 1
count(t::TextNode) = 0
count(el::Elem) = el.count
function count(v::NodeVector)
    s = 0
    for x in v
        s +=_count(x)
    end
    s
end

key(n::Elem) = hasproperties(n) ? get(n.properties, :key, nothing) : nothing
key(n::TextNode) = nothing

# A document type
immutable DocVariant{ns}
    elements::Vector{Symbol}
end

# constructors

Elem(ns, name, props, children) =
    Elem{Symbol(ns) , Symbol(name)}(props, children)

Elem(ns::Symbol, name::Symbol, children=EmptyNode; kwargs...) =
    Elem(ns, name, kwargs, children)

Elem(name, children=EmptyNode; kwargs...) =
    Elem(:xhtml, name, kwargs, children)

isequal{ns,name}(a::Elem{ns,name}, b::Elem{ns,name}) =
    a === b || (isequal(properties(a), properties(b)) &&
                isequal(children(a), children(b)))
isequal(a::Elem, b::Elem) = false

==(a::TextNode, b::TextNode) = a.text == b.text
=={ns, name}(a::Elem{ns, name}, b::Elem{ns,name}) =
    a === b || (a.properties == b.properties &&
                a.children == b.children)
==(a::Elem, b::Elem) = false

# Combining elements
(<<){ns, tag}(a::Elem{ns, tag}, b::AbstractArray) =
    Elem{ns, tag}(hasproperties(a) ? a.properties : [], append(a.children, b))
(<<){ns, tag}(a::Elem{ns, tag}, b::Node) =
    Elem{ns, tag}(hasproperties(a) ? a.properties : [], push(a.children, b))

# Manipulating properties
function recmerge(a, b)
    c = Dict{Any, Any}(a)
    for (k, v) in b
        if isa(v, Associative) && haskey(a, k) && isa(a[k], Associative)
            c[k] = recmerge(a[k], v)
        else
            c[k] = b[k]
        end
    end
    c
end

attrs(; kwargs...) = @compat Dict(:attributes => Dict(kwargs))
props(; kwargs...) = kwargs

(&){ns, name}(a::Elem{ns, name}, itr) =
    Elem{ns, name}(hasproperties(a) ?
        recmerge(a.properties, itr) : itr , children(a))

withchild{ns, name}(f::Function, elem::Elem{ns, name}, i::Int) = begin
    cs = children(elem)
    cs′ = assoc(cs, i, f(cs[i]))
    Elem(ns, name, hasproperties(elem) ? elem.properties : [], cs′)
end

withlastchild(f::Function, elem::Elem) =
    withchild(f, elem, length(children(elem)))

include("diff.jl")
include("parse.jl")

include("jsonfmt.jl")
include("hooks.jl")
include("writers.jl")

function showchildren(io, elems, indent_level)
    length(elems) == 0 && return
    write(io, "\n")
    l = length(elems)
    for i=1:l
        show(io, elems[i], indent_level+1)
        i != l && write(io, "\n")
    end
end

function showindent(io, level)
    for i=1:level
        write(io, "  ")
    end
end

@compat function Base.show(io::IO, el::TextNode, indent_level=0)
    showindent(io, indent_level)
    show(io, el.text)
end
function showprops(io, dict)
    write(io, "{")
    write(io, ' ')
    for (k,v) in dict
        print(io, k)
        write(io, '=')
        show(io, v)
        write(io, ' ')
    end
    write(io, "}")
end

@compat function Base.show{ns, tag}(io::IO, el::Elem{ns, tag}, indent_level=0)
    showindent(io, indent_level)
    write(io, "(")
    if namespace(el) != :xhtml
        write(io, namespace(el))
        write(io, ":")
    end
    write(io, tag)
    if hasproperties(el)
        write(io, " ")
        showprops(io, properties(el))
    end
    showchildren(io, children(el), indent_level)
    write(io, ")")
end

function __init__()
    if isdefined(Main, :IJulia)
        include(joinpath(dirname(@__FILE__), "ijulia.jl"))
    end
    try
        load_js_runtime()
    catch
    end

end

end # module


text:
struct Point{T<:Real}
  x::T
  y::T
end
Point(x::T) where {T<:Complex} = Point(real(x), imag(x))
Point(p::Point{T}) where {T<:Real} = Point(p.x, p.y)
function Point(x::AbstractVector{T}) where {T}
  @assert length(x) == 2
  return Point(x[1], x[2])
end

Base.:+(a::Point, b::Point) = Point(a.x + b.x, a.y + b.y)
Base.:-(a::Point, b::Point) = Point(a.x - b.x, a.y - b.y)
Base.:*(a::Point, f::Vector{T}) where {T<:Number} = Point(a.x * f[1], a.y * f[2])
Base.:/(a::Point{}, f::T) where {T<:Number} = Point(a.x / f, a.y / f)

Vector(p::Point) = [p.x, p.y]

middle(a::Point, b::Point) = Point(0.5*(a.x + b.x), 0.5*(a.y + b.y))

angle(A::Point) = atan(A.y, A.x)

function angle(A::Point, B::Point)
  ϕ = angle(B) - angle(A)
  ϕ < -π && return ϕ + 2π
  ϕ > π && return ϕ - 2π
  return ϕ
end

max_stars_repo_name: jwscook/WindingNumbers

max_issues_repo_name: jwscook/WindingNumbers

max_forks_repo_name: jwscook/WindingNumbers


text:
""" [intended for internal use only]
Enforce small numbers to be zero
Helper function used for testing."""
function set_small_to_zero(a::Number)
	# Note: this tolerates up to approximately 1e-8
   	return(a + 1 ≈ 1 ? 0 : a)
end

"""Squared mahalanobis distance for covariance estimator S and center µ"""
function mahalanobis(arrX::AbstractVector{Float64},
					 µ::AbstractVector{Float64}, S::AbstractArray{Float64,2})
    return((arrX - µ)'*inv(S)*(arrX - µ))
end

# Allows the computation to be mapped on rows of a DataFrame
mahalanobis(x::DataFrameRow, µ::AbstractVector{Float64}, S::AbstractArray{Float64,2}) =
    mahalanobis(convert(Vector, x), µ, S)

"""Squared hellinger distance for covariance estimators S and centers µ"""
function hellinger(µ1::AbstractVector{Float64}, S1::AbstractArray{Float64,2},
				   µ2::AbstractVector{Float64}, S2::AbstractArray{Float64,2})
    S = (S1 + S2)/2
    # Singularity sometimes lead to negative determinant values
    # (computational artifacts) and values need to be set to zero
    dS1 = set_small_to_zero(det(S1))
    dS2 = set_small_to_zero(det(S2))
    dS = set_small_to_zero(det(S))
    H = 1 - (dS1^(1/4))*(dS2^(1/4))/(dS^(1/2))*
    	exp((-1/8)*(µ1-µ2)'*inv(S)*(µ1-µ2))
    return(H)
end


text:
# Bug in sparsetocycles fixed in commit
# b03248da3d9f106682e7ecc8b057cfc119fa2613
# in PermPlain
p = PermList([9,8,7,6,5,4,3,2,1])
@test cycles(psparse(cycles(p))) == cycles(p)

# another sparsetocycles bug
p = PermList([3,2,7,9,8,4,1,6,5])
q = PermList([1,4,3,5,9,7,2,8,6])
sp = psparse(p)
sq = psparse(q)
sprod = sp * sq
@test isperm(sprod)
@test cyclelengths(sprod) == cyclelengths(p * q)
#cycles(sprod)


text:
# include("../src/LarSurf.jl")

using LarSurf
import SparseArrays.spzeros
import SparseArrays.dropzeros!
using Plasm, SparseArrays
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

# threshold = 4000
# pth = Pio3d.datasets_join_path("medical/orig/sample-data/nrn4.pklz")
# segmentation = LarSurf.generate_slope([5,6,7])
segmentation = LarSurf.generate_slope([15,16,17])

# data3d = datap["data3d"]
# segmentation = data3d .> threshold
# data_size = LarSurf.size_as_array(size(data3d))
#
# segmentation = zeros(Int8, 5, 6, 7)
# segmentation[2:5,2:5,2:6] .= 1
# Plasm.view(Plasm.numbering(.6)((V,[VV, EV, filteredFV])))

filteredFV, Flin, V, model = LarSurf.get_surface_grid(segmentation)
(VV, EV, FV, CV) = model
Plasm.View((V,[VV, EV, filteredFV]))


text:
"""
# module Force

# Examples

```jldoctest
julia>
```
"""
module Force

using Distances: euclidean

using VCSMD.Geometry

export lennard_jones_scalar,
    lennard_jones_vector

function lennard_jones_scalar(ε::Float64, σ::Float64)::Function
    function (r::Float64)
        x = σ / r
        24 * ε / σ * (2 * x^13 - x^7)
    end
end

function lennard_jones_vector(ε::Float64, σ::Float64)::Function
    f = lennard_jones_scalar(ε, σ)

    function (ri::Point{N}, rj::Point{N}) where {N}
        r = euclidean(rj, ri)
        (rj - ri) / r * f(r)
    end
end

end


text:
using Test
using ReferenceFrameRotations
using LinearAlgebra
using StaticArrays

import Base: isapprox

# Define the function `isapprox` for `EulerAngleAxis` to make the comparison
# easier.
function isapprox(x::EulerAngleAxis, y::EulerAngleAxis; keys...)
    a = isapprox(x.a, y.a; keys...)
    v = isapprox(x.v, y.v; keys...)

    a && v
end

# Define the function `isapprox` for `EulerAngles` to make the comparison
# easier.
function isapprox(x::EulerAngles, y::EulerAngles; keys...)
    a1 = isapprox(x.a1, y.a1; keys...)
    a2 = isapprox(x.a2, y.a2; keys...)
    a3 = isapprox(x.a3, y.a3; keys...)
    r  = x.rot_seq == y.rot_seq

    a1 && a2 && a3 && r
end

# Available rotations.
const rot_seq_array = [:XYX,
                       :XYZ,
                       :XZX,
                       :XZY,
                       :YXY,
                       :YXZ,
                       :YZX,
                       :YZY,
                       :ZXY,
                       :ZXZ,
                       :ZYX,
                       :ZYZ]

# Number of samples.
const samples = 1000

@time @testset "Direction Cosine Matrices" begin
    include("./test_dcm.jl")
end
println("")

@time @testset "Euler Angles" begin
    include("./test_euler_angles.jl")
end
println("")

@time @testset "Euler Angle and Axis" begin
    include("./test_euler_angle_axis.jl")
end
println("")

@time @testset "Quaternions" begin
    include("./test_quaternions.jl")
end
println("")

@time @testset "DCM <=> Euler Angle and Axis" begin
    include("./test_conversion_dcm_euler_angle_axis.jl")
end
println("")

@time @testset "DCM <=> Euler Angles" begin
    include("./test_conversion_dcm_euler_angles.jl")
end
println("")

@time @testset "DCM <=> Quaternions" begin
    include("./test_conversion_dcm_quaternions.jl")
end
println("")

@time @testset "Euler Angles <=> Euler Angle and Axis" begin
    include("./test_conversion_euler_angles_euler_angle_axis.jl")
end
println("")

@time @testset "Euler Angles <=> Euler Angles" begin
    include("./test_conversion_euler_angles_euler_angles.jl")
end
println("")

@time @testset "Euler Angles <=> Quaternion" begin
    include("./test_conversion_euler_angles_quaternions.jl")
end
println("")

@time @testset "Euler Angle Axis <=> Quaternion" begin
    include("./test_conversion_euler_angle_axis_quaternions.jl")
end
println("")

@time @testset "Kinematics using DCMs" begin
    include("./test_kinematics_dcm.jl")
end
println("")

@time @testset "Kinematics using Quaternions" begin
    include("./test_kinematics_quaternions.jl")
end
println("")

@time @testset "Compose rotations" begin
    include("./test_compose_rotations.jl")
end
println("")


text:
# verify(proposition, proof, verifier) -> Bool
# verify(simulator) -> Bool

# The proposition of course can be verified knowing the secret (secret key or randomization factors) alone:
# verify(proposition, secret) -> Bool

# The challnege is last step obtained executing finite state machine. Does not need to contain responses.
# verify(proposition, proof, challenge) -> Bool

# prove(proposition, secret, verifier) -> Simulator

# shuffle(𝔀, g, pk) -> Tuple{Proposition, Secret}
# shuffle(𝔀, g, pk, verifier) -> Simulator

# decrypt(A, g, sk) -> Proposition # secret key is the one which is possible to be used
# decrypt(A, g, sk, verifier) -> Simulator

# Derives a public key for a given secret key
# Proving knowldge of sk such that: x = g^sk
# power(g, sk) -> Proposition
# power(g, sk, verifier) -> Simulator

using CryptoGroups: Group, ElGamal, Enc, Dec, order, modulus

abstract type Proposition end
abstract type Proof end
abstract type Verifier end

function step end
function challenge end

struct Simulator
    proposition::Proposition # Proposition type is the one which selects a type of prover being used
    proof::Proof
    verifier::Verifier
end

struct Shuffle{G <: Group} <: Proposition
    g::G
    pk::G
    𝐞::ElGamal{G}
    𝐞′::ElGamal{G}

    function Shuffle{G}(g::G, pk::G, 𝐞::ElGamal{G}, 𝐞′::ElGamal{G}) where G <: Group
        @assert length(𝐞) == length(𝐞′)
        new(g, pk, 𝐞, 𝐞′)
    end

    Shuffle(g::G, pk::G, 𝐞::ElGamal{G}, 𝐞′::ElGamal{G}) where G <: Group = Shuffle{G}(g, pk, 𝐞, 𝐞′)
end

struct ShuffleSecret
    𝛙::Vector{<:Integer}
    𝐫′::Vector{<:Integer}
end

struct PoSProof{G <: Group} <: Proof
    𝐜::Vector{G}
    𝐜̂::Vector{G}
    t::Tuple{G, G, G, Tuple{G, G}, Vector{G}}
    s::Tuple{BigInt, BigInt, BigInt, BigInt, Vector{BigInt}, Vector{BigInt}}
end

import Base: ==

==(x::PoSProof, y::PoSProof) = x.𝐜 == y.𝐜 && x.𝐜̂ == y.𝐜̂ && x.t == y.t && x.s == y.s

struct PoSChallenge
    𝐡::Vector{<:Group} # Independent set of generators
    𝐮::Vector{BigInt} # PoS commitment challenge
    c::BigInt # Last bit of a challenge
end

###
function verify(proposition::Shuffle, secret::ShuffleSecret)

    (; 𝐞, 𝐞′, g, pk) = proposition
    (; 𝛙, 𝐫′) = secret

    enc = Enc(pk, g)

    return enc(𝐞, 𝐫′)[𝛙] == 𝐞′
end

function verify(proposition::Shuffle, sk::Integer)

    (; 𝐞, 𝐞′, g, pk) = proposition

    @assert g^sk == pk

    dec = Dec(sk)

    return sort(dec(𝐞)) == sort(dec(𝐞′))
end

function gen_shuffle(enc::Enc, e::ElGamal, r::Vector{T}) where T <: Integer

    e_enc = enc(e, r)
    ψ = sortperm(e_enc)

    sort!(e_enc)

    (; g, pk) = enc

    proposition = Shuffle(g, pk, e, e_enc)
    secret = ShuffleSecret(ψ, r)

    return proposition, secret
end

function gen_commitment(g::G, 𝐡::Vector{G}, b::Vector, r::Integer) where G <: Group

    com = g^r * prod(𝐡 .^ b)

    return com
end

function gen_perm_commitment(g::G, 𝐡::Vector{G}, 𝛙::Vector, 𝐫::Vector) where G <: Group

    N = length(𝛙)

    𝐜 = Vector{G}(undef, N)

    for i in 1:N
        j = 𝛙[i]
        𝐜[j] = g^𝐫[j] * 𝐡[i]
    end

    return 𝐜
end

function gen_commitment_chain(g::Group, c0::T, 𝐮::Vector, 𝐫::Vector) where T

    N = length(𝐮)

    𝐜 = Vector{T}(undef, N)

    𝐜[1] = g^𝐫[1] * c0^𝐮[1]

    for i in 2:N
        𝐜[i] = g^𝐫[i] * 𝐜[i-1]^𝐮[i]
    end

    return 𝐜
end

∑(𝐱, q) = mod(sum(𝐱), q) ### Need to improve
∏(𝐱) = prod(𝐱)
∏(f, 𝐱) = prod(f, 𝐱)

using Random: RandomDevice

function gen_roprg(ρ::AbstractVector{UInt8})

    rohash = Hash("sha256")
    prghash = Hash("sha256")
    roprg = ROPRG(ρ, rohash, prghash)

    return roprg
end

gen_roprg() = gen_roprg(rand(RandomDevice(), UInt8, 32))

function prove(proposition::Shuffle{G}, secret::ShuffleSecret, verifier::Verifier; roprg = gen_roprg()) where G <: Group

    (; 𝛙, 𝐫′) = secret
    (; g, pk, 𝐞, 𝐞′) = proposition

    v1 = step(verifier, proposition) # So I could keep a proposition in the coresponding state machine in the end
    𝐡, h = challenge(v1)

    # Would make more sense for length(proposition) == length(secret)
    @assert length(𝛙) == length(𝐞)

    N = length(𝛙)
    q = order(g)

    n = bitlength(q)

    𝐫 = rand(roprg(:𝐫), n, N) # n is part of the sampler here
    𝐫̂ = rand(roprg(:𝐫̂), n, N)
    𝛚 = rand(roprg(:𝛚), n, 4)
    𝛚̂ = rand(roprg(:𝛚̂), n, N)
    𝛚̂′ = rand(roprg(:𝛚̂′), n, N)

    𝐜 = gen_perm_commitment(g, 𝐡, 𝛙, 𝐫)

    v2 = step(v1, 𝐜)
    𝐮 = challenge(v2)

    𝐮′ = 𝐮[𝛙]

    𝐜̂ = gen_commitment_chain(g, h, 𝐮′, 𝐫̂)

    𝐯 = Vector{BigInt}(undef, N)
    𝐯[N] = 1
    for i in N-1:-1:1
        𝐯[i] = 𝐮′[i+1] * 𝐯[i+1]
    end

    r̄ = ∑(𝐫, q)
    r̂ = ∑(𝐫̂ .* 𝐯, q)
    r̃ = ∑(𝐫 .* 𝐮, q)
    r′ = ∑(𝐫′ .* 𝐮, q)

    t₁ = g^𝛚[1]
    t₂ = g^𝛚[2]
    t₃ = g^𝛚[3] * ∏(𝐡 .^ 𝛚̂′)

    enc = Enc(pk, g)
    t₄ = enc(-𝛚[4]) * ∏(𝐞′ .^ 𝛚̂′)

    𝐭̂ = Vector{G}(undef, N)
    𝐭̂[1] = g^𝛚̂[1] * h^𝛚̂′[1]
    for i in 2:N
        𝐭̂[i] = g^𝛚̂[i] * 𝐜̂[i-1]^𝛚̂′[i]
    end

    y = (𝐞, 𝐞′, 𝐜, 𝐜̂, pk)
    t = (t₁, t₂, t₃, t₄, 𝐭̂)

    v3 = step(v2, 𝐜̂, t)
    c = challenge(v3)

    s₁ = mod(𝛚[1] + c * r̄, q)
    s₂ = mod(𝛚[2] + c * r̂, q)
    s₃ = mod(𝛚[3] + c * r̃, q)
    s₄ = mod(𝛚[4] + c * r′, q)

    𝐬̂ = mod.(𝛚̂ .+ c .* 𝐫̂, q) ### What can I do if I have a 0 as one of the elements?
    𝐬′ = mod.(𝛚̂′ .+ c .* 𝐮′, q) ### What to do if 𝐬′ is 0?

    # It would be bad if any s point to 0
    @assert s₁ != 0
    @assert s₂ != 0
    @assert s₃ != 0
    @assert s₄ != 0
    @assert !(0 in 𝐬̂)
    @assert !(0 in 𝐬′)

    s = (s₁, s₂, s₃, s₄, 𝐬̂, 𝐬′) # Do I need to ensure that `s` are without 0 elements

    proof = PoSProof(𝐜, 𝐜̂, t, s)

    simulator = Simulator(proposition, proof, verifier)

    return simulator
end

function verify(proposition::Shuffle, proof::PoSProof, verifier::Verifier)

    v1 = step(verifier, proposition)

    (; 𝐜) = proof
    v2 = step(v1, 𝐜)

    (; 𝐜̂, t) = proof
    v3 = step(v2, 𝐜̂, t)

    chg = PoSChallenge(v3)
    return verify(proposition, proof, chg)
end

function verify(proposition::Shuffle, proof::PoSProof, challenge::PoSChallenge; verbose=false)

    (; g, pk, 𝐞, 𝐞′) = proposition
    (; 𝐜, 𝐜̂, t, s) = proof
    (; 𝐡, 𝐮, c) = challenge
    h = 𝐡[1]

    (s₁, s₂, s₃, s₄, 𝐬̂, 𝐬′) = s
    (t₁, t₂, t₃, t₄, 𝐭̂) = t

    q = order(g)
    N = length(𝐞)

    c̄ = ∏(𝐜) / ∏(𝐡)
    u = mod(∏(𝐮), q)

    ĉ = 𝐜̂[N] / h^u
    c̃ = ∏(𝐜 .^ 𝐮)

    #@infiltrate
    e′ =  ∏(𝐞 .^ 𝐮)

    t₁′ = c̄^(-c) * g^s₁
    t₂′ = ĉ^(-c) * g^s₂
    t₃′ = c̃^(-c) * g^s₃ * ∏(𝐡 .^ 𝐬′) # 𝐬′ is 0!

    enc = Enc(pk, g)
    t₄′ = e′^(-c) * enc(-s₄) * ∏(𝐞′ .^ 𝐬′)

    𝐭̂′ = Vector(undef, N)

    𝐭̂′[1] = 𝐜̂[1]^(-c) * g^𝐬̂[1] * h^𝐬′[1]

    for i in 2:N
        𝐭̂′[i] = 𝐜̂[i]^(-c) * g^𝐬̂[i] * 𝐜̂[i-1]^𝐬′[i]
    end

    report = Report()

    report &= "t₁", t₁ == t₁′
    report &= "t₂", t₂ == t₂′
    report &= "t₃", t₃ == t₃′
    report &= "t₄", t₄ == t₄′

    report &= "𝐭̂", 𝐭̂ .== 𝐭̂′

    if verbose || isvalid(report) == false
        println(report)
    end

    return isvalid(report)
end

verify(simulator::Simulator) = verify(simulator.proposition, simulator.proof, simulator.verifier)

function shuffle(𝐞::ElGamal{G}, g::G, pk::G; roprg = gen_roprg()) where G <: Group

    # Need to abstract this into a function argument
    q = order(g)
    N = length(𝐞)

    n = bitlength(q)

    𝐫′ = rand(roprg(:𝐫′), n, N)

    enc = Enc(pk, g)

    return gen_shuffle(enc, 𝐞, 𝐫′) # I may also refactor it as shuffle.
end

shuffle(𝐞::ElGamal{G}, enc::Enc; roprg = gen_roprg()) where G <: Group = shuffle(𝐞, enc.g, enc.pk; roprg)

function shuffle(𝐞::ElGamal{G}, g::G, pk::G, verifier::Verifier; roprg = gen_roprg()) where G <: Group
    proposition, secret = shuffle(𝐞, g, pk; roprg)
    return prove(proposition, secret, verifier; roprg)
end

shuffle(𝐞::ElGamal{G}, enc::Enc, verifier::Verifier; roprg = gen_roprg()) where G <: Group = shuffle(𝐞, enc.g, enc.pk, verifier; roprg)


text:
# This file includes the simulation of three node network simulation

import Cairo, Fontconfig
using GraphPlot
using Graphs
using SimpleWeightedGraphs
using LinearAlgebra
using DynamicalNetworks
using DifferentialEquations
using PyPlot
pygui(true)
close("all")

# Helper functions
getwave(x, i, d) = getindex.(x, (i - 1) * d + 1)

# Node settings
nodetype = Lorenz           # Node type
P = diagm([1.0, 0.0, 0])    # State coupling matrix
ηt = 10.0                   # Threshold for maximum Lyapuno exponents

# Topology
n = 3
graph = SimpleWeightedGraph(cycle_graph(n))
add_edge!(graph, 2, 3, 5)
gplot(graph, nodelabel = 1:nv(graph), layout = circular_layout)

# Coupling strength
W = laplacian_matrix(graph)
Q = [
    1 0 0
    0 0 1
    0 1 0
]
M = eigvecs(Q)
G = inv(M) * (W ./ diag(W)) * M
B = G[1:1, 1:1]
D = G[2:3, 2:3]
λt = eigvals(B)
λs = eigvals(D)
λt[1] > λs[2] || error("Synchronization is not possible")
ϵr = ηt ./ [λt[1], λs[2]]
ϵ = ϵr[1] + abs(ϵr[2] - ϵr[1]) * 0.5

# Callback
nbits = 10
tbit = 100.0
bits = rand(Bool, nbits)
ti = 0.0
tf = nbits * tbit
tswitch = ti:tbit:tf-tbit
srcnode1, dstnode1 = 2, 3     # Control link nodes
srcnode2, dstnode2 = 1, 2     # Control link nodes
high = 5.0          # Control coupling strength
low = 1.0           # Control coupling strength
condition(u, t, integrator) = integrator.t ∈ tswitch
function action!(integrator)
    net, bits = integrator.sol.prob.p
    if popfirst!(bits)
        add_edge!(net.graph, srcnode1, dstnode1, high)
        add_edge!(net.graph, srcnode2, dstnode2, low)
    else
        add_edge!(net.graph, srcnode1, dstnode1, low)
        add_edge!(net.graph, srcnode2, dstnode2, high)
    end
end
callback = DiscreteCallback(condition, action!)

# Network
σ = 0.0  # Noise strength
H = getnoisemat(graph) * σ
net = SDENetwork(nodetype, ϵ, graph, H, P)
sol = solvenet!(net, (ti, tf), (copy(bits),), solverkwargs = (saveat = 0.001,), probkwargs = (callback = callback, tstops = tswitch))
t = sol.t
x = sol.u

# Plots
d = dimension(nodetype)
subplot(221)
plot(t, abs.(getwave(x, 1, d) - getwave(x, 2, d)))
axvline.(tswitch, linestyle = "dashed", color = "red")
xlabel("Time")
ylabel("e12")
title(join(string.(Int.(bits))))
subplot(222)
plot(t, abs.(getwave(x, 2, d) - getwave(x, 3, d)))
axvline.(tswitch, linestyle = "dashed", color = "red")
xlabel("Time")
ylabel("e23")
title(join(string.(Int.(bits))))
subplot(223)
plot(t, getwave(x, 1, d))
xlabel("Time")
ylabel("x1")
subplot(224)
plot(getindex.(x, 1), getindex.(x, 2))
xlabel("x1")
ylabel("x2")
tight_layout()

figure()
plot(λt, zeros(length(λt)), marker = "*", label = L"\lambda_t")
plot(λs, ones(length(λs)), marker = "o", label = L"\lambda_s")
axvline.([λt[1], λs[2]], linestyle = "dashed", color = "red")
axvline(ηt / ϵ, linestyle = "dashed", color = "green")
xlabel(L"\lambda")
title("Eigenvalue distribution")
legend()


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
# Date: September 2020
# Copyright: Georgia Institute of Technology, 2020

export loss_supervised

function loss_supervised(
    Net::NetworkConditionalHINT,
    X::AbstractArray{Float32,4},
    Y::AbstractArray{Float32,4}
)

    Zx, Zy, logdet = Net.forward(X, Y)
    CUDA.reclaim()
    z_size = size(Zx)

    f = sum(logpdf(0f0, 1f0, Zx))
    f = f + sum(logpdf(0f0, 1f0, Zy))
    f = f + logdet*z_size[4]

    ΔZx = -gradlogpdf(0f0, 1f0, Zx)/z_size[4]
    ΔZy = -gradlogpdf(0f0, 1f0, Zy)/z_size[4]

    ΔX, ΔY = Net.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
    CUDA.reclaim()
    GC.gc()

    return -f/z_size[4], ΔX, ΔY
end


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

    StandardCurveResult.jl

    Author: Tom Price
    Date: Dec 2018

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

abstract type StandardCurveResult end

struct TargetResultEle <: StandardCurveResult
    target_id   ::Int_T
    slope       ::Float_T
    offset      ::Float_T
    efficiency  ::Float_T
    r2          ::Float_T
end
const EMPTY_TRE = TargetResultEle(0, fill(NaN_T, 4)...)

struct GroupResultEle <: StandardCurveResult
    well        ::Vector{Int_T}
    target_id   ::Int_T
    cq_mean     ::Float_T
    cq_sd       ::Float_T
    qty_mean    ::Float_T
    qty_sd      ::Float_T
end
const EMPTY_GRE = GroupResultEle([], 0, fill(NaN_T, 4)...)

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
## Юнит тест проверяет корректность решения прямой задачи.
## Алгоритм описан в /docs/src/direct/adjoint_check.md
## Возвращает решение, аналитическое решение, сетку по X, сетку по T.

using NonLinearReactionAdvectionDiffusionWithFrontData;
using NonLinearReactionAdvectionDiffusionWithFrontData: phidetermination, Φ;
using NonLinearReactionAdvectionDiffusionWithFrontData: f1, f2;
using NonLinearReactionAdvectionDiffusionWithFrontData: heterogeneity, adjointRP, ∂ARP_∂y;

using ForwardDiff;

## Сначала, сгенирируем экспериментальные данные, на увеличенном числе узлов.
u_l(t) = -8
u_r(t) =  4
qf(x) = 4*sin(3 * π * x);       # Коэффициент линейного усиления, который в обратной
                                ## задаче необходимо определить, но при генерации априорной
                                ## информации мы задаем некоторый коэффициент, который,
                                ## собственно, после имея априорную информацию и будем определять.
ε = 0.2;                        # Малый параметр при старшей производной
a, b = 0, 1;                    # Область по X
t₀, T = 0.0, 0.28;              # Область по T
N, M = 250, 500;                # Увеличенное Кол-во разбиений по X, T
h = (b-a)/N;                    # шаг по X
τ = (T-t₀)/M;                   # шаг по T
Xₙ = [a  + n*h for n in 0:N];   # Сетка по Х
Tₘ = [t₀ + m*τ for m in 0:M];   # Сетка по Т
qₙ =      qf.(Xₙ);              # Сеточные значения коэффициента линейного усиления
ulₘ=    u_l.(Tₘ);               # Сеточные значения левого  ГУ
urₘ=    u_r.(Tₘ);               # Сеточные значения правого ГУ
u₀ = u_init.(Xₙ);               # Начальные условия
##
u, XX, TP = solve(u₀, Xₙ, N, Tₘ, M, ε, ulₘ, urₘ, qₙ);

## Генерация априорной информации
ϕl      = phidetermination(qₙ, ulₘ, Xₙ, N, Tₘ, M);                  # Левый вырожденный корень
ϕr      = phidetermination(qₙ, urₘ, Xₙ, N, Tₘ, M, reverseX = true); # Правый вырожденный корень
ϕ       = Φ(ϕl, ϕr, N, M);                                          # Серединный корень
f1_data = f1(ϕ, u, Xₙ, N, M);                                       # Положение переходного слоя
f2_data = f2(f1_data, u, Xₙ, N, M);                                 # Значение функции на переходном слое

## Решение сопряженной задачи

## Теперь переходим непосредственно к проверке
## Модельная функция
raw"""
    g(n::Int, X = Xₙ, N::Int,
        m::Int, T = Tₘ, M::Int) -> Real

Модельная функция ``(1-2t)\sin(\pi x)``.
`n` — номер узла в сетке по X. `m` — номер шага в сетке по T.
`X` — сетка по X, `N` — кол-во **интервалов** сетки.
`Tₘ` — сетка по T, `M` — кол-во **интервалов** сетки.
Захватывает переменную `T` конечное время моделирования из окружения.
"""
function g(n::Int, X, N::Int,
            m::Int, Tₘ, M::Int)
    t = Tₘ[m]
    x = X[n]
    return (1 - 2t/T)*sin(π*x)
end

## Выберем априорный параметр аппроксимации дельта-функции
ω = 0.0001;
raw"""
    g_d(n::Int, Xₙ::Vector, N::Int,
        m::Int, Tₘ::Vector, M::Int,
        ε::Real, qₙ::Vector,
        u::Matrix, f1::Vector, f2::Vector) -> Real

Вычисляет невязку, т.е. рез ультат подстановки ``g`` в постановку сопряженной задачи.
`n` — номер узла в сетке по X. `m` — номер шага в сетке по T.
`X` — сетка по X, `N` — кол-во **интервалов** сетки.
`Tₘ` — сетка по T, `M` — кол-во **интервалов** сетки.
"""
function g_d(n::Int, Xₙ::Vector, N::Int,
                m::Int, Tₘ::Vector, M::Int,
                ε::Real, qₙ::Vector,
                u::Matrix, f1::Vector, f2::Vector,
               w::Real)
    x = Xₙ[n];
    t = Tₘ[m];
    out  = 2/T * sin(π * x) - ( - ε * π^2 * (1 - 2t/T) * sin(π * x)) + π * (1 - 2t/T) * cos(π * x) * u[n,m] + qₙ[n] * (1 - 2t/T) * sin(π * x) - heterogeneity(n, m, Xₙ[2:N], N, Tₘ, M, u, f1, f2, w)
    return out;
end

Uₙₘ = u;                                        # Сохраним матрицу решения прямой задачи
ψl = [0.0 for i in 1:M+1];                      # Левые  ГУ для модельной функции
ψr = [0.0 for i in 1:M+1];                      # Правые ГУ для модельной функции
ψ₀ = [g(n, Xₙ, N, M+1, Tₘ, M) for n in 1:N+1];  # Начальные условия для модельной функции
## Внимание! сопряженная задача — ретроспективная
## ``u₀ = g(x, T)``
## Модельное решение найденное с помощью известного аналитического решения
## Для его генерации используем инвертированную сетку по времени, чтобы
## массив `ψ_model` и `ψ` имели одно и тоже направление хода времени.
ψ_model = [ g(n, Xₙ, N, m, Tₘ, M) for n in 1:N+1, m in M+1:-1:1];

## Создадим функцию, которая будет вычислять вектор правой части с добавлением невязки
raw"""
    RP(y, m, Xₙ, N, Tₘ, M, ε, qₙ, u, f1, f2) -> Vector

## Return
    Вектор размера `N-1`, сеточные значения правой части уравнения,
    для которого `g(x,t)` будет являться решением.

!!! warning
    Функция не предназначена для самостоятельного использования,
    она передается в качестве аргумента в `solve`, внутри которой
    для нее сформуются аргументы нужной длины, для которых не
    потребуется смещений индексов.
"""
function RP(y, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, u, f1, f2, w)
    arp = adjointRP(y, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, u, f1, f2, w)
    d = [ g_d(n, Xₙ, N, m, Tₘ, M, ε, qₙ, u, f1, f2, w) for n in 1:N-1]
    return arp .- d
end
## Якобиан сконструируем тем же способом, что и внутри пакета в `src/adjoint.jl`.
j(y, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, U, f1, f2, w) = ForwardDiff.jacobian( z -> adjointRP(z, m, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, U, f1, f2, w), y)

## С использованием автоматического дифференцирования
ψ = solve_adjoint(ψ₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, RP, j, w = ω)
@test all(isapprox.(ψ_model, ψ, atol = 0.025))

## С использованием трехдиагонального якобиана
ψ = solve_adjoint(ψ₀, Xₙ, N, Tₘ, M, ε, ψl, ψr, qₙ, Uₙₘ, f1_data, f2_data, RP, ∂ARP_∂y, w = ω)
@test all(isapprox.(ψ_model, ψ, atol = 0.025))

(ψ, ψ_model, Xₙ, Tₘ)
nothing #hide


text:
using Test
using CompScienceMeshes

rectangle = meshrectangle(1.0, 1.0, 1.0);

@test vertices(rectangle,1) == point(0.0, 0.0, 0.0)
@test vertices(rectangle,2) == point(0.0, 1.0, 0.0)
@test vertices(rectangle,3) == point(1.0, 0.0, 0.0)

@test numvertices(rectangle) == 4
@test numcells(rectangle) == 2
@test dimension(rectangle) == 2

faces = skeleton(rectangle, 2)
edges = skeleton(rectangle, 1)
verts = skeleton(rectangle, 0)

@test numcells(verts) == 4
@test numcells(edges) == 5
@test numcells(faces) == 2

Λ  = connectivity(verts, edges)
Σᵀ = connectivity(edges, faces)

@test size(Λ)  == (5,4)
@test size(Σᵀ) == (2,5)

@test norm(Σᵀ*Λ, Inf) ≈ 0

bnd = boundary(rectangle)

@test dimension(bnd) == 1
@test numvertices(bnd) == 4
@test numcells(bnd) == 4

vtoe, num_vtoe = vertextocellmap(edges)
vtof, num_vtof = vertextocellmap(faces)

@test size(vtoe) == (4,3)
@test size(vtof) == (4,2)

@test size(num_vtoe) == (4,)
@test size(num_vtof) == (4,)

file = joinpath(dirname(@__FILE__), "sphere.in")
sphere = readmesh(file)
@test numvertices(sphere) == 335
@test numcells(sphere) == 666

m = meshrectangle(1.0,1.0,1.0)
rotate!(m, [1,0,0]*(π/2))

tol = sqrt(eps())
@test norm(vertices(m,1) - point(0.0,0.0,0.0)) < tol
@test norm(vertices(m,2) - point(0.0,0.0,1.0)) < tol
@test norm(vertices(m,3) - point(1.0,0.0,0.0)) < tol
@test norm(vertices(m,4) - point(1.0,0.0,1.0)) < tol

## test fliporientation
V = [
    point(0,0,0),
    point(1,0,0),
    point(0,1,0),
]

F = [
    index(1,2,3)
]

m = Mesh(V,F)
n = fliporientation(m)

@test m.faces[1] == index(1,2,3)
@test n.faces[1] == index(2,1,3)


text:
export FGMRES

function FGMRES(A,
                b::AbstractVector{T},
                x0::AbstractVector{T};
                m::Integer=10,
                maxiter::Number=Inf,
                precond=nothing,
                tol::AbstractFloat=1e-6,
                outputfreq::Integer=0)::Tuple{Vector{T},Vector{Float64}} where {T<:Number}

    n = size(A,2)
    size(A,1)==size(A,2) || throw(ArgumentError("A must be square"))
    size(b,1)==n || throw(ArgumentError("A and b must have compatible dimensions, got size of A $(size(A,1)), size of b $(length(b))"))
    size(b,2)==1 || throw(ArgumentError("b must be a vector"))
    length(x0)==0 || length(x0)==n || throw(ArgumentErrror("x0 must be empty or a vector with size compatible with A"))
    if length(x0)==0 || norm(x0)<1e-10
        x = zeros(T,size(A,2))
        r = copy(b)
    else
        x = x0
        r = b - A*x
    end

    normr0 = norm(r)
    ej = zeros(T,m+1)
    ej[1] = 1
    if precond==nothing
        P = x->x
        prec_dirac = true
    elseif typeof(precond)<:Function
        P = precond
        prec_dirac = false
    elseif typeof(precond)<:joAbstractOperator
        P = x->precond*x
        prec_dirac = false
    end

    hst = Array{Float64,1}()
    push!(hst,1.0)
    res = Array{Float64,1}()
    push!(res,1.0)
    it_counter = 1
    V = Matrix{T}(undef, n,m+1)
    if prec_dirac
        Z = Matrix{T}(undef, n,0)
    else
        Z = Matrix{T}(undef, n,m+1)
    end
    H = zeros(T,m+1,m)
    y = zeros(T,m)
    w = zeros(T,n)

    while hst[end] > tol && div(it_counter,m)<maxiter

        beta = norm(r)
        @. V[:,1] = r/beta

        innerit = 1
        for k=1:m
            if prec_dirac
                w .= A*(V[:,k])
            else
                Z[:,k] .= (P(V[:,k]))::Vector{T}
                w .= A*(Z[:,k])
            end
            for j=1:k
                H[j,k] = dot((@view V[:,j]),w)
                w .-= H[j,k].*(@view V[:,j])
            end
            H[k+1,k] = norm(w)

            if k<m
                @. V[:,k+1] = w/H[k+1,k]
            end

            # Solve min_y || H[1:k+1,1:k]-beta*ej[1:k+1] ||_2
            # We won't bother with dot-assignment here because
            # the matrices are so small
            y[1:k] = H[1:k+1,1:k]\(beta*ej[1:k+1])

            push!(hst,norm(beta*ej[1:k+1]-H[1:k+1,1:k]*y[1:k])/normr0)
            innerit+=1

            if hst[end] < tol
                break
            end
            it_counter += 1
        end
        if prec_dirac
            for k=1:innerit-1
                x .+= y[k].*(@view V[:,k])
            end
        else
            for k=1:innerit-1
                x .+= y[k].*(@view Z[:,k])
            end
        end

        r .= b - A*x
        if outputfreq > 0 && mod(div(it_counter,m),outputfreq)==0
            @printf("it %3d res %3.3e\n",div(it_counter,m), norm(r)/normr0)
        end
        push!(res,norm(r))
        if norm(r)/normr0<tol
            break
        end
    end

    return (x,res)
end


text:
Presentation.Slide(
"""

# Presentation.jl

- Introduction

- Installation

- REPL

"""
)

Presentation.Slide(
"""

## Introduction

- Allows presenting Julia code from the terminal

- Allows defined predefined code to run

"""
)

Presentation.Slide(
"""

## Installation

- Run `Pkg.add("https://github.com/kdheepak/Presentation.jl")`

- Add `using Presentation` to `~/.juliarc.jl`

- Create slides in `slideshow.jl`

- Run `julia -i slideshow.jl`

"""
)

Presentation.Slide(
"""

## REPL

- Move forward slides by using Ctrl+f

- Move back slides by using Ctrl+b

- Input predefined expressions on a slide by using Ctrl+e

    - Slide(..., [
      :(1+1)
      :(println("hello world"))
      quote
          function add(x, y)
              x+y
          end
      end
      ]
    )

""", [
      :(1+1)
      :(println("hello world"))
      quote
          function add(x, y)
              x+y
          end
      end
     ]
)


text:
__precompile__(true)
module DataBench
using IndexedTables, PooledArrays, NamedTuples, RCall, JuliaDB
import Base.ht_keyindex
# using FastGroupBy

export R_bench, createIndexedTable, run_juliadb_bench, run_juliadb_bench_pmap
export gen_string_vec_fixed_len, gen_string_vec_var_len, gen_string_vec_id_fixed_len
export createSynDataFrame

include("julia/blogpost/indexedtable.jl_bench.jl")
include("julia/r_data_table_group_by/rbench.jl")
include("julia/r_data_table_group_by/DataFrames.jl")

end # module


text:
@testset "MISC" begin

	@test_throws ErrorException("Must select one convention (S options. Run gmthelp(velo) to learn about them)") velo!("",mat2ds([0. -8 0 0 4 6 0.5; -8 5 3 3 0 0 0.5], ["4x6", "3x3"]), region=(-15,10,-10,10))

	println("	MISC")
	G = nearneighbor(rand(100,3) * 150, R="0/150/0/150", I=1, N=4, S=10, r=true);
	G = GMT.mat2grid(G.z; reg=0, hdr=[G.range; G.registration; G.inc]);
	G1 = gmt("grdmath -R-2/2/-2/2 -I0.5 X Y MUL");
	G2 = G1;
	G3 = G1 + G2;
	G3 = G1 + 1;
	G3 = 1 + G1;
	G3 = G1 - G2;
	G3 = G1 - 1;
	G3 = G1 * G2;
	G3 = G1 * 2;
	G3 = 2 * G1;
	G3 = G1 ^ 2;
	G3 = -G1;
	G3 = G1 / G2;
	G3 = G1 / 2;
	G1 = mat2grid([0.0 1; 2 3]);
	G2 = mat2grid([4 5; 6 7; 8 9]);
	@test_throws ErrorException("The HDR array must have 9 elements") mat2grid(rand(4,4), reg=0, hdr=[0. 1 0 1 0 1]);
	@test_throws ErrorException("Grids have different sizes, so they cannot be added.") G1 + G2;
	@test_throws ErrorException("Grids have different sizes, so they cannot be subtracted.") G1 - G2;
	@test_throws ErrorException("Grids have different sizes, so they cannot be multiplied.") G1 * G2;
	@test_throws ErrorException("Grids have different sizes, so they cannot be divided.") G1 / G2;
	G1 = GMT.mat2grid(rand(4,4));
	G2 = GMT.mat2grid(rand(Float32,4,4), G1);
	GMT.mat2img(rand(UInt8,8,8), G1);
	G2 = mat2grid(rand(Float32,4,4), mat2img(rand(UInt16,32,32),x=[220800 453600], y=[3.5535e6 3.7902e6]));
	G2 = GMT.mat2grid(rand(Int32,4,4));
	G2 = GMT.mat2grid(rand(4,4));
	G1 .* G2;
	getindex(G1,1:2);
	setindex!(G1, [-1 -1],1:2)
	size(G1)

	GMT.WrapperPluto("aaa")
	Base.BroadcastStyle(typeof(G1))
	getindex(G1,1)

	GMT.find4similar(G1)
	GMT.find4similar(G1,0)
	GMT.find4similar(1)
	GMT.find4similar(())
	GMT.find4similar((0,1))
	GMT.find4similar([],0)
	I = mat2img(rand(UInt8,4,4,3))
	GMT.find4similar(I,0)
	size(I)
	Base.BroadcastStyle(typeof(I))
	getindex(I,1);
	setindex!(I, [101 1],1:2)
	I .+ UInt8(0)

	GMT.GMTdataset();
	isempty(GMT.GMTcpt())
	size(GMT.GMTcpt())
	isempty(GMT.GMTps())
	GMT.GMTdataset(rand(2,2), "lixo");
	GMT.GMTdataset(rand(Float32, 2,2), ["aiai"])
	GMT.GMTdataset(rand(Float32, 2,2), "aiai")
	GMT.GMTdataset(rand(Float32, 2,2))
	D = mat2ds(GMT.fakedata(4,4), x=:ny, color=:cycle, multi=true)
	D[1].text = ["lixo", "l", "p", "q"];
	GMT.find4similar(D[1],0)
	getindex(D[1],1);
	setindex!(D[1], 1,1)
	Base.BroadcastStyle(typeof(D[1]))
	display(D);
	plot(D, legend=true, Vd=dbg2)
	mat2ds(rand(5,4), x=:ny, color=:cycle, hdr=" -W1");
	mat2ds(rand(5,4), x=1:5, hdr=[" -W1" "a" "b" "c"], multi=true);
	@test_throws ErrorException("The header vector can only have length = 1 or same number of MAT Y columns") mat2ds(rand(2,3), hdr=["a" "b"]);

	GMT.mat2grid(rand(Float32, 10,10), reg=1);
	GMT.mat2grid(1, hdr=[0. 5 0 5 1 1])
	GMT.num2str(rand(2,3));
	text_record([-0.4 7.5; -0.4 3.0], ["a)", "b)"]);
	text_record(["aa", "bb"], "> 3 5 18p 5i j");
	text_record(["> 3 5 18p 5i j", "aa", "bb"]);
	text_record(Array[["aa", "bb"],["cc", "dd", "ee"]]);
	text_record([["aa", "bb"],["cc", "dd", "ee"]]);

	# TEST THE API DIRECTLY (basically to improve coverage under GMT6)
	PS = plot(rand(3,2), ps=1);
	API = GMT.GMT_Create_Session("GMT", 2, GMT.GMT_SESSION_NOEXIT + GMT.GMT_SESSION_EXTERNAL + GMT.GMT_SESSION_COLMAJOR);
	GMT.ps_init(API, PS, 0);
	gmt("destroy")

	# Test ogr2GMTdataset
	API = GMT.GMT_Create_Session("GMT", 2, GMT.GMT_SESSION_NOEXIT + GMT.GMT_SESSION_EXTERNAL + GMT.GMT_SESSION_COLMAJOR);
	gmtread("lixo.gmt");		# This "lixo.gmt" was created in test_avatars.jl
	if (GMTver > v"6.1.1")
		GMT.gmt_ogrread(API, "lixo.gmt", C_NULL);		# Keep testing the old gmt_ogrread API.
	end
	GMT.GMT_Get_Default(API, "API_VERSION", "        ");
	D = gmtconvert([1.0 2 3; 2 3 4], a="2=lolo+gPOINT");	# There's a bug in GMT for this. No data points are printed
	gmtwrite("lixo.gmt", D)
	@test gmtconvert([1.0 2 3; 2 3 4], binary_out="3f", write="a.bin", Vd=2) == "gmtconvert  > a.bin -bo3f"
	@test gmtconvert([1.0 2 3; 2 3 4], binary_out="3f", append="a.bin", Vd=2) == "gmtconvert  >> a.bin -bo3f"
	#GMT.ogr2GMTdataset(GMT.gmt_ogrread(API, "lixo.gmt"));
	rm("lixo.gmt")

	if (GMTver >= v"6.1")
		check = UInt8[zeros(9,9) ones(9,9) ones(9,9).*2; ones(9,9).*3 ones(9,9).*4 ones(9,9).*5; ones(9,9).*6 ones(9,9).*7 ones(9,9).*8];
		C = makecpt(range=(0,9,1));
		I = mat2img(check, cmap=C);
		rgb = GMT.ind2rgb(I);
		image_alpha!(I, alpha_ind=5);
		image_alpha!(I, alpha_vec=round.(UInt32,rand(6).*255));
		image_alpha!(I, alpha_band=round.(UInt8,rand(27,27).*255))
		image_cpt!(I, C)
		GMT.transpcmap!(I, true)
		GMT.transpcmap!(I, false)
		image_cpt!(I, clear=true)
	end

	GMT.linspace(1,1,100);
	GMT.logspace(1,5);
	GMT.fakedata(50,1);
	GMT.meshgrid(1:5, 1:5, 1:5);
	fields(7);
	fields(rand(2,2))
	tic();toc()
	@test_throws ErrorException("`toc()` without `tic()`") toc()

	# MB-System
	println("	MB-System")
	mbgetdata("aa", Vd=2)
	mbimport("aa", Vd=2)
	mbsvplist("aa", Vd=2)
	mblevitus(Vd=2)

	println("	EXAMPLES")
	# EXAMPLES
	plot(1:10,rand(10), lw=1, lc="blue", marker="square", markeredgecolor=:white, size=0.2, markerfacecolor="red", title="Hello World", xlabel="Spoons", ylabel="Forks")

	x = range(0, stop=2pi, length=180);	seno = sin.(x/0.2)*45;
	coast(region="g", proj="A300/30/6c", axis="g", resolution="c", land="navy")
	plot!(collect(x)*60, seno, lw=0.5, lc="red", marker="circle", markeredgecolor=0, size=0.05, markerfacecolor="cyan")

	G = GMT.peaks();
	show(G);
	show(mat2ds(rand(2,13), multi=true));
	grdcontour(G, cont=1, annot=2, axis="a")
	cpt = makecpt(T="-6/8/1");      # Create the color map
	grdcontour(G, axis="a", color=cpt, pen="+c", fmt=:png, savefig="lixo")
	D = grdcontour(G, cont=[-2,0,5], dump=true);

	GMT.geodetic2enu(-81.998,42.002,1000,-82,42,200);

	println("	DRAWING")
	circle(0,0,1,first=true,Vd=dbg2)
	cross(0,0,1,Vd=dbg2)
	custom(0,0, "bla", 1,Vd=dbg2)
	diamond(0,0,1,Vd=dbg2)
	hexagon(0,0,1,Vd=dbg2)
	itriangle(0,0,1,Vd=dbg2)
	letter(0,0, 1, "A", "Helvetica", "CM",Vd=dbg2)
	minus(0,0,1,Vd=dbg2)
	pentagon(0,0,1,Vd=dbg2)
	plus(0,0,1,Vd=dbg2)
	square(0,0,1,Vd=dbg2)
	star(0,0,1,Vd=dbg2)
	triangle(0,0,1,Vd=dbg2)
	ydash(0,0,1,Vd=dbg2)
	box(0,0,1,1,Vd=dbg2)
	rect(0,0,1,1,Vd=dbg2)
	ellipseAz(0,0, 0, 1, 1,Vd=dbg2)
	rotrect(0,0, 0, 1, 1,Vd=dbg2)
	rotrectAz(0,0, 0, 1, 1,Vd=dbg2)
	roundrect(0,0, 1, 1, 1,Vd=dbg2)
	ellipse(300,201,0, 200, 50, first=true, units=:points, fill=:purple, pen=1)
	circle(305,185,56, fill=:black, figname="lixo.ps")
end


text:
struct EffectiveChargeDensity{T, N, S} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S}
end

@inline size(ρ::EffectiveChargeDensity{T, N, S}) where {T, N, S} = size(ρ.data)
@inline length(ρ::EffectiveChargeDensity{T, N, S}) where {T, N, S} = length(ρ.data)
@inline getindex(ρ::EffectiveChargeDensity{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ρ.data, I...)
@inline getindex(ρ::EffectiveChargeDensity{T, N, S}, i::Int) where {T, N, S} = getindex(ρ.data, i)
@inline getindex(ρ::EffectiveChargeDensity{T, N, S}, s::Symbol) where {T, N, S} = getindex(ρ.grid, s)

function EffectiveChargeDensity(fss::PotentialSimulationSetup{T, N, S})::EffectiveChargeDensity{T, N, S} where {T, N, S}
    return EffectiveChargeDensity{T, N, S}( fss.ρ, fss.grid )
end

function NamedTuple(ρ::EffectiveChargeDensity{T}) where {T <: SSDFloat}
    return (
        grid = NamedTuple(ρ.grid),
        values = ρ.data * internal_voltage_unit,
    )
end
Base.convert(T::Type{NamedTuple}, x::EffectiveChargeDensity) = T(x)

function EffectiveChargeDensity(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    EffectiveChargeDensity{T, N, S}( ustrip.(uconvert.(internal_voltage_unit, nt.values)), grid)
end
Base.convert(T::Type{EffectiveChargeDensity}, x::NamedTuple) = T(x)


text:
# The functions dogravdlnMh() and dothermdlnMh() are needed only for "PlotFigure2.jl"
using MatterPower
using HaloMF
using DoubleExponentialFormulas

function dogravdlnMh(
   pk,
   z::Real,
   Ωm::Real,
   lnMh::Real,
   Ωcb = Ωm;
   virial = false,
   t10MF = false,
)
   ρc = 2.775e11 # in units of h^2 M⊙/Mpc^3
   GN = 3 / 8π / 2998^2 / ρc # = 4.786e-20 Mpc/M⊙ is Newton's constant
   if virial # Determine Δm used by the mass function
      E2 = Ωm * (1 + z)^3 + 1 - Ωm # flat Universe
      Ωmz = Ωm * (1 + z)^3 / E2
      Δvir = 18 * π^2 + 82 * (Ωmz - 1) - 39 * (Ωmz - 1)^2
      Δm = Δvir / Ωmz
   else
      Δm = 200
   end
   spl = dndlnMh(pk, z, Δm, Ωcb, lnMh, t10MF = t10MF)
   # %% Compute ρgrav = (1/2) * \int dlnM dn/dlnM Ag*GM^2/R, in units of h^2 M⊙/Mpc^3
   # c(M) relation is due to Equation (2) and "NFW, F, 0-2" of Table 1 in
   # Duffy et al., MNRAS, 390, L64 (2008)
   if virial # For Δ = Δvir(z)
      RΔh = cbrt(exp(lnMh) * 3 / 4π / (ρc * E2 * Δvir)) # in h^-1 Mpc
      A0, B0, C0 = 7.85, -0.081, -0.71 # Avir, Bvir, Cvir in Table 1 of Duffy et al.
   else # For Δ = Δm = 200
      RΔh = cbrt(exp(lnMh) * 3 / 4π / (ρc * Ωcb * Δm)) / (1 + z) # in h^-1 Mpc
      A0, B0, C0 = 10.14, -0.081, -1.01 # Amean, Bmean, Cmean in Table 1 of Duffy et al.
   end
   c = A0 * (exp(lnMh) / 2e12)^B0 * (1 + z)^C0
   dρdlnMh = 0.5 * spl * GN * exp(2 * lnMh) / RΔh * Ag(c)
   halfW = -dρdlnMh / ρc / Ωcb
   Ωgrav = halfW * Ωcb
end

function dothermdlnMh(
   pk,
   z::Real,
   Ωm::Real,
   h0::Real,
   lnMh::Real,
   Ωcb = Ωm;
   massbias = 1.266,
   t10MF = false,
   xmax = 6.0,
)
   ρc = 2.775e11 # in units of h^2 M⊙/Mpc^3
   ρceVcm3 = 1.05375e4 # in units of h^2 eV/cm^3
   Δc = 500 # 500 times the critical density of the Universe
   E2 = Ωm * (1 + z)^3 + 1 - Ωm  # flat Universe
   Δm = Δc * E2 / Ωm / (1 + z)^3 # Δm for mass function
   spl = dndlnMh(pk, z, Δm, Ωcb, lnMh, t10MF = t10MF)
   # %% Pressure profile integral
   # Planck's universal pressure profile
   # Reference: Planck Collaboration, A&A, 550, A131 (2013)
   c, γ, α, β, P0 = [1.81, 0.31, 1.33, 4.13, 6.41]
   upp(x) =
      x^2 * P0 * (0.7 / h0)^1.5 * (c * x)^(-γ) * (1 + (c * x)^α)^((γ - β) / α)
   uppint, err = quadde(upp, 0, xmax)
   # %% Compute ρth = \int dlnM dn/dlnM \int dV Pe, in units of h^2 eV/cm^3
   αp = 0.12
   E2 = Ωm * (1 + z)^3 + 1 - Ωm # for a flat Universe
   RΔh = cbrt(exp(lnMh) * 3 / 4π / (ρc * E2 * Δc)) # in h^-1 Mpc
   Pe =
      1.65 * # eV/cm^3
      (1 / 0.7)^2 * # h^2 cancels h^2 in the critical density in the denominator
      E2^(4 / 3) *
      (exp(lnMh) / massbias / 3e14 / 0.7)^(2 / 3 + αp)
   dρdlnMh = spl * Pe * 4π * RΔh^3 / massbias
   Y = 0.24 # helium abundance
   Ωtherm = (8 - 5Y) / (4 - 2Y) * dρdlnMh * uppint / ρceVcm3
end

"""
    dndlnMh(pk, z, Δm, Ωm, lnMh; t10MF = false)

This function returns a value of the comoving halo number density per log mass interval in units of h^3/Mpc^3 at a specified value of the log(halo mass).

# Arguments
- `pk`(k): a function which returns a linear matter power spectrum with the argument k being the comoving wavenumber.
   - `pk` is in units of (Mpc/h)^3.
- `z::Real`: redshift.
- Δm: halo overdensity with respect to the mean mass density of the Universe
- `Ωm::Real`: present-day matter density parameter. When excluding the neutrino contribution, use Ωcb (baryon + cold dark matter density parameter).
- `lnMh::Real`: log(halo mass)
   - halo mass is in units of M⊙/h.

# Optional keyword arguments
- `t10MF::Bool=false`: if `true`, use `tinker10MF` for the halo multiplicity function. If `false` (the default), use `tinker08MF`.
"""
function dndlnMh(pk, z, Δm, Ωm, lnMh; t10MF = false)
   ρc = 2.775e11 # in units of h^2 M⊙/Mpc^3
   Rh = cbrt(exp(lnMh) * 3 / 4π / ρc / Ωm) # in units of Mpc/h
   σ2 = sigma2(pk, Rh)
   dlnσ2dlnRh = Rh * dsigma2dR(pk, Rh) / σ2
   lnν = 2 * log(1.6865) - log(σ2)
   if t10MF
      MF = tinker10MF(lnν, z, Δm)
   else
      MF = tinker08MF(lnν, z, Δm)
   end
   dndlnMh = -dlnσ2dlnRh * MF / 4π / Rh^3 # in units of h^3 Mpc^-3
end

"""
   Ag(c:Real; method = "poly4")

Integral of the squared Fourier transform of an NFW density profile, ``(1/π)∫_0^∞ dx |u(x,c)|^2``.

*Reference*: Equation (TBD) of Chiang, Makiya, Komatsu & Ménard (in prep)

# Arguments
- `c::Real`: concentration parameter.

# Optional keyword arguments
- `method::String="poly4"`: approximation method.
   - `poly4` (the default) or other values for 4th order polynomial fit.
   - `poly3` for 3rd order polynomial fit.
   - `exact` for numerical integration.
"""
function Ag(c::Real; method = "poly4")
   if method == "poly3" # 3rd order polynomial fit
      y = c - 5
      Ag = 1.0184 + 5.389e-2 * y - 1.514e-3 * y^2 + 3.845e-5 * y^3
   elseif method == "exact" # exact integration
      f = log(1 + c) - c / (1 + c)
      u2(x) =
         (
            sin(x / c) * (sinint((1 / c + 1) * x) - sinint(x / c)) +
            cos(x / c) * (cosint((1 / c + 1) * x) - cosint(x / c)) -
            sin(x) / (1 / c + 1) / x
         )^2
      res, err = quadde(u2, 0, 100)
      Ag = res / π / f^2
   else # 4th order polynomial fit
      y = c - 5
      Ag =
         1.0202 + 5.376e-2 * y - 1.842e-3 * y^2 + 9.489e-5 * y^3 -
         2.352e-6 * y^4
   end
end


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_eccentricities(obj::Antenna, arg0::Frequency)
    return jcall(obj, "getEccentricities", Vector3D, (Frequency,), arg0)
end

function get_frequencies(obj::Antenna)
    return jcall(obj, "getFrequencies", List, ())
end

function get_phase_center_variation(obj::Antenna, arg0::Frequency, arg1::Vector3D)
    return jcall(obj, "getPhaseCenterVariation", jdouble, (Frequency, Vector3D), arg0, arg1)
end

function get_sinex_code(obj::Antenna)
    return jcall(obj, "getSinexCode", JString, ())
end

function get_type(obj::Antenna)
    return jcall(obj, "getType", JString, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end


text:
"""
SARSA-λ learner (core)
"""

mutable struct SarsaL
    nr_states::Int
    nr_actions::Int
    α::Float64
    α_eff::Float64
    γ::Float64
    λ::Float64
    initvalue::Float64
    QSA::Array{Float64, 2}
    eligTrace::Array{Float64, 2}
    rpe::Array{Float64, 1}
    use_replacing_trace::Bool
    do_debug_print::Bool
end
function SarsaL(; nr_states = 10, nr_actions = 4, α = 0.1, γ = .9,
    λ = 0.9, initvalue = 0.,
    QSA = zeros(Float64, nr_states, nr_actions) .+ initvalue,
    use_replacing_trace = true,
    do_debug_print = false)

    SarsaL(nr_states, nr_actions, α, α, γ, λ, initvalue, QSA,
        zeros(Float64, nr_states, nr_actions), [0.],
        use_replacing_trace, do_debug_print)
end
export SarsaL
function do_debug_print_final(learner::Union{SarsaL, SarsaLModulated})
    println("Final:")
    for s in 1:size(learner.QSA,1)
        @show learner.QSA[s, :]
    end
    println("--------------")
end

""" Update: main function """
function update!(learner::Union{SarsaL, SarsaLModulated}, state, action, next_state, next_action, reward,
                isendofexperiment, modulatingfactor)

    if learner.use_replacing_trace # replacing trace
        learner.eligTrace[state,action] = 1.
    else # (accumulating) eligibility trace
        learner.eligTrace[state,action] += 1.
    end

    getlearningrate!(learner, modulatingfactor)

    if reward>0 #goal state, reward, NO NEXT ACTION (in our paradigm)
        learner.rpe[1] = (reward - learner.QSA[state, action])
        learner.QSA[:, :] .+= learner.α_eff * learner.eligTrace * learner.rpe[1] #e_trace for the current state is 1
        learner.eligTrace = zeros(Float64, learner.nr_states, learner.nr_actions) # reset ET at end of episode
    else #non goal state
        if !isendofexperiment
            learner.rpe[1] = reward + learner.γ * learner.QSA[next_state, next_action] - learner.QSA[state, action]
            learner.QSA[:, :] .+= learner.α_eff * learner.eligTrace * learner.rpe[1]
            learner.eligTrace[:, :] .*= learner.λ * learner.γ
        else
            learner.rpe[1] = 0. # end of experiment
        end
    end
    # @show state, action, next_state, reward
    # @show learner.rpe[1]
    # @show learner.eligTrace[state, 1], learner.eligTrace[state, 2]
    # @show learner.QSA[state, 1], learner.QSA[state, 2]
    # @show learner.QSA
    if learner.do_debug_print
        do_debug_print_final(learner)
    end
end
export update!
function getlearningrate!(learner::Union{SarsaL, REINFORCE, REINFORCEbaseline}, modulatingfactor)
    learner.α_eff = learner.α
end


text:
function get_operator(process::AbstractBlackScholesProcess, grid::Vector{Float64}, residualTime::Float64, ::Bool)
  # if timeDependent
  #   return BSMTermOperator(grid, process, residualTime)
  # else
  return BSMOperator(grid, process, residualTime)
  # end
end

function BSMOperator(grid::Vector{Float64}, process::AbstractBlackScholesProcess, residualTime::Float64)
  tridiag = TridiagonalOperator(length(grid))
  logGrid = LogGrid(grid)
  cc = PdeConstantCoeff(process, residualTime, state_variable(process).value)
  generate_operator!(cc, residualTime, logGrid, tridiag)

  return tridiag
end


text:
for fn_fits in fns
    local file = FITS(fn_fits)
    local l1 = read(file[2],"wavelength")
    local f1 = read(file[2],"spectrum")
    local v1 = read(file[2],"uncertainty")
    fn_jld2 = replace(fn_fits,"fits"=>"jld2")
    println(fn_jld2)
    @save fn_jld2 λ=l1 flux=f1 var=v1
end

@time begin
    for fn_fits in fns
    fn_jld2 = replace(fn_fits,"fits"=>"jld2")
    local f = jldopen(fn_jld2, "r")
    local l1 = f["λ"]
    local f1 = f["flux"]
    local v1 = f["var"]
    println("sum(l) = ", sum(view(l1[:,10:20],.!isnan.(l1[:,10:20]))),
    " sum(f) = ", sum(view(f1[:,10:20],.!isnan.(f1[:,10:20]))),
    " sum(v) = ", sum(view(v1[:,10:20],.!isnan.(v1[:,10:20]))) )
end
end

@time begin for fn_fits in fns
    local file = FITS(fn_fits)
    local l1 = read(file[2],"wavelength")
    local f1 = read(file[2],"spectrum")
    local v1 = read(file[2],"uncertainty")
    println("sum(l) = ", sum(view(l1[:,10:20],.!isnan.(l1[:,10:20]))),
    " sum(f) = ", sum(view(f1[:,10:20],.!isnan.(f1[:,10:20]))),
    " sum(v) = ", sum(view(v1[:,10:20],.!isnan.(v1[:,10:20]))) )
end
end

@time begin
    for fn_fits in fns
    fn_jld2 = replace(fn_fits,"fits"=>"jld2")
    local f = h5open(fn_jld2, "r")
    local l1 = readmmap(f["λ"])
    local f1 = readmmap(f["flux"])
    local v1 = readmmap(f["var"])
    println("sum(l) = ", sum(view(l1[:,10:20],.!isnan.(l1[:,10:20]))),
    " sum(f) = ", sum(view(f1[:,10:20],.!isnan.(f1[:,10:20]))),
    " sum(v) = ", sum(view(v1[:,10:20],.!isnan.(v1[:,10:20]))) )
end
end


text:
module q3pkg

ExtractMthEnd(daily) = floor(daily./100)
end
 # module