text:
using Parameters
# This file sets up the data-loading infrastructure.
export @gid_str, load_glacier

## Note some types are also defined in data-structs.jl

"""
$(TYPEDEF)

Types of different data to be dealt with
"""
abstract type DataKind end
# glacier outline
abstract type OutlineData <: DataKind end
# misc model parameters
abstract type ParaData <: DataKind end
# DEM 2D
abstract type DEMData <: DataKind end
# ice surface velocity
abstract type IVData <: DataKind end
# Surface mass balance
abstract type BdotData <: DataKind end
# dh/dt data
abstract type DhdtData <: DataKind end
# basal sliding fraction
abstract type FslData <: DataKind end
# ice temperature
abstract type TempData <: DataKind end
# ice thickness (e.g. radar or nothing)
abstract type ThicknessData <: DataKind end
# ice flux at terminus
abstract type TerminusFluxData <: DataKind end

"Standard variable names for each DataKind"
datakind2varname = Dict(OutlineData=>:outline,
                        #ParaData=>:para,
                        DEMData=>:dem,
                        IVData=>:iv,
                        BdotData=>:bdot,
                        DhdtData=>:dhdt,
                        FslData=>:fsl,
                        TempData=>:temp,
                        ThicknessData=>:thick)

"""
Returns a dict of default options to individual loaders.  Add opts
here if they are universally used.  For others, modify the dict
returned from this function.
"""
function def_dataset_opts()
    out = Dict{DataType,Any}()
    out[DataKind] = Dict{Symbol,Any}()
    for st in subtypes(DataKind)
        out[st] = Dict{Symbol,Any}()
    end
    out[DEMData][:dist_land_maxdist] = 3e3
    # Set windows to volume-area scaling thickness estimate within below bounds:
    out[ParaData][:pm] = Dict{Symbol,Any}(
        :window_dem_smooth => (0,200),
        :window =>  (0,2000),
        )
    out[ParaData][:pn] = Dict{Symbol,Any}()
    out[ParaData][:pp] = Dict{Symbol,Any}()
    return out
end

function def_dataset_files()
    out = Dict{DataType,Vector{String}}()
    for st in subtypes(DataKind)
        out[st] = [""]
    end
    return out
end
function def_dataset_dateranges()
    out = Dict{DataType,Tuple{Date,Date}}()
    for st in subtypes(DataKind)
        out[st] = (Date(),Date())
    end
    return out
end

# Loading parameters, used to customize the loading process
"""
$(TYPEDEF)

Misc parameters used during loading/synthetic-creation of datasets used for
**customization** of the loading process
Any parameters which are used inside `load_glacier` (with the
exception of some physical parameters) should be in here.

Particular sets of glaciers should overload the constructor of this to
return a suitable pl. Example

    LoadPara(gid::ITMIXGlacier, use_glogem=false; kwargs...) = ...

The flow goes: `LoadPara` instance is used to create the DataSet
instances, which in terms make up the `dt` DataTable instance.  The
`dt` then holds all info needed for data-loading/data-creation:
- Input/output directories and file-paths of data-files.
- The `dataset_opts` present here will be inserted into the `DataSet` instance.
  The options in `dataset_opts[DataKind]` will be added to all via a merge.

    code_root::String=normpath(dirname(dirname(@__FILE__))) # root directory of this package
    data_root::String=normpath(joinpath(code_root, "..", "data")) # root directory of data
    cache_dir::String=joinpath(data_root, "cache") # directory to write cache-files
    update_cache::Bool=false # if true, recalculate all cached files.
    target_proj::String="from DEM" # in what projection to re-project data-sets.  Default "from DEM" uses the projection of the DEM.
    # The following fields mirror the DataSet fields:
    dataset_opts::Dict{DataType,Any}=def_dataset_opts() # holds options for load! functions.
    dataset_LOADERS::Dict{DataType,Symbol}=Dict{DataType,Symbol}() # Dict between the two type parameters for the DataSet instances.
    dataset_files::Dict{DataType,Vector{String}}=def_dataset_files()
    dataset_dateranges::Dict{DataType,Tuple{Date,Date}}=def_dataset_dateranges()
"""
@with_kw struct LoadPara<:APara @deftype String
    gid::GlacierID
    code_root::String=normpath(dirname(@__DIR__)) # root directory of this package
    data_root::String=normpath(joinpath(code_root, "scripts", "data")) # root directory of data
    cache_dir::String=joinpath(data_root, "cache") # directory to write cache-files
    update_cache::Bool=false # if true, recalculate all cached files.
    target_proj::String="from DEM" # in what projection to re-project data-sets.  Default "from DEM" uses the projection of the DEM.
    # The following fields mirror the DataSet fields:
    dataset_opts::Dict{DataType,Any}=def_dataset_opts() # holds options for load! functions.
    dataset_LOADERS::Dict{DataType,Symbol}=Dict{DataType,Symbol}() # Dict between the two type parameters for the DataSet instances.
    dataset_files::Dict{DataType,Vector{String}}=def_dataset_files()
    dataset_dateranges::Dict{DataType,Tuple{Date,Date}}=def_dataset_dateranges()
end

"""
$(TYPEDEF)
$(FIELDS)

Holds information needed to load or create a data-set, which can be a
single data-set or a set of related one (e.g. a DEM and its mask.
Although there the outline comes in too, which is separate...).

Parameterized on the DataKind it holds and the LOADER
it needs to be loaded.

Common opts:
- proj -- expected projection
"""
@with_kw struct DataSet{DK<:DataKind, LOADER}
    "File names, if applicable."
    files::Vector{String}
    "Extra options for loader functions.  Unknown options are (usually) ignored."
    opts::Dict{Symbol,Any}
    "Date range of measurements.  (Date(),Date()) if unknown (yet)."
    daterange::Tuple{Date,Date}
    "File-cache directory"
    datacache_dir::String # File-cache directory.  Note that there may
                          # be several cache-files associated with a
                          # DataSet. Say several masks for a DEM containing
                          # several glaciers.
    "If true, recalculate the cache files"
    datacache_dir_update::Bool
    "If true, do not use the RAM cache"
    datacache_ram_update::Bool
    function DataSet{DK,L}(files, opts, daterange=(Date(),Date()),
                           datacache_dir="", datacache_dir_update=false,
                           datacache_ram_update=false) where {DK,L}
        new{DK,L}(files, opts, daterange, datacache_dir, datacache_dir_update,
                  datacache_ram_update)
    end
end
Base.convert(::Type{DataSet{DK,L}}, ds::DataSet) where {DK,L} =
    DataSet{DK,L}(ds.files, ds.opts, ds.daterange, ds.datacache_dir, ds.datacache_dir_update, ds.datacache_ram_update)
Base.show{DK,LS}(io::IO, ds::DataSet{DK,LS}) =
    println(io, "DataSet{$DK,$LS}: $(ds.files), $(ds.daterange)")

# for caching to RAM use hash:
function Base.hash(ds::DataSet{DK,L}, h::UInt) where {DK,L}
    h += hash(DK,h)
    h += hash(L,h)
    h += hash(ds.files,h)
    h += hash(ds.opts,h)
    h += hash(ds.daterange,h)
    h += hash(ds.datacache_dir,h)
    h += hash(ds.datacache_dir_update,h)
    return h
end

function Base.:(==){DK1,LS1,DK2,LS2}(ds1::DataSet{DK1,LS1}, ds2::DataSet{DK2,LS2})
    DK1==DK2 || return false
    LS1==LS2 || return false
    for fl in fieldnames(ds1)
        getfield(ds1,fl)==getfield(ds2,fl) || return false
    end
    return true
end

"""
A table which holds where/how to load/produce the data for a
particular glacier and a particular run.  (i.e. a glacier may have
several datasets)

Fields:

    gid::GlacierID
    outline
    para
    dem
    iv
    bdot
    dhdt
    fsl
    temp
    thickness

"""
@with_kw_noshow struct DataTable
    gid::GlacierID
    outline::DataSet{OutlineData}
    para::DataSet{ParaData}
    dem::DataSet{DEMData}
    iv::DataSet{IVData}
    bdot::DataSet{BdotData}
    dhdt::DataSet{DhdtData}
    fsl::DataSet{FslData}
    temp::DataSet{TempData}
    thickness::DataSet{ThicknessData}
    terminus_flux::DataSet{TerminusFluxData}
end
function Base.show(io::IO, dt::DataTable)
    println(io, "Data Table for Glacier: $(dt.gid)")
    for fl in fieldnames(dt)
        show(io, getfield(dt, fl))
    end
    println("")
end
"Return all known datakinds in the order used in DataTable"
datakinds() = [_get_first_para(ST)for ST in DataTable.types[2:end]]
_get_first_para(::Type{DataSet{T}}) where {T} = T

"""
$(SIGNATURES)
$(METHODLIST)

Constructs a `DataTable` for a glacier with `GlacierID` and `pl`
loading options.  It will choose reasonable defaults for a glacier for
all the data inputs.  The output, a `DataTable` can later be updated to
choose different data sources for some fields.
"""
function make_datatable(gid::GlacierID, pl::LoadPara, dt=nothing)
    @assert gid==pl.gid "Not matching gid: $gid, $(pl.gid)"
    @unpack dataset_LOADERS, dataset_files,
            dataset_opts, dataset_dateranges = pl
    args = []
    for DK in datakinds()
        push!(args, DataSet{DK,dataset_LOADERS[DK]}(
            dataset_files[DK],
            merge(dataset_opts[DataKind],dataset_opts[DK]),
            dataset_dateranges[DK],
            pl.cache_dir,
            pl.update_cache,
            pl.update_cache
        ))
    end
    return DataTable(gid, args...), pl
end

# """
# Hold data for a model run.
# """
# struct ModelRunData
#     datatable::DataTable
#     outfile::String # path to a JLD file
#     comments::Dict
# end
# function ModelRunData(gid, outfile, comments=Dict())
#     datatable = getdatatable(gid)
#     ModelRunData(datatable, outfile, comments)
# end

#########
# Loader functions
#########
# of form

"""
    load!(ds::DataSet{OutlineData, LS}, gid::GlacierID, target_proj)
    load!(ds::DataSet{DEMData, LS}, gid::GlacierID, target_proj, outline, pp::Phys)
    load!(ds::DataSet{OtherData, LS}, gid::GlacierID, target_proj, outline, cropbox,
          pp::Phys)

Data loader/synthesizing functions.  Overload `load!` functions for
needed combinations of `DataKind` and `LOADER` (some symbol, often
only one combo makes sense).  The `ds.opts` hold options for specific
loaders.

The `load!` function modifies `ds` in that it stores the raw(ish)
loaded data into `ds.datacache`.  This allows caching data like DEMs
which maybe applicable to several glaciers.  If the `ds.datacache` is
already filled then the reading from disc step is skipped.

The primary variable it returns should be the loaded data tailored to
the glacier at hand.
"""
function load! end

"""
    load!(::DataSet{OutlineData, :dummy}, gid::GlacierID)

Loading the outline determines which glacier to run.

Return
- outline: the outline of the glacier including land-islands.
- outcrops: nunataks, or similar outcrops which are not included in
  outline.  Their insides are all treated as land for masking
  purposes.  (This is a bit redundant with the outline, but some data
  comes in this format.)  Return (
- iscalving::Bool -- true for calving glacier (tide-water, fresh-water, hanging glacier)
- proj -- projection as Proj4 string (or something else)
"""
function load!(::DataSet{OutlineData, :dummy}, gid::GlacierID)
    if target_proj=="from DEM"
        proj = use_Outline_proj() # if needed, it get re-projected later
    end

    return outline::Outline, iscalving::Bool, proj::String
end

"""
    load!(::DataSet{ParaData}, gid::GlacierID, outline::Outline, isicecap::Bool)

Load parameters for forward model.  Probably this generic
implementation is fine.  This default loader returns the default
para-structs which are modifed by `ds.opts[:pp]`, etc.  Can
additionally be overloaded for `gid`.

The `outline` and `isicecap` input are used to automatically set the
smoothing windows to approximate mean ice thickness (estimated with
volume area scaling).

Return
- pp::Phys
- pm::MPara
- pn::Num
TODO: add inverse parameters too.
"""
function load!(ds::DataSet{ParaData}, gid::GlacierID, outline::Outline, isicecap::Bool)
    opts = ds.opts
    if isa(opts[:pm][:window],Tuple) || isa(opts[:pm][:window_dem_smooth],Tuple)
        # set bed-DEM smoothing window with volume-area
        hmean = volume_area_mean_h(area(outline), isicecap=isicecap)
        if isa(opts[:pm][:window],Tuple)
            mi,ma = opts[:pm][:window]
            opts[:pm][:window] = min(ma, max(mi,1.5*hmean) )
        end
        if isa(opts[:pm][:window_dem_smooth],Tuple)
            mi,ma = opts[:pm][:window_dem_smooth]
            opts[:pm][:window_dem_smooth] = min(ma, max(mi,hmean) )
        end
    end
    pp = Phys(; opts[:pp]...)
    pm = MPara(; opts[:pm]...)
    pn = Num(; opts[:pn]...)
    return pp::Phys, pm::MPara, pn::Num
end

"""
    load!(ds::DataSet{DEMData, :dummy}, outline,
               pp::Phys)

ds.opts:
- grid : if specified use a different grid from the grid of the DEM as
  computational grid.  Given in the `target_proj`.
- cropbox -- used for cropping the DEM and other 2D data.  Needs to
  contain all of the outline. `F[xmin, xmax, ymin, ymax]`, defaults
  being calculated during DEM-load with some space around the outline
  (needed for interpolations).

Return:
- dem::Gridded
- glaciermask::Matrix{Bool} -- true on glacier
- landmask::Matrix{Bool} -- true on solid ground (i.e. false on glaciers, ice shelves, water)
- dist_land::Matrix{F} -- distance to closest land
- dist_land_maxdist -- maximum distance to which to calculate above
- cropbox::Vector{Float64} -- box around region of interest (contains glacier + some buffer)
                              `F[xmin, xmax, ymin, ymax]`
"""
function load!(::DataSet{DEMData, :dummy}, outline,
               pp::Phys)
    return dem::Gridded, glaciermask::Matrix{Bool}, landmask::Matrix{Bool},
           dist_land::Matrix{F}, dist_land_maxdist::F, cropbox::Vector{Float64},
           dem_proj::String
end
function load!(::DataSet{IVData, :dummy}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    return iv::SomeData, ivmask::Matrix{Bool}
end
function load!(::DataSet{ThicknessData, :dummy}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    return thickness::SomeData
end
function load!(::DataSet{BdotData, :dummy}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    return bdot::SomeData
end
function load!(::DataSet{DhdtData, :dummy}, dem, glaciermask, target_proj, outline, cropbox,
              iscalving, pp::Phys)
    return dhdt::SomeData
end
function load!(::DataSet{FslData, :dummy}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    return fsl::SomeData
end
function load!(::DataSet{TempData, :dummy}, dem, glaciermask, target_proj, outline, cropbox,
              pp::Phys)
    return temp::SomeData
end
# scalar loaders
"""
    load!(::DataSet{TerminusFluxData}, iscalving, pp::Phys)

Ice flux at terminus.
"""
function load!(::DataSet{TerminusFluxData}, iscalving, pp::Phys)
    error()
end

"""
     special_processing!(gid::GlacierID, outline, dem, glaciermask, landmask,
                             iv, ivmask, thickness, bdot, dhdt, fsl, temp,
                             pl::LoadPara, pp::Phys, pm::MPara, pn::Num)

Allows to do some special processing on a glacier by glacier basis by
overloading this function.
"""
function special_processing!(gid::GlacierID, outline, dem, glaciermask, landmask,
                             iv, ivmask, thickness, bdot, dhdt, fsl, temp,
                             pl::LoadPara, pp::Phys, pm::MPara, pn::Num)
    nothing
end

"""
$(SIGNATURES)

Load/synthesize all data for a model run of a glacier using either:
- a `DataTable`
- a `DataTable` and an already "loaded" `DataTable`
- a glacier ID (TODO: implement)
as input.  For the latter the default dataset is used.  (It
potentially modifies the `_datacache_ram` field of the `DataSet`s, but
this is an internal fields, thus no !."

Input
- glacier-id or datatable
- load-para

Return
- Glacier
- pp, pm, pn
"""
function load_glacier end

function load_glacier(gid::GlacierID, pl::LoadPara)
    dt, pl = make_datatable(gid, pl)
    load_glacier(dt,pl)
end

function load_glacier(datatable::DataTable, pl::LoadPara)
    @assert datatable.gid==pl.gid  "Not matching gid: $gid, $(pl.gid)"
    @unpack gid, outline, para, dem, iv, bdot, dhdt,
            fsl, temp, thickness, terminus_flux = datatable

    @unpack target_proj = pl # TODO do cache dir

    # TODO: maybe reverse loading order: first load DEM with grid, then
    # outline.  With outline come the cropbox & masks, which would be more
    # consistent.  Make DEM & IV masks their separate loader?  dist_land
    # is also generated with the masks.

    outline, iscalving, isicecap, outline_proj = load!(outline, gid)

    # parameters, updated using the optional input parameters
    pp, pm, pn = load!(para, gid, outline, isicecap)

    # note some variables get re-bound:
    (dem, glaciermask, landmask, dist_land,
     dist_land_maxdist, cropbox, dem_proj) = load!(dem, outline, pp)

    @assert !any(dem.v[glaciermask].==FILL) "The gl.dem has FILL values inside the glacier mask, probably mask those points"
    # projections
    target_proj = target_proj=="from DEM" ? dem_proj : target_proj
    @assert dem_proj==outline_proj==target_proj "Need to have same projection"
    # dem, glaciermask, landmask, dist_land, cropbox, target_proj =
    #     reproject_DEM(dem, glaciermask, landmask, dist_land, cropbox, dem_proj, target_proj)
    # outline, outcrops = reproject_outline(outline, outcrops, outline_proj, target_proj)

    ## TODO update about projections
    iv, ivmask = load!(iv, dem, glaciermask, target_proj, outline, cropbox, pp)
    thickness = load!(thickness, dem, glaciermask, target_proj, outline, cropbox, pp)
    terminus_flux = load!(terminus_flux, dem, glaciermask, iscalving, pp)
    bdot = load!(bdot, dem, glaciermask, target_proj, outline, cropbox, pp)
    dhdt = load!(dhdt, bdot, terminus_flux, dem, glaciermask, target_proj,
                 outline, cropbox, iscalving, pp)
    fsl = load!(fsl, dem, glaciermask, target_proj, outline, cropbox, pp)
    temp = load!(temp, dem, glaciermask, target_proj, outline, cropbox, pp)

    special_processing!(gid, outline, dem, glaciermask, landmask,
                        iv, ivmask, thickness, bdot, dhdt, fsl, temp,
                        pl, pp, pm, pn)
    alpha, elemean, elemax, elemin, elemedian = proc_dem(dem, glaciermask)

    return Glacier(gid,
                   iscalving, isicecap,
                   target_proj,
                   dem,
                   outline, cropbox,
                   glaciermask, landmask,
                   dist_land, dist_land_maxdist,
                   alpha, elemean, elemax, elemin, elemedian,
                   bdot, dhdt,
                   temp, fsl,
                   terminus_flux,
                   thickness, iv, ivmask,
                   Dict{Symbol,Any}()),
           pp,pm,pn,pl
end

# function load_glacier(datatable::DataTable, datatable_other::DataTable, pl::LoadPara)
#     @unpack gid, outline, para, dem, iv, bdot, dhdt, fsl, temp, thickness = datatable
#     ot = datatable_other
#     for fl in fieldnames(datatable)
#         fl==:gid && continue
#         if compare_nocache(getfield(datatable, fl), getfield(datatable_other, fl))
#             setfield!(datatable, fl) = getfield(datatable_other, fl)
#         end
#     end
#     load_glacier(datatable, pl)
# end

################
# All together
################

"""
    init_forward(gid::GlacierID, verbose=true; kwargs_for_LoadPara...)

Initialize forward model.  Returns

-  gl,gb,pp,pm,pn,pl
"""
function init_forward(gid::GlacierID, verbose=true; kwargs_for_LoadPara...)
    verbose && print("Loading glacier $gid ... ")
    pl = LoadPara(gid; kwargs_for_LoadPara...)
    t = @elapsed gl,pp,pm,pn,pl = load_glacier(gid,pl)
    if verbose
        println("$(round(t,3)) sec")
        print("Making bands ... ")
    end
    pn.test && Base.Test.@inferred make_bands(gl, pp, pm, pn)
    t = @elapsed gb,pm = make_bands(gl, pp, pm, pn)
    verbose && println("$(round(t,3)) sec")
    return gl,gb,pp,pm,pn,pl
end

#################
# Helper functions see data-loading-misc.jl
################


text:
# install script for Linux & MacOS machines.

using Pkg

function install()
    checkIneFEMpart()
    intro()
    skipline()
    checkOK()
    skipline()
    checkJuliaVersion()
    skipline()
    addStartupFile()
    skipline()
    install_fgt()
    skipline()
    addModules()
    skipline()
    finished()
end

function checkIneFEMpart()
    if !(basename(pwd()) == "eFEMpart")
        error("current directory is not eFEMpart")
    end
end

function intro()
    println("You have chosen to install eFEMpart,")
    println("a package for implementing finite element")
    println("code with particle simulators in Julia")
end

"""
    input(prompt::String="")::String

Read a string from STDIN. The trailing newline is stripped.

The prompt string, if given, is printed to standard output without a
trailing newline before reading input.
"""
function input(prompt::String)
    print(prompt)
    return chomp(readline())
end

function checkOK()
    println("This install script does the following:")
    println("   - Adds eFEMpart packageto your Julia LOAD_PATH")
    println("   - Downloads and installs figtree package")
    println("   - Adds the following dependencies:")
    println("      + JLD")
    println("      + Plots")
    println("      + PyPlot")
    println("      + PyCall")
    println("      + LaTeXStrings")
    println("      + Parameters")

    keepasking = true
    while keepasking
        OK = input("Is this OK? (y/n) ")

        if OK == "y"
            keepasking = false
        elseif OK == "n"
            error("User did not want to install")
        else
            println("Sorry, that isn't an option.")
        end
    end

    println()
    println("Great! Installation will begin now.")
    println("This could take several minutes...")
end

function skipline()
    println()
end

function checkJuliaVersion()
    print("Checking correct Julia version...")
    if VERSION>v"1.0.0"
        println("OK")
    else
        println("FAILED")
        error("Julia version ",VERSION," needs to be at least 1.0")
    end
end

function addStartupFile()
    print("Adding eFEMpart to your Julia PATH...")
    eFEMpartPATH = string(pwd(),"/src/")
    newLoadPath = string("push!(LOAD_PATH,\"",eFEMpartPATH,"\")")

    # create startup.jl file
    mkpath("$(homedir())/.julia/config")
    touch("$(homedir())/.julia/config/startup.jl")

    fr = open("$(homedir())/.julia/config/startup.jl","r")

    # check that LOAD_PATH doesn't already contain what you're about to add
    alreadyThere = false
    lines = readlines(fr)
    for i=1:length(lines)
        if occursin(newLoadPath, lines[i])
            alreadyThere = true
        end
    end
    close(fr)

    # add startup file
    if !(alreadyThere)
        fa = open("$(homedir())/.julia/config/startup.jl","a")
        write(fa,"\n",newLoadPath)
        close(fa)
    end

    # load startup file -- so that you don't have to restart julia
    include("$(homedir())/.julia/config/startup.jl")

    println("OK")
end

function addModules()
    # JLD
    println("Adding JLD package...")
    Pkg.add("JLD");

    # IterativeSolvers
    println("Adding IterativeSolvers package...")
    Pkg.add("IterativeSolvers");

    # Preconditioners
    println("Adding Preconditioners package...")
    Pkg.add("Preconditioners");

    # BenchmarkTools
    println("Adding BenchmarkTools package...")
    Pkg.add("BenchmarkTools");

    # Plots
    println("Adding Plots package...")
    Pkg.add("Plots");

    # PyPlot
    println("Adding PyPlot package...")
    Pkg.add("PyPlot");

    # PyCall
    println("Adding PyCall package...")
    Pkg.add("PyCall");

    # LaTeXStrings
    println("Adding LaTeXStrings package...")
    Pkg.add("LaTeXStrings");

    # Parameters
    println("Adding Parameters package...")
    Pkg.add("Parameters");
end

function install_fgt()
    println("installing figtree...")
    #run(`sudo mkdir -p src/part`)
    mkpath("src/part")
    cd("src/part")
    # download zip
    downloadurl = "https://sourceforge.net/projects/figtree/files/latest/download"
    run(`curl -OL $(downloadurl)`);

    # unzip
    run(`unzip download`);
    #run(`unzip download && mv figtree download`)

    # compile
    cd("figtree-0.9.3")
    run(`make`);

    # move relevant *.so file
    run(`mv lib/libfigtree.so ../libfigtree.so`)
    run(`mv lib/libann_figtree_version.so ../libann_figtree_version.so`)

    cd("..")

    # remove unnecessary package
    run(`rm -r figtree-0.9.3/ download`)

    # add *.so folder to library path
    pwdir = pwd()
    bashstart  = "$(homedir())/.bashrc"
    tcshstart = "$(homedir())/.tcshrc"

    # add LD_LIBRARY_PATH to bash startup
    f=open(bashstart,"a")
    println(f)
    println(f,"export LD_LIBRARY_PATH=$(pwdir):\$LD_LIBRARY_PATH")
    close(f)

    # re-run bash file to load LD_LIBRARY_PATH
    cmd = `bash -c "source $(bashstart)"`
    run(cmd)

    # add LD_LIBRARY_PATH to tcsh startup
    #f = open(tcshstart,"a")
    #println(f)
    #println(f,"export LD_LIBRARY_PATH=$(pwdir)")
    #close(f)

    # re-run tcsh file to load LD_LIBRARY_PATH
    #cmd = `tcsh -c "source $(tcshstart)"`
    #run(cmd)

    cd("../..")

    println("succesfully installed figtree!")
end
function finished()
    skipline()
    println("All components installed correctly!")
    println("You can now use eFEMpart!")
end

install()

ext: jl

max_stars_repo_name: pseastham/eFEMpart

max_issues_repo_name: pseastham/eFEMpart


text:
# Benchmark on reduction

# macros

macro bench_reduc1(ptname, name, rp, f, a)
    quote
        println("    on $($name) ...")

        ($f)($a)
        t0 = @elapsed for i in 1 : ($rp)
            ($f)($a)
        end

        ($f)($a, 1)
        t1 = @elapsed for i in 1 : ($rp)
            ($f)($a, 1)
        end

        ($f)($a, 2)
        t2 = @elapsed for i in 1 : ($rp)
            ($f)($a, 2)
        end

        push!($ptname, ($name, [t0, t1, t2] * 1000.))
    end
end

macro bench_reduc2(ptname, name, rp, f, a, b)
    quote
        println("    on $($name) ...")

        ($f)($a, $b)
        t0 = @elapsed for i in 1 : ($rp)
            ($f)($a, $b)
        end

        ($f)($a, $b, 1)
        t1 = @elapsed for i in 1 : ($rp)
            ($f)($a, $b, 1)
        end

        ($f)($a, $b, 2)
        t2 = @elapsed for i in 1 : ($rp)
            ($f)($a, $b, 2)
        end

        push!($ptname, ($name, [t0, t1, t2] * 1000.))
    end
end

# data preparation

a2 = randn(1000, 1000)
b2 = rand(1000, 1000) .+ 0.5
u = randn(1000)

#################################################
#
#  Benchmark on old functions
#
#################################################

typealias DIMS Union(Int, (Int, Int))

_sum(x::Array) = sum(x)
_sum(x::Array, d::DIMS) = sum(x, d)

_maximum(x::Array) = maximum(x)
_maximum(x::Array, d::DIMS) = maximum(x, d)

_minimum(x::Array) = minimum(x)
_minimum(x::Array, d::DIMS) = minimum(x, d)

_sumabs(x::Array) = sum(abs(x))
_sumabs(x::Array, d::DIMS) = sum(abs(x), d)

_maxabs(x::Array) = maximum(abs(x))
_maxabs(x::Array, d::DIMS) = maximum(abs(x), d)

_minabs(x::Array) = minimum(abs(x))
_minabs(x::Array, d::DIMS) = minimum(abs(x), d)

_sumsq(x::Array) = sum(abs2(x))
_sumsq(x::Array, d::DIMS) = sum(abs2(x), d)

_dot(x::Array, y::Array) = sum(x .* y)
_dot(x::Array, y::Array, d::DIMS) = sum(x .* y, d)

_sumabsdiff(x::Array, y::Array) = sum(abs(x - y))
_sumabsdiff(x::Array, y::Array, d::DIMS) = sum(abs(x - y), d)

_maxabsdiff(x::Array, y::Array) = maximum(abs(x - y))
_maxabsdiff(x::Array, y::Array, d::DIMS) = maximum(abs(x - y), d)

_minabsdiff(x::Array, y::Array) = minimum(abs(x - y))
_minabsdiff(x::Array, y::Array, d::DIMS) = minimum(abs(x - y), d)

_sumsqdiff(x::Array, y::Array) = sum(abs2(x - y))
_sumsqdiff(x::Array, y::Array, d::DIMS) = sum(abs2(x - y), d)

_entropy(x::Array) = -sum(x .* log(x))
_entropy(x::Array, d::DIMS) = -sum(x .* log(x), d)

_sumxlogy(x::Array, y::Array) = sum(x .* log(y))
_sumxlogy(x::Array, y::Array, d::DIMS) = sum(x .* log(y), d)

_var(x::Array) = var(x)
_var(x::Array, d::DIMS) = var(x, d)

_std(x::Array) = std(x)
_std(x::Array, d::DIMS) = std(x, d)

_logsumexp(x::Array) = (u = maximum(x); log(sum(exp(x .- u))) .+ u)
_logsumexp(x::Array, d::DIMS) = (u = maximum(x, d); log(sum(exp(x .- u))) .+ u)

function _softmax(x::Array)
    u = maximum(x)
    r = exp(x .- u)
    r / sum(r)
end

function _softmax(x::Array, d::Int)
    u = maximum(x, d)
    r = exp(x .- u)
    r ./ sum(r, d)
end

# benchmark

const oldperf = Array((ASCIIString, Vector{Float64}), 0)

println("Benchmark results on Base methods:")

@bench_reduc1 oldperf "sum" 10 _sum a2
@bench_reduc1 oldperf "maximum" 10 _maximum a2
@bench_reduc1 oldperf "minimum" 10 _minimum a2

@bench_reduc1 oldperf "sumabs" 10 _sumabs a2
@bench_reduc1 oldperf "maxabs" 10 _maxabs a2
@bench_reduc1 oldperf "minabs" 10 _minabs a2
@bench_reduc1 oldperf "sumsq" 10 _sumsq a2

@bench_reduc2 oldperf "dot" 10 _dot a2 b2
@bench_reduc2 oldperf "sumabsdiff" 10 _sumabsdiff a2 b2
@bench_reduc2 oldperf "maxabsdiff" 10 _maxabsdiff a2 b2
@bench_reduc2 oldperf "minabsdiff" 10 _minabsdiff a2 b2
@bench_reduc2 oldperf "sumsqdiff" 10 _sumsqdiff a2 b2

@bench_reduc1 oldperf "entropy" 10 _entropy b2
@bench_reduc2 oldperf "sumxlogy" 10 _sumxlogy b2 b2
@bench_reduc1 oldperf "var" 10 _var a2
@bench_reduc1 oldperf "std" 10 _std a2
@bench_reduc1 oldperf "logsumexp" 10 _logsumexp a2
@bench_reduc1 oldperf "softmax" 10 _softmax a2

push!(oldperf, ("varm", [NaN, NaN, NaN]))

#################################################
#
#  Benchmark on new functions
#
#################################################

using NumericExtensions
const newperf = Array((ASCIIString, Vector{Float64}), 0)

println("Benchmark results in New methods:")

@bench_reduc1 newperf "sum" 10 sum a2
@bench_reduc1 newperf "maximum" 10 maximum a2
@bench_reduc1 newperf "minimum" 10 minimum a2

@bench_reduc1 newperf "sumabs" 10 sumabs a2
@bench_reduc1 newperf "maxabs" 10 maxabs a2
@bench_reduc1 newperf "minabs" 10 minabs a2
@bench_reduc1 newperf "sumsq" 10 sumsq a2

@bench_reduc2 newperf "dot" 10 dot a2 b2
@bench_reduc2 newperf "sumabsdiff" 10 sumabsdiff a2 b2
@bench_reduc2 newperf "maxabsdiff" 10 maxabsdiff a2 b2
@bench_reduc2 newperf "minabsdiff" 10 minabsdiff a2 b2
@bench_reduc2 newperf "sumsqdiff" 10 sumsqdiff a2 b2

@bench_reduc1 newperf "entropy" 10 entropy b2
@bench_reduc2 newperf "sumxlogy" 10 sumxlogy b2 b2
@bench_reduc1 newperf "var" 10 var a2
@bench_reduc1 newperf "std" 10 std a2
@bench_reduc1 newperf "logsumexp" 10 logsumexp a2
@bench_reduc1 newperf "softmax" 10 softmax a2

println("    on varm ...")
varm(a2, 1.0)
t0 = @elapsed for i = 1:10; varm(a2, 1.0); end
varm(a2, u, 1)
t1 = @elapsed for i = 1:10; varm(a2, u, 1); end
varm(a2, u, 2)
t2 = @elapsed for i = 1:10; varm(a2, u, 2); end
push!(newperf, ("varm", [t0, t1, t2] * 1000.))

#################################################
#
#  Organize and print results
#
#################################################

function organize_results(name, raw)
    tab = BenchmarkTable(name, ["full-reduction", "colwise-reduction", "rowwise-reduction"])
    for e in raw
        rname = e[1]
        row = e[2]
        add_row!(tab, rname, row)
    end
    tab
end

function compute_gains(name, tab0::BenchmarkTable, tab1::BenchmarkTable)
    # compare colnames and rownames
    @assert tab0.colnames == tab1.colnames
    @assert tab0.rownames == tab1.rownames

    gtab = BenchmarkTable(name, tab0.colnames)
    m = nrows(tab0)
    for i in 1 : m
        add_row!(gtab, tab0.rownames[i], tab0.rows[i] ./ tab1.rows[i])
    end
    gtab
end

oldtable = organize_results("Matrix reduction using Base methods (millisec)", oldperf)
newtable = organize_results("Matrix reduction using New methods (millisec)", newperf)

gaintable = compute_gains("Speed gain of New vs. Base", oldtable, newtable)

println("\nResults:")
println("***********************\n")

println(oldtable)
println()

println(newtable)
println()

println(gaintable)
println()


text:
using Test
using Gravityfield

let
    @test V(1.0m) == (4π / 3)m³
    @test ρ(m = (4π / 3)kg, r = 1m) == 1kg/m³
    @test ρ(r = 1m, m = (4π / 3)kg) == 1kg/m³
    @test ρ_earth |> ustrip |> typeof == Float64
end

@testset "Gravity acceleration" begin
    @test g(1, 0) == 0.0m^3*kg^-1*s^-2
    @test g(1., 1.) == γ
    @test g(m_earth, r_earth) ≈ 10.0391664m/s^2
    ω_earth_sun = 2π / yr
    a_earth_sun = ω_earth_sun^2 * d_sun_earth
    @test g(m_sun, d_sun_earth) / a_earth_sun ≈ 1.000 atol = 0.001
end

@testset "Force vector on 1 due to 2" begin
    p1 = [0.,0.,0.]m
    p2 = [0.,1.,0.]m
    m1 = 316227kg
    m2 = 129099kg
    f1_2 = m1* uvec(p1, p2) * g(m2, len(p2-p1))
    @test f1_2 |> length == 3
    @test ustrip.(f1_2) |> typeof == Vector{Float64}
    @test f1_2 ≈ [0.0, 2.7247474, 0.]N
end
@testset "Euclidean vector length" begin
    @test len(-1m) == 1m
    @test len([1]m) == 1m
    @test len([2,0]m) == 2m
    @test len([1m,1m]) == sqrt(2)m
end
@testset "Vector between bodies or points" begin
    p1 = [0.,0.,0.]m
    p2 = [0.,1.,0.]m
    s = Sphere(p1, m_earth, r_earth)
    @test rvec(s,s) == [0., 0., 0.]m
    @test len(s,s) == 0m
    @test uvec(s,s) == [0., 0., 0.]
    @test rvec(p2,s) == [0., -1.0, 0.]m
    @test len(p2,s) == 1m
    @test uvec(p2,s) == [0., -1.0, 0.]

    s2 = Sphere(p2, m_earth, r_earth)
    @test rvec(s,s2) == [0., 1.0, 0.]m
    @test len(s,s2) == 1m
    @test uvec(s,s2) == [0., 1.0, 0.]

    s3 = Sphere([-d_sun_earth], 1, 1)
    s4 = Sphere([0.0m], 1, 1)
    @test rvec(s3,s4) == [d_sun_earth]
    @test len(s4,s3) == d_sun_earth
    @test uvec(s3,s4) == [1]
    @test uvec(s4,s3) == [-1]
end

@testset "Gravity acceleration vector btw. bodies or point and body" begin
    p1 = [0.,0.,0.]m
    s1 = Sphere(p = p1, m = 1kg, r = 1m)
    @test gvec(s1,s1) == [0., 0., 0.]m/s²
    @test gvec([-1, 0., 0.]m, s1) == [γ*kg/m², 0m/s², 0m/s²]
    s2 = Sphere([0.0]m, m_earth, r_earth)
    @test gvec([r_earth], s2) == [-10.039166439909296]m/s²
    @test gvec([r_earth/2], s2) == [-10.039166439909296]m/s²/2
end

@testset "Volume" begin
    @test V(1m) == (4π/3)m³
    @test V(r_sun) * ρ_sun / m_sun ≈ 0.999 atol = 0.001
    @test V(r_sun) isa Quantity
end

@testset "Sphere constructors" begin
    s1 = Sphere(p = [0]m, r = r_earth)
    s2 = Sphere([0]m, m_earth, r_earth)
    @test s1.m == s2.m
    @test s1.r == s2.r
    s3 = Sphere()
    @test s3.m ≈ s2.m
    @test s3.r ≈ s2.r
    s4 = earthlikeSphere_r_known(r = 10.E6m)
    s5 = earthlikeSphere_m_known(m = ρ_earth * V(10.E6m))
    @test s4.m ≈ s5.m
    @test s4.r ≈ s5.r
end

@testset "Position vectors" begin
    @test [1, 2, 3]m isa Pos
    @test !isa([1 2; 3 4]m, Pos)
    @test !isa([1m, 2, 3], Pos)
    @test !isa([1, sin, 3], Pos)
    @test !isa([[1, 2, 3, 4,5,6]], Pos)
    @test [1, 2, 3, 4,5,6]m isa Pos
    @test !isa([[1, 2, 3, 4,5,6]], Pos)
    @test !isa([1, 2, 3, 4,5,6]kg, Pos)
end

@testset "Gravity acceleration vector in system of bodies" begin
    earth = Sphere()
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0m, 0m], m = m_moon)
    sun = Sphere([-d_sun_earth, 0m, 0m], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    bodies = [sun, earth, moon]
    @test gvec(earth.p, bodies)[1] < 0m/s² # Sun stronger than moon in centre of earth
    @test gvec([r_earth/2, 0m, 0m], bodies)[1] < -5m/s² # Halfway to centre of earth, earth dominates
    xgrid, ygrid = xygrids( xmin = -d_moon_earth/2.,
        xmax = d_moon_earth * 3/2,
        ymin =-d_moon_earth/2,
        ymax = d_moon_earth * 3/2, n = 100)
    gmagn = gfieldvals(nest(xgrid, ygrid, moon.r),  bodies)
    @test foldr(+, gmagn)  ≈  222.5649149449058m/s²
    @test gvec(earth, bodies) == gvec(earth.p, bodies)
end

@testset "Nested position vectors" begin
    earth = Sphere()
    moon = earthlikeSphere_m_known(p = [d_moon_earth, 0m, 0m], m = m_moon)
    sun = Sphere([-d_sun_earth, 0m, 0m], m_sun, (3m_sun / (4π * ρ_sun))^(1/3))
    bodies = [sun, earth, moon]
    p1 = Pos([1., 2.]m)
    p2 = Pos([-1., -2.]m)
    vp = Npos([p1, p2])
    @test vp == [[1., 2.], [-1., -2.]]m
    mp = hcat(vp, vp)
    @test mp isa Npos
    pu1 = [1., 2.]
    pu2 = [-1., -2.]
    mup = hcat([pu1, pu2], [pu1, pu2])
    @test mp == mup*m
    @test nest(1m:1m:3m, 2m, 3m) == [[1, 2, 3], [2, 2, 3], [3, 2, 3]]m
    @test nest(1m,2m:1m:3m, 4m) == [[1, 2, 4], [1, 3, 4]]m
    @test nest(1m,2m, 3m:1m:4m) == [[1, 2, 3], [1, 2, 4]]m
    @test nest(1m:1m:2m, 1m:1m:2m) == [[1,1], [2,2]]m
    @test nest(1m,2m:1m:3m) == [[1,2], [1,3]]m
    @test nest(1m:1m:2m, 3m) == [[1,3], [2,3]]m
    @test nest(range(1.5m, 3.5m, length = 2)) == [[1.5], [3.5]]m
    @test gfieldvals(nest(1m,2m:1m:3m,4m), bodies)[1] ≈ 0.0059000661748m/s²
    @test nest([1m, 2m, 3m], [4, 5, 6]m, [7, 8, 9]m) ==  [[1, 4, 7], [2, 5, 8], [3, 6, 9]]m
    np = nest([1, 10]m, [2, 20]m, [3, 30]m, [4, 40]m, [5, 50]m, [6, 60]m)
    @test np ==  [[1, 2, 3, 4, 5, 6]m, [10, 20, 30, 40, 50, 60]m]
    @test slice_1(np) == [1, 10]m
    @test slice_2(np) == [2, 20]m
    @test slice_3(np) == [3, 30]m
    @test slice_4(np) == [4, 40]m
    @test slice_5(np) == [5, 50]m
    @test slice_6(np) == [6, 60]m
    @test gvec([[0m,0m,0m]], bodies)[1][1] ≈ -0.005898468352763m/s²
end

@testset "Bodies constructors" begin
    @test eltype(bodies_3d()) == Sphere
    @test eltype(bodies_2d()) == Sphere
    @test eltype(bodies_1d()) == Sphere
end

@testset "Pretty-printing" begin
    function sout(p)
        iob=IOBuffer()
        show(IOContext(iob, :color=>true), p)
        take!(iob)|>String
    end
    p = [1,2]m
    @test sout(p) == "[1.0, 2.0]\e[36mm\e[39m"
    function slongout(p)
        iob=IOBuffer()
        show(IOContext(iob, :color=>true), :"text/plain", p)
        take!(iob)|>String
    end
    @test slongout(p) == "Pos{Float64}([1.0, 2.0]\e[36mm\e[39m)"
    vp = [p, p]
    @test sout(vp) == "[[1.0, 2.0], [1.0, 2.0]]\e[36mm\e[39m"
    @test slongout(vp) == "2-element Npos{Pos{Float64}}:\n[[1.0, 2.0], [1.0, 2.0]]\e[36mm\e[39m"
    p1 = [1,1]m
    p2 = [2,2]m
    p3 = [3,3]m
    p4 = [4,4]m
    mp = hcat([p1, p2], [p3, p4])
    @test mp[2,1] == p2
    @test sout(mp) == "[[1.0, 1.0] [3.0, 3.0]; [2.0, 2.0] [4.0, 4.0]]\e[36mm\e[39m"
    @test slongout(mp) == "2x2 Npos{Pos{Float64}}:\n[[1.0, 1.0] [3.0, 3.0]; [2.0, 2.0] [4.0, 4.0]]\e[36mm\e[39m"
    @test sout(bodies_3d()) == "Sphere[Sphere(p = [-1.496e11, 0.0, 0.0]\e[36mm\e[39m, m = 1.989e30 kg, r = 6.95733e8 m), Sphere(p = [0.0, 0.0, 0.0]\e[36mm\e[39m, m = 5.97e24 kg, r = 6.3e6 m), Sphere(p = [3.844e8, 0.0, 0.0]\e[36mm\e[39m, m = 7.34767e22 kg, r = 1.45456e6 m)]"
    @test sout(bodies_2d()) == "Sphere[Sphere(p = [-1.496e11, 0.0]\e[36mm\e[39m, m = 1.989e30 kg, r = 6.95733e8 m), Sphere(p = [0.0, 0.0]\e[36mm\e[39m, m = 5.97e24 kg, r = 6.3e6 m), Sphere(p = [3.844e8, 0.0]\e[36mm\e[39m, m = 7.34767e22 kg, r = 1.45456e6 m)]"
    @test sout(bodies_1d()) == "Sphere[Sphere(p = [-1.496e11]\e[36mm\e[39m, m = 1.989e30 kg, r = 6.95733e8 m), Sphere(p = [0.0]\e[36mm\e[39m, m = 5.97e24 kg, r = 6.3e6 m), Sphere(p = [3.844e8]\e[36mm\e[39m, m = 7.34767e22 kg, r = 1.45456e6 m)]"
    @test slongout(bodies_3d()) == "3-element Array{Sphere,1}:\n  p = [-1.496e11, 0.0, 0.0]\e[36mm\e[39m m = 1.989e30 kg r = 6.95733e8 m\n  p = [0.0, 0.0, 0.0]\e[36mm\e[39m m = 5.97e24 kg r = 6.3e6 m\n  p = [3.844e8, 0.0, 0.0]\e[36mm\e[39m m = 7.34767e22 kg r = 1.45456e6 m"
    @test slongout(bodies_2d()) == "3-element Array{Sphere,1}:\n  p = [-1.496e11, 0.0]\e[36mm\e[39m m = 1.989e30 kg r = 6.95733e8 m\n  p = [0.0, 0.0]\e[36mm\e[39m m = 5.97e24 kg r = 6.3e6 m\n  p = [3.844e8, 0.0]\e[36mm\e[39m m = 7.34767e22 kg r = 1.45456e6 m"
    @test slongout(bodies_1d()) == "3-element Array{Sphere,1}:\n  p = [-1.496e11]\e[36mm\e[39m m = 1.989e30 kg r = 6.95733e8 m\n  p = [0.0]\e[36mm\e[39m m = 5.97e24 kg r = 6.3e6 m\n  p = [3.844e8]\e[36mm\e[39m m = 7.34767e22 kg r = 1.45456e6 m"
end

nothing

max_stars_repo_name: hustf/Gravityfield

max_issues_repo_name: hustf/Gravityfield

max_forks_repo_name: hustf/Gravityfield


text:
mutable struct BedgraphHeader{T} #TODO: determine what and how this will be.
    data::T
end

function Base.convert(::Type{String}, header::BedgraphHeader{Vector{String}}) :: String

    str = ""
    for line in header.data
        str = string(str, line, '\n')
    end

    return str
end

function generateBasicHeader(records::Vector{Record}; bump_forward=true) :: BedgraphHeader{Vector{String}} #Note: we assume that records are sorted by chrom and left position.

    chrom = records[1].chrom

    pos_start = records[1].first
    pos_end = records[end].last

    if bump_forward
        pos_start = pos_start + 1
        pos_end = pos_end + 1
    end

    return BedgraphHeader(["browser position $chrom:$pos_start-$pos_end", "track type=bedGraph"])
end

generateBasicHeader(chrom::String, pos_start::Int, pos_end::Int; bump_forward=true) = generateBasicHeader([Record(chrom, pos_start, pos_end, 0)], bump_forward=bump_forward)

function _readHeader(io) :: Vector{String}
    position(io) == 0 || seekstart(io)

    header = String[]
    line = readline(io)

    while !eof(io) && !isLikeRecord(line) # TODO: seek more rebust check.
        push!(header, line)
        line = readline(io)
    end

    return header
end

function Base.read(io::IO, ::Type{BedgraphHeader{Vector{String}}}) :: BedgraphHeader{Vector{String} }
    return BedgraphHeader(_readHeader(io))
end

function Base.write(io::IO, header::BedgraphHeader{Vector{String}})
    return Base.write(io, convert(String, header))
end


text:
struct QR{TQ,TR}
    Q::TQ
    R::TR
end
struct QnoR{TQ}
    Q::TQ
end
struct RnoQ{TR}
    R::TR
end
struct QRBuffers{TB,TH}
    β::TB
    h::TH
end
QRBuffers(A) = QRBuffers(zeros(eltype(A), min(size(A)...)), zeros(eltype(A), size(A,1)))
struct QRCompact{T,S}
    A::T
    β::S
end
QRCompact(A) = QRCompact(A, zeros(eltype(A), min(size(A)...)))

function qrfact!(A, qrb::QRBuffers=QRBuffers(A))
    qrfact!(A, qrb.β, qrb.h)
    QRCompact(A, qrb.β)
end
function qrfact!(qrc::QRCompact, h=zeros(eltype(qrc.A), size(qrc.A,1)))
    qrfact!(qrc.A, qrc.β, h)
    qrc
end

function QR(A, ::Type{Void})
    m, n = size(A)
    qr = QR(similar(A), zeros(eltype(A), min(m,n), n))
    qr
end
function QR!(qr::QR, qrc::QRCompact)
    accumulate!(qr.Q, qrc.A, qrc.β)
    getR!(qr.R, qrc.A)
    qr
end
function QR(qrc::QRCompact)
    m, n = size(qrc.A)
    qr = QR(similar(qrc.A), zeros(eltype(qrc.A), min(m,n), n))
    QR!(qr, qrc)
    qr
end

function QnoR!(qnor::QnoR, qrc::QRCompact)
    accumulate!(qnor.Q, qrc.A, qrc.β)
    qnor
end
function QnoR(qrc::QRCompact)
    qnor = QnoR(similar(qrc.A))
    QnoR!(qnor, qrc)
    qnor
end

function RnoQ!(rnoq::RnoQ, qrc::QRCompact)
    getR!(rnoq.R, qrc.A)
    rnoq
end
function RnoQ(qrc::QRCompact)
    m, n = size(qrc.A)
    rnoq = RnoQ(zeros(eltype(qrc.A), min(m,n), n))
    RnoQ!(rnoq, qrc)
    rnoq
end

const qr! = qrfact!


text:
# This file was generated, do not modify it. # hide
Xnew = DataFrame([hpn.^i for i in 1:10])
yy5 = predict(mtm, Xnew)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy5, lw=3)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

savefig("assets/literate/ISL-lab-5-g4.svg") # hide


text:
"""
        AnimalBehavior
"""
module AnimalBehavior

using Turing, StructArrays, Distributions, Random, PrettyTables
using StatsFuns: softmax!, logsumexp
using ForwardDiff: ForwardDiff
using DataFrames: DataFrames

import Base: rand, convert, show
import Turing: sample, loglikelihood, dic

# Models
include("models/evolutions.jl")
include("models/observations.jl")
include("models/macros.jl")

# Inference
include("inference/check_types.jl")
include("inference/sampling.jl")
include("inference/Posterior.jl")
include("inference/criteria.jl")
include("inference/posterior_sampling.jl")
include("inference/summarystats.jl")

# Simulation
include("simulation/simulate.jl")
include("simulation/interface.jl")
include("simulation/Simulation.jl")

export  @evolution,
        @observation,
        @model, # from Turing
        sample, # from Turing
        posterior,
        sample_hyperparams,
        sample_latent,
        expectation,
        dic,
        waic,
        bic,
        simulate,
        delta_rule!,
        epsilon_argmax,
        epsilon_greedy!,
        softmax!, # from StatsFuns
        ucb!
end

max_stars_repo_name: sqwayer/AnimalBehavior

max_issues_repo_name: sqwayer/AnimalBehavior


text:
using ConfigEnv
using PolygonIO
using TypedTables
using Test

dotenv()

const API_KEY = ENV["API_KEY"]
const tabular_opts = PolyOpts(API_KEY, Table)
const regular_opts = PolyOpts(API_KEY, nothing)

@testset "Reference API" begin
    # tickets test
    @test tickers(tabular_opts, "bitcoin") |> length == 10
    @test tickers(regular_opts, "bitcoin") |> length == 10

    # ticker_types test
    @test ticker_types(tabular_opts) |> length == 2
    @test ticker_types(regular_opts) |> length == 2

    # ticker_details test
    @test ticker_details(tabular_opts, "AAPL") |> length == 1
    @test ticker_details(regular_opts, "AAPL") |> length >= 28

    # ticker_details_vX next version test
    @test ticker_details_vX(tabular_opts, "AAPL", "2019-07-31") |> length == 1
    @test ticker_details_vX(regular_opts, "AAPL", "2019-07-31") |> length == 19

    # ticker_news test
    @test ticker_news(tabular_opts, "AAPL") |> length == 10
    @test ticker_news(regular_opts, "AAPL") |> length == 10

    # markets test
    @test markets(regular_opts) |> length >= 7
    @test markets(tabular_opts) |> length >= 7

    # locales test
    @test locales(regular_opts) |> length >= 19
    @test locales(tabular_opts) |> length >= 19

    # stock_splits test
    @test stock_splits(regular_opts, "AAPL") |> length >= 4
    @test stock_splits(tabular_opts, "AAPL") |> length >= 4

    # stock_dividends test
    @test stock_dividends(regular_opts, "AAPL") |> length >= 65
    @test stock_dividends(tabular_opts, "AAPL") |> length >= 65

    # stock_financials test
    @test stock_financials(regular_opts, "AAPL") |> length == 5
    @test stock_financials(tabular_opts, "AAPL") |> length == 5

    # market_holidays test
    @test market_holidays(regular_opts) |> length >= 1
    @test market_holidays(tabular_opts) |> length >= 1

    # market_status test
    @test market_status(regular_opts) |> length >= 6
    @test market_status(tabular_opts) |> length >= 1

    # stock_exchanges test
    @test stock_exchanges(regular_opts) |> length >= 30
    @test stock_exchanges(tabular_opts) |> length >= 30

    # condition_mappings test
    @test condition_mappings(regular_opts, "trades") |> length >= 1
    @test condition_mappings(tabular_opts, "quotes") |> length >= 1

    # stock_exchanges test
    @test crypto_exchanges(regular_opts) |> length >= 5
    @test crypto_exchanges(tabular_opts) |> length >= 5
end

@testset "Stock API" begin
    # trades test
    @test stocks_trades(tabular_opts, "AAPL", "2020-10-14") |> length == 10
    @test stocks_trades(regular_opts, "AAPL", "2020-10-14") |> length == 10

    # quotes_nbbo test
    @test stocks_quotes_nbbo(tabular_opts, "AAPL", "2020-10-14") |> length == 10
    @test stocks_quotes_nbbo(regular_opts, "AAPL", "2020-10-14") |> length == 10

    # last_trade_symbol test
    @test stocks_last_trade_symbol(tabular_opts, "AAPL") |> length == 1
    @test stocks_last_trade_symbol(regular_opts, "AAPL") |> length >= 1

    # last_quote_symbol test
    @test stocks_last_quote_symbol(regular_opts, "AAPL") |> length >= 1
    @test stocks_last_quote_symbol(tabular_opts, "AAPL") |> length == 1

    # daily_open_close test
    @test stocks_daily_open_close(tabular_opts, "AAPL", "2020-10-14") |> length == 1
    @test stocks_daily_open_close(regular_opts, "AAPL", "2020-10-14") |> length >= 10

    # grouped_daily_bars test
    @test stocks_grouped_daily_bars(tabular_opts, "2020-10-14"; adjusted=true) |> length >= 8000
    @test stocks_grouped_daily_bars(regular_opts, "2020-10-14"; adjusted=false) |> length >= 8000

    # previous_close test
    @test stocks_previous_close(tabular_opts, "AAPL") |> length == 1
    @test stocks_previous_close(regular_opts, "AAPL") |> length == 1

    # aggregates_bars test
    @test stocks_aggregates_bars(tabular_opts, "AAPL") |> length == 1
    @test stocks_aggregates_bars(regular_opts, "AAPL") |> length == 1

    # snapshot_all_tickers test
    @test stocks_snapshot_all_tickers(tabular_opts, "AAPL,AMZN") |> length == 2
    @test stocks_snapshot_all_tickers(regular_opts, "AAPL,AMZN") |> length == 2

    # snapshot_ticker test
    @test stocks_snapshot_ticker(tabular_opts, "AAPL") |> length == 9
    @test stocks_snapshot_ticker(regular_opts, "AAPL") |> length == 9

    # snapshot_gainers_losers test
    @test stocks_snapshot_gainers_losers(tabular_opts, "losers") |> length >= 20
    @test stocks_snapshot_gainers_losers(regular_opts, "gainers") |> length >= 20
end

@testset "Crypto API" begin
    # last_trade_crypto_pair test
    @test last_trade_crypto_pair(regular_opts) |> length >= 5
    @test last_trade_crypto_pair(tabular_opts) |> length >= 1

    # crypto_daily_open_close test
    @test crypto_daily_open_close(regular_opts, "BTC", "USD", "2020-10-14"; adjusted=true) |> length >= 7
    @test crypto_daily_open_close(tabular_opts, "BTC", "USD", "2020-10-14"; adjusted=false) |> length >= 7

    # historic_crypto_trades test
    @test historic_crypto_trades(regular_opts, "BTC", "USD", "2020-10-14"; limit=100) |> length >= 7
    @test historic_crypto_trades(tabular_opts, "BTC", "USD", "2020-10-14"; limit=100) |> length >= 7

    # crypto_grouped_daily_bars test
    @test crypto_grouped_daily_bars(regular_opts, "2020-10-14"; adjusted=true) |> length >= 1
    @test crypto_grouped_daily_bars(tabular_opts, "2020-10-14"; adjusted=false) |> length >= 1

    # crypto_previous_close test
    @test crypto_previous_close(regular_opts, "X:BTCUSD"; adjusted=true) |> length == 1
    @test crypto_previous_close(tabular_opts, "X:BTCUSD"; adjusted=false) |> length == 1

    # crypto_aggregates_bars test
    @test crypto_aggregates_bars(regular_opts, "X:BTCUSD", 1, "day", "2020-10-14", "2020-10-14"; adjusted=true) |> length == 1
    @test crypto_aggregates_bars(tabular_opts, "X:BTCUSD", 1, "day", "2020-10-14", "2020-10-14"; adjusted=false) |> length == 1

    # crypto_snapshot_all_tickers test
    @test crypto_snapshot_all_tickers(regular_opts) |> length >= 1
    @test crypto_snapshot_all_tickers(tabular_opts) |> length >= 1

    # crypto_snapshot_ticker test
    @test crypto_snapshot_ticker(regular_opts, "X:BTCUSD") |> length >= 1
    @test crypto_snapshot_ticker(tabular_opts, "X:BTCUSD") |> length >= 1

    # crypto_snapshot_ticker_full_book
    @test crypto_snapshot_ticker_full_book(regular_opts, "X:BTCUSD") |> length >= 1
    @test crypto_snapshot_ticker_full_book(tabular_opts, "X:BTCUSD") |> length >= 1

    # crypto_snapshot_gainers_losers test
    @test crypto_snapshot_gainers_losers(regular_opts, "gainers") |> length >= 1
    @test crypto_snapshot_gainers_losers(tabular_opts, "losers") |> length >= 1
end

@testset "Forex API" begin
    # historic_forex_ticks test
    @test historic_forex_ticks(regular_opts, "AUD", "USD", "2020-10-14"; limit=100) |> length >= 1
    @test historic_forex_ticks(tabular_opts, "AUD", "USD", "2020-10-14"; limit=100) |> length >= 1

    # real_time_currency_conversion test
    @test real_time_currency_conversion(regular_opts, "AUD", "USD") |> length >= 8
    @test real_time_currency_conversion(tabular_opts, "AUD", "USD") |> length >= 1

    # last_quote_currency_pair test
    @test last_quote_currency_pair(regular_opts, "AUD", "USD") |> length >= 1
    @test last_quote_currency_pair(tabular_opts, "AUD", "USD") |> length >= 1

    # forex_grouped_daily_bars test
    @test forex_grouped_daily_bars(regular_opts, "2020-10-14"; adjusted=true) |> length >= 1
    @test forex_grouped_daily_bars(tabular_opts, "2020-10-14"; adjusted=false) |> length >= 1

    # forex_previous_close test
    @test forex_previous_close(regular_opts, "C:EURUSD"; adjusted=true) |> length >= 1
    @test forex_previous_close(tabular_opts, "C:EURUSD"; adjusted=false) |> length >= 1

    # forex_aggregates_bars test
    @test forex_aggregates_bars(regular_opts, "C:EURUSD", 1, "day", "2020-10-14", "2020-10-14"; adjusted=true) |> length >= 1
    @test forex_aggregates_bars(tabular_opts, "C:EURUSD", 1, "day", "2020-10-14", "2020-10-14"; adjusted=false) |> length >= 1

    # forex_snapshot_ticker test
    @test forex_snapshot_ticker(regular_opts, "C:EURUSD") |> length >= 1
    @test forex_snapshot_ticker(tabular_opts, "C:EURUSD") |> length >= 1

    # forex_snapshot_all_tickers test
    @test forex_snapshot_all_tickers(regular_opts) |> length >= 1
    @test forex_snapshot_all_tickers(tabular_opts) |> length >= 1

    # forex_snapshot_gainers_losers test
    @test forex_snapshot_gainers_losers(regular_opts, "gainers") |> length >= 1
    @test forex_snapshot_gainers_losers(tabular_opts, "losers") |> length >= 1
end


text:
using PyPlot
using Optim

function PlotNoise_VsR(en,in;YLims = [1e-10,1e-8])
    k_b = 1.38e-23
    R_L_List = 1:1:1e5
    noise = zeros(length(R_L_List),1)
    Thnoise = zeros(length(R_L_List),1)
    NoiseZ = en/in

    for i in 1:length(R_L_List)
        noise[i] = Noise(R_L_List[i],en,in)
        Thnoise[i] = ThermalNoise(R_L_List[i])
    end
    display(PyPlot.plot(R_L_List,noise[:],label="Total Noise"))
    PyPlot.plot(R_L_List,Thnoise[:],"r:",label="Thermal Noise Only")
    PyPlot.plot(NoiseZ,Noise(NoiseZ,en,in),"k^",label="Noise Matched Location",ms=10)
    ylim(YLims[1],YLims[2])
    legend()
    xlabel("Termination Resistance(Ω)")
    ylabel("Input Noise (V/√Hz)")
    xscale("log")
    yscale("log")
end

function PlotNoiseFigure_VsR(en,in;YLims = [1 10])
    k_b = 1.38e-23
    R_L_List = 1:1:1e5
    noise = zeros(length(R_L_List),1)
    NF = zeros(length(R_L_List),1)
    Thnoise = zeros(length(R_L_List),1)
    NoiseZ = en/in

    for i in 1:length(R_L_List)
        noise[i] = Noise(R_L_List[i],en,in)
        Thnoise[i] = ThermalNoise(R_L_List[i])

    end
    NF = 1 .+ noise[:] ./ Thnoise[:]

    display(PyPlot.plot(R_L_List,NF,label="Total Noise"))

    ylim(YLims[1],YLims[2])
    legend()
    xlabel("Termination Resistance(Ω)")
    ylabel("Noise Figure (NF)")

end

"""
This function takes in a vector of measured noise (in V/√Hz) and the list of resistors that were used for termination, and estimates the preamp's en,in
    Inputs: RVals, vector of resistances
            MeasNoise, vector of the same length of measured noise values V/√Hz at the output.

"""
function DetermineNoise(RVals::Vector,MeasNoise::Vector)
    ## This function takes in a vector of measured noise (in V/√Hz) and the
    ## List of resistors that were used for termination, and estimates the preamp's en,in
    Error(NoiseParams) = 1e6*sum((MeasNoise - Noise(RVals,NoiseParams[1],NoiseParams[2])).^2)
    OptimResults = optimize(Error,[0.0,0.0],g_tol=1e-13)
    OptimResults.minimizer[1]
    PlotNoise_VsR(OptimResults.minimizer[1],OptimResults.minimizer[2])
    PyPlot.plot(RVals,MeasNoise,"ro",label = "Measured",ms=5)
    legend()
    return (OptimResults.minimizer)
end

function Noise(R_L,en,in;T=300)
    k_b = 1.38e-23
     sqrt.(4*k_b*T*R_L) .+ en .+ in*R_L
 end

function ThermalNoise(R_L)
    k_b = 1.38e-23
     sqrt.(4*k_b*300*R_L)
end


text:
include("../LF_common/m_LF1d.jl")
include("../LF_common/m_LF3d.jl")
include("../LF_common/m_Gvectors.jl")

include("../LF_common/apply_Laplacian.jl")
include("../LF_common/solve_poisson_cg.jl")
include("linsolve_cg_v2.jl")
include("linsolve_cr.jl")
include("solve_poisson_FFT.jl")

function test_main( NN::Array{Int64} )

  AA = [0.0, 0.0, 0.0]
  BB = [16.0, 16.0, 16.0]

  # Initialize LF
  LF = init_LF3d_p( NN, AA, BB, verbose=true )

  # Box dimensions
  Lx = BB[1] - AA[1]
  Ly = BB[2] - AA[2]
  Lz = BB[3] - AA[3]

  # Center of the box
  x0 = Lx/2.0
  y0 = Ly/2.0
  z0 = Lz/2.0

  # Parameters for two gaussian functions
  sigma1 = 0.75
  sigma2 = 0.50
  Npoints = LF.Nx * LF.Ny * LF.Nz

  rho = zeros( Npoints )
  phi = zeros( Npoints )

  # Initialization of charge density
  for ip = 1 : Npoints
    r = norm( LF.lingrid[:,ip] - [ x0, y0, z0 ] )
    rho[ip] = exp( -r^2 / (2.0*sigma2^2) ) / (2.0*pi*sigma2^2)^1.5 -
              exp( -r^2 / (2.0*sigma1^2) ) / (2.0*pi*sigma1^2)^1.5
  end

  deltaV = LF.LFx.h * LF.LFy.h * LF.LFz.h

  @printf("#Test norm charge: %f\n", sum(rho)*deltaV)
  @printf("#Solving Poisson equation:\n")
  #phi = solve_poisson_cg( LF, -4.0*pi*rho, 1000, verbose=true, TOL=5e-9 )

  phi = zeros(Npoints)
  linsolve_cg_v2!( LF, phi, -4.0*pi*rho, verbose=true, TOL=5e-9 )
  #linsolve_cr!( LF, phi, -4.0*pi*rho, verbose=true, TOL=5e-9 )

  #phi = linsolve_bicgstab( LF, -4.0*pi*rho, 200, verbose=true )

  Gvec    = GvectorsT( NN, diagm( [Lx,Ly,Lz] ) )
  phi_FFT = solve_poisson_FFT( Gvec, rho )

  @printf("LF ΔV = %18.10f\n", deltaV)
  @printf("Gv ΔV = %18.10f\n", Gvec.Ω/prod(NN))

  # Calculation of Hartree energy
  Unum = 0.5*sum( rho .* phi ) * deltaV
  Unum_FFT = 0.5*sum( rho .* phi_FFT ) * deltaV
  Uana = ( ( 1.0/sigma1 + 1.0/sigma2 ) / 2.0 - sqrt(2.0)/sqrt(sigma1^2 + sigma2^2) ) / sqrt(pi)
  @printf("#CG, analytic, diff: %20.10f %20.10f %20.10f\n", Unum, Uana, abs(Unum-Uana))
  @printf("#PW, analytic, diff: %20.10f %20.10f %20.10f\n", Unum_FFT, Uana, abs(Unum_FFT-Uana))
end

#@code_native test_main([3,3,3])
@time test_main([25,25,25])


text:
#=

MuSSE equations

Ignacio Quintero Mächler

t(-_-t)

September 16 2017

=#

"""
    make_musse(t,u,p,du)

MuSSE equation function for integration.
"""
function make_musse(k::Int64)

  Qf! = makeQ(k)

  # preallocate matrices & vectors
  Q ::Array{Float64,2} = zeros(k, k)
  Qs::Array{Float64,1} = zeros(k)
  Qe::Array{Float64,1} = zeros(k)

  # estimate indices previously
  qsi::UnitRange{Int64} = 1:k
  qei::UnitRange{Int64} = (k+1):2k

  function f(du, u, p, t)

    Qf!(Q, p)
    BLAS.gemv!('N', 1.0, Q, u[qsi], 0.0, Qs)
    BLAS.gemv!('N', 1.0, Q, u[qei], 0.0, Qe)

    for i in Base.OneTo(k)
      du[i]   =        - (p[i] + p[k+i])*u[i]   + 2.0*p[i]*u[i+k]*u[i]   + Qs[i]
      du[i+k] = p[k+i] - (p[i] + p[k+i])*u[i+k] +     p[i]*u[i+k]*u[i+k] + Qe[i]
    end

  end

  return f
end

"""
    makeQ(k::Int64)

Transform instantaneous transition rates into the instantaneous
transition matrix `Q`.
"""
function makeQ(k::Int64)

  pidx::UnitRange{Int64}  = (2k+1):(k + k*k)
  Didx::StepRange{Int64}  = 1:(k+1):(k*k)
  Qidx::Array{Int64,1}    = setdiff(1:k*k, Didx)
  pitr::Base.OneTo{Int64} = Base.OneTo(lastindex(pidx))
  kitr::Base.OneTo{Int64} = Base.OneTo(k)

  function f(Q::Array{Float64,2}, p::Array{Float64,1})

    @inbounds begin
      # assign parameters
      @simd for i in pitr
        Q[Qidx[i]] = p[pidx[i]]
      end

      #assign 0s
      for i in Didx
        Q[i] = 0.0
      end

      #row sums
      for j = kitr, i = kitr
        if i == j
          continue
        end
        Q[Didx[i]] -= Q[i,j]
      end
    end

  return nothing
  end
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
function compress_then_write(b::AbstractVector{T}, io) where {T}
    compress_then_write(Vector(b), io)
end

# the generic dispatch for T where is isbits
function compress_then_write(b::Vector{T}, io) where {T}
    bbc = Blosc.compress(b)
    res = length(bbc)
    write(io, bbc)
    return (len = res, type = T)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
type HestonModel{CalibratedModelType} <: Model{CalibratedModelType}
  modT::CalibratedModelType
  theta::ConstantParameter
  kappa::ConstantParameter
  sigma::ConstantParameter
  rho::ConstantParameter
  v0::ConstantParameter
  process::HestonProcess
  common::ModelCommon
end

function HestonModel(process::HestonProcess)
  theta = ConstantParameter([process.theta], PositiveConstraint())
  kappa = ConstantParameter([process.kappa], PositiveConstraint())
  sigma = ConstantParameter([process.sigma], PositiveConstraint())
  rho = ConstantParameter([process.rho], BoundaryConstraint(-1.0, 1.0))
  v0 = ConstantParameter([process.v0], PositiveConstraint())

  hm = HestonModel(CalibratedModelType(), theta, kappa, sigma, rho, v0, process, ModelCommon())

  generate_arguments!(hm)

  return hm
end

get_theta(hm::HestonModel) = hm.theta.data[1]
get_kappa(hm::HestonModel) = hm.kappa.data[1]
get_sigma(hm::HestonModel) = hm.sigma.data[1]
get_rho(hm::HestonModel) = hm.rho.data[1]
get_v0(hm::HestonModel) = hm.v0.data[1]

generate_arguments!(hm::HestonModel) =
  hm.process = HestonProcess(hm.process.riskFreeRate, hm.process.dividendYield, hm.process.s0,
              get_v0(hm), get_kappa(hm), get_theta(hm), get_sigma(hm), get_rho(hm))


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior.Gherkin: issuccessful, parsescenario!, Given, When, Then, ByLineParser, ScenarioStep, ParseOptions

@testset "Data tables          " begin
    @testset "A Scenario with a data table; The data table is associated with the step" begin
        text = """
        Scenario: Data tables
            When some action
            Then some tabular data
                | header 1 | header 2 |
                | foo 1    | bar 1    |
                | foo 2    | bar 2    |
                | foo 3    | bar 3    |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps[2].datatable == [
            ["header 1", "header 2"],
            ["foo 1", "bar 1"],
            ["foo 2", "bar 2"],
            ["foo 3", "bar 3"],
        ]
    end

    @testset "A data table with comments in the elements; Comments are ignored" begin
        text = """
        Scenario: Data tables
            When some action
            Then some tabular data
                | header 1 | header 2 |
                | foo 1    | bar 1    |
                | foo 2    | bar 2    |
                # Comment
                | foo 3    | bar 3    |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps[2].datatable == [
            ["header 1", "header 2"],
            ["foo 1", "bar 1"],
            ["foo 2", "bar 2"],
            ["foo 3", "bar 3"],
        ]
    end

    @testset "A data table with blank lines; Blank lines are ignored" begin
        text = """
        Scenario: Data tables
            When some action
            Then some tabular data
                | header 1 | header 2 |
                | foo 1    | bar 1    |

                | foo 2    | bar 2    |

                | foo 3    | bar 3    |
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps[2].datatable == [
            ["header 1", "header 2"],
            ["foo 1", "bar 1"],
            ["foo 2", "bar 2"],
            ["foo 3", "bar 3"],
        ]
    end

    @testset "The data table step is followed by another step; The data table still has four rows" begin
        text = """
        Scenario: Data tables
            When some action
            Then some tabular data
                | header 1 | header 2 |
                | foo 1    | bar 1    |
                | foo 2    | bar 2    |
                | foo 3    | bar 3    |
             And some other step
        """

        byline = ByLineParser(text)
        result = parsescenario!(byline)

        @test issuccessful(result)
        scenario = result.value

        @test scenario.steps[2].datatable == [
            ["header 1", "header 2"],
            ["foo 1", "bar 1"],
            ["foo 2", "bar 2"],
            ["foo 3", "bar 3"],
        ]
    end
end


text:
macro test_elements(name, μ, ele, pos, vel)
    quote
        @testset $name begin
            μ = $μ
            a, e, i, Ω, ω, ν = $ele
            pos = $pos
            vel = $vel

            pos₁, vel₁ = cartesian(a, e, i, Ω, ω, ν, μ)
            a₁, e₁, i₁, Ω₁, ω₁, ν₁ = keplerian(pos, vel, μ)

            @testset "Cartesian" begin
                @testset for i = 1:3
                    @test pos₁[i] ≈ pos[i]
                    @test vel₁[i] ≈ vel[i]
                end
            end

            @testset "Keplerian" begin
                @test a ≈ a₁
                @test e ≈ e₁ atol=sqrt(eps())
                @test i ≈ i₁
                @test Ω ≈ Ω₁
                @test ω ≈ ω₁
                @test ν ≈ ν₁
            end
        end
    end
end

macro test_elements(name, μ, ele)
    quote
        @testset $name begin
            μ = $μ
            a, e, i, Ω, ω, ν = $ele

            pos₁, vel₁ = cartesian(a, e, i, Ω, ω, ν, μ)
            a₁, e₁, i₁, Ω₁, ω₁, ν₁ = keplerian(pos₁, vel₁, μ)

            @test a ≈ a₁
            @test e ≈ e₁ atol=sqrt(eps())
            @test i ≈ i₁
            @test Ω ≈ Ω₁
            @test ω ≈ ω₁
            @test ν ≈ ν₁
        end
    end
end

@testset "Two Body" begin
    @testset "Orbital Elements" begin
        @testset "Perifocal" begin
            p = 1.13880762905224e7
            ecc = 0.7311
            ν = 0.44369564302687126
            μ = 3.9860047e14

            pos = [6194863.12535486, 2944437.90016286, 0.0]
            vel = [-2539.71254827, 9668.69568539, 0.0]
            pos₁, vel₁ = perifocal(p, ecc, ν, μ)
            @testset for i = 1:3
                @test pos₁[i] ≈ pos[i]
                @test vel₁[i] ≈ vel[i]
            end
        end
        # Reference values from Orekit
        @test_elements("Elliptic",
                       3.9860047e14,
                       [24464560.0, 0.7311, 0.122138, 1.00681, 3.10686, 0.44369564302687126],
                       [-0.107622532467967e+07, -0.676589636432773e+07, -0.332308783350379e+06],
                       [0.935685775154103e+04, -0.331234775037644e+04, -0.118801577532701e+04])
        # Reference values from poliastro
        @test_elements("Circular",
                       3.986004418e14,
                       [6778136.6, 0.0, deg2rad(15), deg2rad(20), 0.0, deg2rad(30)],
                       [4396398.60746266, 5083838.45333733,  877155.42119322],
                       [-5797.06004014,  4716.60916063,  1718.86034246])
        @test_elements("Circular Orekit",
                       3.9860047e14,
                       # Orekit uses an argument of pericenter of 3.10686 here
                       [24464560.0, 0.0, 0.122138, 1.00681, 0.0, 0.048363])
        @test_elements("Hyperbolic Orekit",
                       3.9860047e14,
                       [-24464560.0, 1.7311, 0.122138, 1.00681,  3.10686, 0.12741601769795755])
        @test_elements("Equatorial",
                       3.9860047e14,
                       [24464560.0, 0.7311, 0.0, 0.0, 3.10686, 0.44369564302687126])
        @test_elements("Circular-Equatorial",
                       3.9860047e14,
                       [24464560.0, 0.0, 0.0, 0.0, 0.0, 0.44369564302687126])
    end
    @testset "Anomalies" begin
        anomalies = -π : 0.1π : π
        anomalies2π = 0.0 : 0.1π : 2π
        @testset "Elliptic" begin
            eccentricities = 0.1:0.1:0.9
            @testset "E <-> ν" begin
                @testset for E = anomalies, ecc = eccentricities
                    ν = transform(elliptic, eccentric_anomaly, true_anomaly, E, ecc)
                    E₁ = transform(elliptic, true_anomaly, eccentric_anomaly, ν, ecc)
                    @test E₁ ≈ E
                end
                @testset for E = anomalies2π, ecc = eccentricities
                    ν = transform(elliptic, eccentric_anomaly, true_anomaly, E, ecc)
                    E₁ = transform(elliptic, true_anomaly, eccentric_anomaly, ν, ecc)
                    @test E₁ ≈ E
                end
            end
            @testset "E <-> M" begin
                @testset for E = anomalies, ecc = eccentricities
                    M = transform(elliptic, eccentric_anomaly, mean_anomaly, E, ecc)
                    E₁ = transform(elliptic, mean_anomaly, eccentric_anomaly, M, ecc)
                    @test E₁ ≈ E
                end
                @testset for E = anomalies2π, ecc = eccentricities
                    M = transform(elliptic, eccentric_anomaly, mean_anomaly, E, ecc)
                    E₁ = transform(elliptic, mean_anomaly, eccentric_anomaly, M, ecc)
                    @test E₁ ≈ E
                end
            end
            @testset "M <-> ν" begin
                @testset for M = anomalies, ecc = eccentricities
                    ν = transform(elliptic, mean_anomaly, true_anomaly, M, ecc)
                    M₁ = transform(elliptic, true_anomaly, mean_anomaly, ν, ecc)
                    @test M₁ ≈ M
                end
                @testset for M = anomalies2π, ecc = eccentricities
                    ν = transform(elliptic, mean_anomaly, true_anomaly, M, ecc)
                    M₁ = transform(elliptic, true_anomaly, mean_anomaly, ν, ecc)
                    @test M₁ ≈ M
                end
            end
        end
        @testset "Hyperbolic" begin
            eccentricities = 1.1:0.1:10.0
            @testset "E <-> ν" begin
                @testset for E = anomalies, ecc = eccentricities
                    ν = transform(hyperbolic, eccentric_anomaly, true_anomaly, E, ecc)
                    E₁ = transform(hyperbolic, true_anomaly, eccentric_anomaly, ν, ecc)
                    @test E₁ ≈ E
                end
                @testset for E = anomalies2π, ecc = eccentricities
                    ν = transform(hyperbolic, eccentric_anomaly, true_anomaly, E, ecc)
                    E₁ = transform(hyperbolic, true_anomaly, eccentric_anomaly, ν, ecc)
                    @test E₁ ≈ E
                end
            end
            @testset "E <-> M" begin
                @testset for E = anomalies, ecc = eccentricities
                    M = transform(hyperbolic, eccentric_anomaly, mean_anomaly, E, ecc)
                    E₁ = transform(hyperbolic, mean_anomaly, eccentric_anomaly, M, ecc)
                    @test E₁ ≈ E
                end
                @testset for E = anomalies2π, ecc = eccentricities
                    M = transform(hyperbolic, eccentric_anomaly, mean_anomaly, E, ecc)
                    E₁ = transform(hyperbolic, mean_anomaly, eccentric_anomaly, M, ecc)
                    @test E₁ ≈ E
                end
            end
            @testset "M <-> ν" begin
                @testset for M = anomalies, ecc = eccentricities
                    ν = transform(hyperbolic, mean_anomaly, true_anomaly, M, ecc)
                    M₁ = transform(hyperbolic, true_anomaly, mean_anomaly, ν, ecc)
                    @test M₁ ≈ M
                end
                @testset for M = anomalies2π, ecc = eccentricities
                    ν = transform(hyperbolic, mean_anomaly, true_anomaly, M, ecc)
                    M₁ = transform(hyperbolic, true_anomaly, mean_anomaly, ν, ecc)
                    @test M₁ ≈ M
                end
            end
        end
        @testset "Parabolic" begin
            @testset "E <-> ν" begin
                @testset for E = anomalies
                    ν = transform(parabolic, eccentric_anomaly, true_anomaly, E, 1.0)
                    E₁ = transform(parabolic, true_anomaly, eccentric_anomaly, ν, 1.0)
                    @test E₁ ≈ E
                end
                @testset for E = anomalies2π
                    ν = transform(parabolic, eccentric_anomaly, true_anomaly, E, 1.0)
                    E₁ = transform(parabolic, true_anomaly, eccentric_anomaly, ν, 1.0)
                    @test E₁ ≈ E
                end
            end
            @testset "E <-> M" begin
                @testset for E = anomalies
                    M = transform(parabolic, eccentric_anomaly, mean_anomaly, E, 1.0)
                    E₁ = transform(parabolic, mean_anomaly, eccentric_anomaly, M, 1.0)
                    @test E₁ ≈ E
                end
                @testset for E = anomalies2π
                    M = transform(parabolic, eccentric_anomaly, mean_anomaly, E, 1.0)
                    E₁ = transform(parabolic, mean_anomaly, eccentric_anomaly, M, 1.0)
                    @test E₁ ≈ E
                end
            end
            @testset "M <-> ν" begin
                @testset for M = anomalies
                    ν = transform(parabolic, mean_anomaly, true_anomaly, M, 1.0)
                    M₁ = transform(parabolic, true_anomaly, mean_anomaly, ν, 1.0)
                    @test M₁ ≈ M
                end
                @testset for M = anomalies2π
                    ν = transform(parabolic, mean_anomaly, true_anomaly, M, 1.0)
                    M₁ = transform(parabolic, true_anomaly, mean_anomaly, ν, 1.0)
                    @test M₁ ≈ M
                end
            end
        end
    end
    @testset "Kepler" begin
        # Source: Vallado, Fundamentals of Astrodynamics and Applications, 4th edition, p. 94-95
        μ = 3.986004418e5
        r₀exp = [1131.340, -2282.343, 6672.423]
        v₀exp = [-5.64305, 4.30333, 2.42879]
        Δt = 40*60
        r₁exp = [-4219.7527, 4363.0292, -3958.7666]
        v₁exp = [3.689866, -1.916735, -6.112511]
        r₁, v₁ = kepler(μ, r₀exp, v₀exp, Δt)
        @test r₁ ≈ r₁exp rtol=1e-4
        @test v₁ ≈ v₁exp rtol=1e-5
        r₀, v₀ = kepler(μ, r₁, v₁, -Δt)
        @test r₀ ≈ r₀exp
        @test v₀ ≈ v₀exp
        r₁, v₁ = kepler(μ, r₀exp, v₀exp, eps())
        @test r₁ == r₀exp
        @test v₁ == v₀exp
        @test_throws ErrorException kepler(μ, r₀exp, v₀exp, Δt, 1)
    end
end


text:
# Всякие вспомогательные функции

@doc raw"""
    shishkin_mesh(a::Real, b::Real,
                  x_tp::Real, ε::Real,
                  N::Int = 50,
                  C_i::Real = 1.0, K_i::Real = 1.0,
                  C_b::Real = 1.0, K_b::Real = 1.0) --> Vector

Возвращает кусочнораномерную сетку со сгущением на границах и
на переходном слое.

Экскиз:
```md
.....  .  .  .  .  .  .  .  .  .  .......  .  .  .  .  .  .  .  .  .  .  ......
```

# Arguments
- `a::Int`      Левая граница.
- `b::Int`      Правая граница.
- `x_tp::Real`  Положение переходного слоя.
- `ε::Real`     Малый параметр при старшей производной.
- `C_i::Real`   Масштабирующий коэффициент для ширины внутреннего сгущения.
- `K_i::int`    Масштабирующий коэффициент количества интервалов внутреннего сгущения.
- `C_b::Real`   Масштабирующий коэффициент для ширин пограничных сгущений.
- `K_b::Int`    Масштабирующий коэффициент количеств интервалов пограничных сгущений.
- `K::Int`      Кол-во интервалов вне всех сгущений.
"""
function shishkin_mesh(a::Real, b::Real,
                       x_tp::Real, ε::Real,
                       N::Int = 50,
                       C_i::Real = 1.0, K_i::Real = 1.0,
                       C_b::Real = 1.0, K_b::Real = 1.0)

    if ! ( a < x_tp < b )
        throw(ArgumentError("``a, x_tp, b`` должна образовывать строго упорядоченную тройку: $((a, x_tp, b))"))
    end

    N_i = Int(round(K_i * N));
    N_b = Int(round(K_b * N));
    x = zeros(N + N_i + 2*N_b + 1);

    if all( x -> x<10, [N_i, N_b, N] )
        throw(ArgumentError("Количество точек в сгущениях или вне их меньше 10. " *
                            "Используйте другую комбинацию базового количества точек `N` и сгущающих коэффициентоа `K_i, K_b` " *
                            "\n (N, N_b, N_i) == $((N, N_b, N_i))"));
    end

        # Вычисляем долю узлов от N, которые будут лежать слева от внутреннего слоя
        # Проверяем, не пересекаются ли "зоны влияния" внутреннего и пограничного слоёв
        if ((x_tp - abs(C_i*ε*log(ε))) - (a + abs(C_b*ε*log(ε)))) >= 0

            if ((b - abs(C_b*ε*log(ε))) - (x_tp + abs(C_i*ε*log(ε)))) >= 0
                k_left = ((x_tp - abs(C_i*ε*log(ε))) - (a + abs(C_b*ε*log(ε))))/(b - a - 2*abs(C_i*ε*log(ε)) - 2*abs(C_b*ε*log(ε)));
            else
                k_left = 1;
            end
        else
            k_left = 0;
        end
        # Вычисляем число узлов, которые буду лежать слева и справа от внутреннего слоя
        N_left = Int(round(k_left*N));
        N_right = N - N_left;

        # Вычисляем соответствующие шаги

        # Проверяем, не пересекаются ли "зоны влияния" внутреннего и пограничного слоёв
        # В разных случаях действуем по-разному.
        if ((x_tp - abs(C_i*ε*log(ε))) - (a + abs(C_b*ε*log(ε)))) >= 0
            if ((b - abs(C_b*ε*log(ε))) - (x_tp + abs(C_i*ε*log(ε)))) >= 0
                h_int = 2*abs(C_i*ε*log(ε))/(N_i);
                h_left = ((x_tp - abs(C_i*ε*log(ε))) - (a + abs(C_b*ε*log(ε))))/N_left;
                h_right = ((b - abs(C_b*ε*log(ε))) - (x_tp + abs(C_i*ε*log(ε))))/N_right;
                h_bound = abs(C_b*ε*log(ε))/(N_b);
            else
                h_bound = abs(C_b*ε*log(ε))/(N_b);
                h_left = ((x_tp - abs(C_i*ε*log(ε))) - (a + abs(C_b*ε*log(ε))))/N_left;
                h_int = (b - (x_tp - abs(C_i*ε*log(ε))))/(N_i + N_b);
            end
        else
            h_int = ((x_tp + abs(C_i*ε*log(ε))) - a)/(N_b + N_i);
            h_right = ((b - abs(C_b*ε*log(ε))) - (x_tp + abs(C_i*ε*log(ε))))/N_right;
            h_bound = abs(C_b*ε*log(ε))/(N_b);
        end

        x[1] = a;

        # Проверяем, не пересекаются ли "зоны влияния" внутреннего и пограничного слоёв
        # В разных случаях действуем по-разному.

        if ((x_tp - abs(C_i*ε*log(ε))) - (a + abs(C_b*ε*log(ε)))) >= 0
            if ((b - abs(C_b*ε*log(ε))) - (x_tp + abs(C_i*ε*log(ε)))) >= 0
                for n = 1:(N_b)
                    x[n + 1] = x[n] + h_bound;
                end
                for n = (N_b + 1):(N_b + N_left)
                    x[n + 1] = x[n] + h_left;
                end
                for n = (N_b + N_left + 1):(N_b + N_left + N_i)
                    x[n + 1] = x[n] + h_int;
                end
                for n = (N_b + N_left + N_i + 1):(N_b + N + N_i)
                    x[n + 1] = x[n] + h_right;
                end
                for n = (N_b + N + N_i + 1):(N + N_i + 2*N_b)
                    x[n + 1] = x[n] + h_bound;
                end
            else
                for n = 1:(N_b)
                    x[n + 1] = x[n] + h_bound;
                end
                for n = (N_b + 1):(N_b + N_left)
                    x[n + 1] = x[n] + h_left;
                end
                for n = (N_b + N_left + 1):(N + N_i + 2*N_b)
                    x[n + 1] = x[n] + h_int;
                end
            end
        else
            for n = 1:(N_b + N_left + N_i)
                x[n + 1] = x[n] + h_int;
            end
            for n = (N_b + N_left + N_i + 1):(N_b + N + N_i)
                x[n + 1] = x[n] + h_right;
            end
            for n = (N_b + N + N_i + 1):(N + N_i + 2*N_b)
                x[n + 1] = x[n] + h_bound;
            end
        end

    return x
end

@doc raw"""
    strip_borderPoints(a::Vector, N) -> Vector

Функция для внутренного использования.

Входящий массив должен быть размера `N+1`.
Обрезает граничные точки слева и справа.
Возвращает массив размера `N-1`.

# Example
```jldoctest
julia> N = 10; a = collect(1:N+1);

julia> NonLinearReactionAdvectionDiffusionWithFrontData.strip_borderPoints(a, N)
9-element Array{Int64,1}:
  2
  3
  4
  5
  6
  7
  8
  9
 10
```
"""
function strip_borderPoints(a::Vector, N)
    @assert length(a) == N+1 "Размерность входного вектора отличается от $(N+1)"
    return a[2:N]
end

function strip_borderPoints(a::Matrix, N)
    @assert size(a,1) == N+1 "Размерность входного вектора отличается от $(N+1)"
    return a[2:N,:]
end

@doc raw"""
    u_init(x::Real; ε::Real = 0.2, x_tp = 0.15) -> Real

Начальные условия в виде $(x^2 - x -2) -6 \tanh( -3 \xi)$, где $\xi = \frac{x - x_tp}{ε}$.
`x_tp` — положение внутрннего переходного слоя.

!!! note
    Граничные условия для этих начальных условий должны быть заданы как `(-8, 4)`.

!!! note
    Вы вольны устанавливать начальные условия на свое усмотрения, эта функция внесена в модуль для удобства повествования.

# Example
```jldoctest
julia> NonLinearReactionAdvectionDiffusionWithFrontData.u_init.(0:0.1:1)
11-element Array{Float64,1}:
 -7.868156688432881
 -5.900893714323723
  1.650893714323725
  3.6581566884328813
  3.753366656356917
  3.749669571706625
  3.7599835485135173
  3.7899991809276505
  3.839999959220786
  3.909999997969723
  3.999999999898918
```
"""
function u_init(x::Real; ε::Real = 0.2, x_tp = 0.15)
    ξ = (x - x_tp) / ε;
    return ((x^2 - x - 2) - 6 * tanh(-3 * ξ))
end

@doc raw"""
    delta(x, Xₙ, x₀ = 0) -> ::Real

Вычисляет конечно разностную аппроксимацию дельта функции ``\delta(x; x₀)``
на сетке `Xₙ` исходя из ``\int\limits_a^b \delta(x; x₀) dx = 1``:
если ``x, x_0 \in [X_n, X_{n+1}]``, возвращаем ``1/h``, где ``h`` --- шаг сетки.
"""
function delta(x::Real, Xₙ::Vector, x₀::Real = 0)
    if ! (Xₙ[1] <= x <= Xₙ[end])
        throw(DomainError("`x` находится вне `Xₙ`"))
    end
    if ! (Xₙ[1] <= x₀ < Xₙ[end])
        throw(DomainError("`x₀` находится вне полуоткрытого отрезка "*
                          "``[Xₙ[1], Xₙ[end])`` [$(Xₙ[1]), $(Xₙ[end]))"))
    end

    N = length(Xₙ) - 1
    out = 0.0;

    for n in 1:N
        if ( x >= Xₙ[n]) && ( x < Xₙ[n+1] )
            if ( x₀ >= Xₙ[n] ) && ( x₀ < Xₙ[n+1])
                out = 1 / ( Xₙ[n+1] - Xₙ[n] )
            end
        end
    end

    return out;
end

@doc raw"""
    deltaw(n::Int, x₀::Real, Xₙ::Vector, N::Int, w::Real) --> ::Real

Вычисляет ``\delta(x - x₀)``, где `x = Xₙ[n]`.
`n` — номер узла в сетке `Xₙ`(без граничных точек).

Использует более точную конечно разностную аппроксимацию дельта функции,
`w` — априорный параметр, см. [`δw`](@ref).
"""
function deltaw(n::Int, x₀::Real, Xₙ::Vector, N::Int, w::Real)

    @assert length(Xₙ) == N-1
    @assert 1 <= n <= N-1 "`n` вне корректного отрезка"

    out = δw(Xₙ[n] - x₀, w)

    return out;
end

@doc raw"""
    δw(x::Real, ω::Real) -> Real

Линейная аппроксимация ``δ(x)``, `ω` — эвристический парамет.
Подбирается так, чтобы невязка сопряженной задачи
``2 \delta( x - f_1(t)) ( u(x,t) - f_2(t))`` при подстановке
в неё решения `u(x,t;q)` при истинном `q`,
обнулялась почти везде, **но не всюду**.
"""
function δw(x::Real, ω::Real)
    if abs(x/ω) <= 1
        return ( 1 - abs(x/ω))/ω
    else
        return 0.0
    end
end

@doc raw"""
    dparams()

Возвращает стандартный набор параметров для равномерной сетки.

# Return
    `a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀ `
"""
function dparams(;x_tp = 0.1, qfunc = x -> sin(3pi*x),
                 ε = 0.2,
                 a = 0.0, b = 1.0,
                 t_start = 0.0, T_end = 0.36,
                 Nx = 50, Mt = 100,
                 ul_func = x-> -8.0, ur_func = x -> 4.0)

    qf = qfunc;
    h = (b-a)/Nx;
    τ = (T_end-t_start)/Mt;
    Xₙ = [a  + n*h for n in 0:Nx];
    Tₘ = [t_start + m*τ for m in 0:Mt];
    qₙ =    qf.(Xₙ);
    ulₘ=    ul_func.(Tₘ);
    urₘ=    ur_func.(Tₘ);
    u₀ = u_init.(Xₙ, ε = ε, x_tp = x_tp);

    return a, b, t_start, T_end, Nx, Mt, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀

end

@doc raw"""
    dparams_nonuniform()

Возвращает стандартный набор параметров для неравномерной сетки.

# Return
    return a, b, t₀, T, N, M, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀, meshf;
"""
function dparams_nonuniform(;x_tp = 0.1, qfunc = x -> sin(3pi*x),
                 ε = 0.2,
                 a = 0.0, b = 1.0,
                 t_start = 0.0, T_end = 0.36,
                 Nx = 50,
                 Mt = 100,
                 ul_func = x-> -8.0, ur_func = x -> 4.0
                           )

    qf = qfunc;       # Коэффициент линейного усиления
    h = (b-a)/Nx;
    τ = (T_end-t_start)/Mt;
    Tₘ = [t_start + m*τ for m in 0:Mt];   # Сетка по Т
    ulₘ=    ul_func.(Tₘ);
    urₘ=    ur_func.(Tₘ);

    # Замыкание функции формирования сетки по положению переходного слоя
    meshf(x_tp) = shishkin_mesh(a, b, x_tp, ε, Nx, 0.5, 1.0, 0.5, 1.0);

    Xₙ  = meshf(x_tp); ;                 # Сетка по Х
    N   = length(Xₙ) - 1                 # Кол-во интервалов
    qₙ  = qf.(Xₙ);                       # Коэффициент линейного усиления
    u₀  = u_init.(Xₙ, ε=ε, x_tp = x_tp); # Начальные условия

    return a, b, t_start, T_end, N, Mt, ε, Xₙ, Tₘ, qₙ, ulₘ, urₘ, u₀, meshf;
end

@doc raw"""
    heterogeneity_map(XX::Matrix, N::Integer,
                      Tₘ::Vector, M::Integer,
                      u::Matrix,
                      f1_data::Vector, f2_data::Vector,
                      w::Real) --> ::Matrix

Возвращает матрицу неоднородности ``-\delta(x - f_1(t))(u - f_2(t))``.
Следует использовать для контроля выбора параметра аппроксимации ``w`` в
дельта-функции.
"""
function heterogeneity_map(XX::Matrix, N::Integer,
                           Tₘ::Vector, M::Integer,
                           u::Matrix,
                           f1_data::Vector, f2_data::Vector,
                           w::Real)
    out = [ - heterogeneity(n, m, XX[2:N, m], N, Tₘ, M, u[2:N, :], f1_data, f2_data, w) for n in 1:N-1, m in 1:M+1]

    return out
end


text:
function model_referenceregion(
    conv::Function = expconv;
    t::AbstractVector,
    params::NamedTuple,
    crr::AbstractVector,
)
    @extract (rel_kt, kep, kep_rr) params
    ct = rel_kt .* crr .+ (rel_kt * (kep_rr - kep)) .* conv(crr, kep, t)
    return ct
end

function fit_rrm_nls(
    conv::Function = expconv;
    t::AbstractVector,
    crr::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) =
        resolve_fitting_inputs(; t, ca = crr, ct, mask)
    rel_kt, kep, kep_rr = (fill(NaN, volume_size...) for _ = 1:3)
    lls_est = fit_rrm_lls(; t, crr, ct, mask).est
    init_rel_kt, init_kep, init_kep_rr = select(lls_est, (:rel_kt, :kep, :kep_rr))
    model(x, p) = model_referenceregion(
        conv;
        t = x,
        crr,
        params = (rel_kt = p[1], kep = p[2], kep_rr = p[3]),
    )
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        initialvalues = [init_rel_kt[idx], init_kep[idx], init_kep_rr[idx]]
        rel_kt[idx], kep[idx], kep_rr[idx] =
            curve_fit(model, t, ct[idx, :], initialvalues).param
    end
    rel_ve = @. rel_kt * kep_rr / kep
    return (; est = (; rel_kt, rel_ve, kep, kep_rr))
end

function fit_rrm_lls(;
    t::AbstractVector,
    crr::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) =
        resolve_fitting_inputs(; t, ca = crr, ct, mask)
    x1, x2, x3 = (fill(NaN, volume_size...) for _ = 1:3)

    M = zeros(num_timepoints, 3)
    M[:, 1] = crr
    M[:, 2] = cumul_integrate(t, crr, TrapezoidalFast())
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        M[:, 3] = -cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        x1[idx], x2[idx], x3[idx] = M \ ct[idx, :]
    end
    # Apply correction because fit actually returns: [rel_kt, kt/ve_rr, kep]
    rel_kt = x1
    rel_ve = x2 ./ x3
    kep = x3
    kep_rr = x2 ./ x1
    return (; est = (; rel_kt, rel_ve, kep, kep_rr))
end

function fit_crrm_nls(
    conv::Function = expconv;
    t::AbstractVector,
    crr::AbstractVector,
    ct::AbstractArray,
    mask = true,
    kep_rr = 0.0,
)
    (t, ct, mask, num_timepoints, volume_size) =
        resolve_fitting_inputs(; t, ca = crr, ct, mask)
    rel_kt, kep = (fill(NaN, volume_size...) for _ = 1:2)

    if kep_rr <= 0
        rrm_est = fit_rrm_nls(conv; t, crr, ct, mask).est
        viable_est = positive_only_mask(rrm_est)
        kep_rr = interquartile_mean(rrm_est.kep_rr[viable_est])
    end
    lls_est = fit_crrm_lls(; t, crr, ct, kep_rr, mask = mask).est
    init_rel_kt, init_kep = select(lls_est, (:rel_kt, :kep))
    model(x, p) = model_referenceregion(
        conv;
        t = x,
        crr,
        params = (rel_kt = p[1], kep = p[2], kep_rr = kep_rr),
    )
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        initialvalues = [init_rel_kt[idx], init_kep[idx]]
        rel_kt[idx], kep[idx] = curve_fit(model, t, ct[idx, :], initialvalues).param
    end
    rel_ve = @. rel_kt * kep_rr / kep
    return (; est = (; rel_kt, rel_ve, kep, kep_rr))
end

function fit_crrm_lls(;
    t::AbstractVector,
    crr::AbstractVector,
    ct::AbstractArray,
    mask = true,
    kep_rr = 0.0,
)
    (t, ct, mask, num_timepoints, volume_size) =
        resolve_fitting_inputs(; t, ca = crr, ct, mask)
    rel_kt, kep = (fill(NaN, volume_size...) for _ = 1:2)

    if kep_rr <= 0
        rrm_est = fit_rrm_lls(; t, crr, ct, mask).est
        viable_est = positive_only_mask(rrm_est)
        kep_rr = interquartile_mean(rrm_est.kep_rr[viable_est])
    end

    M = zeros(num_timepoints, 2)
    M[:, 1] = crr + kep_rr * cumul_integrate(t, crr, TrapezoidalFast())
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        M[:, 2] = -cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        rel_kt[idx], kep[idx] = M \ ct[idx, :]
    end
    rel_ve = @. rel_kt * kep_rr / kep
    return (est = (rel_kt = rel_kt, rel_ve = rel_ve, kep = kep, kep_rr = kep_rr),)
end

function fit_errm_lls(;
    t::AbstractVector,
    crr::AbstractVector,
    ct::AbstractArray,
    mask = true,
)
    (t, ct, mask, num_timepoints, volume_size) =
        resolve_fitting_inputs(; t, ca = crr, ct, mask)
    x1, x2, x3, x4 = (fill(NaN, volume_size...) for _ = 1:4)

    M = zeros(num_timepoints, 4)
    M[:, 1] = cumul_integrate(t, crr, TrapezoidalFast())
    M[:, 2] = cumul_integrate(t, M[:, 1], TrapezoidalFast())
    M[:, 4] = crr
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        cur_ct = cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        M[:, 3] = -cumul_integrate(t, cur_ct, TrapezoidalFast())
        x1[idx], x2[idx], x3[idx], x4[idx] = M \ cur_ct
    end
    # Apply correction because fit actually returns: [mess, mess, kep, vp/kt_rr]
    a = x1 ./ x4
    b = x2 ./ x4
    root_term = @. a^2 - 4 * b
    @. root_term[root_term<0] = 0
    kep_rr = @. (a - sqrt(root_term)) / 2
    rel_kt = @. x4 * (a - x3 - kep_rr)
    rel_ve = @. rel_kt * kep_rr / x3
    kep = x3
    rel_vp = x4
    return (; est = (; rel_kt, rel_ve, rel_vp, kep, kep_rr))
end

function fit_cerrm_lls(;
    t::AbstractVector,
    crr::AbstractVector,
    ct::AbstractArray,
    mask = true,
    kep_rr = 0.0,
)
    (t, ct, mask, num_timepoints, volume_size) =
        resolve_fitting_inputs(; t, ca = crr, ct, mask)
    x1, rel_vp, kep = (fill(NaN, volume_size...) for _ = 1:3)

    if kep_rr <= 0
        rrm_est = fit_errm_lls(; t, crr, ct, mask).est
        viable_est = positive_only_mask(rrm_est)
        kep_rr = interquartile_mean(rrm_est.kep_rr[viable_est])
    end

    M = zeros(num_timepoints, 3)
    crr_int = cumul_integrate(t, crr, TrapezoidalFast())
    crr_int2 = cumul_integrate(t, crr_int, TrapezoidalFast())
    M[:, 1] = crr_int + kep_rr * crr_int2
    M[:, 2] = crr + kep_rr * crr_int
    for idx in eachindex(IndexCartesian(), mask)
        if mask[idx] == false
            continue
        end
        cur_ct = cumul_integrate(t, ct[idx, :], TrapezoidalFast())
        M[:, 3] = -cumul_integrate(t, cur_ct, TrapezoidalFast())
        x1[idx], rel_vp[idx], kep[idx] = M \ cur_ct
    end
    rel_kt = @. x1 - kep * rel_vp
    rel_ve = @. rel_kt * kep_rr / kep
    return (; est = (; rel_kt, rel_ve, rel_vp, kep, kep_rr))
end

function fit_rrift_with_cerrm(; crr, cp, t, ct, tail_start, kep_rr = 0.0, mask = true)
    cerrm = fit_cerrm_lls(; crr, ct, t, kep_rr, mask).est
    kep_rr = cerrm.kep_rr
    kt_rr = fit_rrift(; t, cp, crr, kep_rr, tail_start)
    ve_rr = kt_rr / kep_rr
    est = relative_to_absolute(cerrm; kt_rr, ve_rr)
    return (est = (; est..., kep_rr, kt_rr, ve_rr))
end

function fit_rrift_with_crrm(; crr, cp, t, ct, tail_start, kep_rr = 0.0, mask = true)
    crrm = fit_crrm_lls(; crr, ct, t, kep_rr, mask).est
    kep_rr = crrm.kep_rr
    kt_rr = fit_rrift(; t, cp, crr, kep_rr, tail_start)
    ve_rr = kt_rr / kep_rr
    est = relative_to_absolute(crrm; kt_rr, ve_rr)
    return (est = (; est..., kep_rr, kt_rr, ve_rr))
end

function fit_rrift(;
    crr::AbstractVector,
    cp::AbstractVector,
    t::AbstractVector,
    tail_start::Int,
    kep_rr::Number,
)
    @assert length(crr) == length(cp) == length(t)
    @assert tail_start < length(crr)
    crr_tail = crr[tail_start:end]
    cp_tail = cp[tail_start:end]
    t_tail = t[tail_start:end]
    numerator =
        crr_tail .- crr_tail[1] .+
        kep_rr .* cumul_integrate(t_tail, crr_tail, TrapezoidalFast())
    denominator = cumul_integrate(t_tail, cp_tail, TrapezoidalFast())
    kt_rr = denominator \ numerator
    return kt_rr
end

function relative_to_absolute(rel_params; kt_rr, ve_rr)
    kt = @. rel_params.rel_kt * kt_rr
    ve = @. rel_params.rel_ve * ve_rr
    kep = rel_params.kep
    if haskey(rel_params, :rel_vp)
        vp = @. rel_params.rel_vp * kt_rr
    else
        vp = zeros(size(kt))
    end
    return (kt = kt, ve = ve, vp = vp, kep = kep)
end


text:
# Script to demonstrate convergence analysis of Ecut values for DFTK.jl

using Atomistic
using AtomsBase
using DFTK
using NBodySimulator
using Unitful
using UnitfulAtomic

setup_threading(n_blas = 4)

N = 8
element = :Ar
box_size = 1.5u"nm" # this number was chosen arbitrarily
reference_temp = 94.4u"K"
thermostat_prob = 0.1 # this number was chosen arbitrarily
Δt = 1e-2u"ps"

initial_bodies = generate_bodies_in_cell_nodes(N, element, box_size, reference_temp)
initial_system = FlexibleSystem(initial_bodies, CubicPeriodicBoundaryConditions(austrip(box_size)))

eq_steps = 20000
eq_thermostat = AndersenThermostat(austrip(reference_temp), thermostat_prob / austrip(Δt))
eq_simulator = NBSimulator(Δt, eq_steps, thermostat = eq_thermostat)
potential = LennardJonesParameters(1.657e-21u"J", 0.34u"nm", 0.765u"nm")

eq_result = @time simulate(initial_system, eq_simulator, potential)

display(@time plot_rdf(eq_result, potential.σ, 0.5))

dftk_potential = DFTKPotential(
    psp = ElementPsp(:Ar, psp = load_psp(list_psp(:Ar, functional = "lda")[1].identifier)),
    lattice = box_size * [[1.0 0 0]; [0 1.0 0]; [0 0 1.0]],
    Ecut = 5u"hartree",
    kgrid = [1, 1, 1],
    damping = 0.7,
    mixing = LdosMixing()
)

display(@time analyze_convergence(get_system(eq_result), dftk_potential, (5:2.5:25)u"hartree"))

;


text:
using OrgTables
using Base.Test

# r,h=readorg("testtable3.org", HeaderT=String)
# @test h==String[]
# @test typeof(r)==Matrix{String}

# r,h=readorg("testtable1.org")
# @test h==String["X","Y","id","r-dists contour","f","contours","dB","bed ele","bmax","Surf DEM","remarks"]
# @test typeof(r)==Matrix{String}
# @test_throws ErrorException readorg("testtable1.org", String, "", 2)
# @test_throws ErrorException readorg("testtable1.org", String, "", 3)

# r,h=readorg("testtable2.org", Float64, NaN)
# @test h==String["X","Y","r-dists contour","f","dB","bed ele","bmax","Surf DEM"]
# @test typeof(r)==Matrix{Float64}
# @test all([x[1]===x[2] for x in zip(r[1,:],[624479.0, 232874.0, 1.0, NaN, NaN, 376.0, NaN, 434.0])])

# r,h=readorg("testtable2.org", String, "", 2)
# @test h==String["X","Y","r-dists contour","f","dB","bed ele","bmax","Surf DEM"]
# @test typeof(r)==Matrix{String}
# @test all([x[1]==x[2] for x in zip(r[:,1],["a624479", "a619161"])])
# @test_throws ErrorException readorg("testtable2.org", String, "", 3)

# r,h=readorg("testtable1.org", Float64, NaN, drop=[3,6,11])
# @test h==String["X","Y","r-dists contour","f","dB","bed ele","bmax","Surf DEM"]
# @test typeof(r)==Matrix{Float64}

# test FileIO.jl inter-op
using FileIO
r,h=load("testtable3.org", HeaderT=String)
@test h==String[]
@test typeof(r)==Matrix{String}

r,h=load("testtable1.org")
@test h==String["X","Y","id","r-dists contour","f","contours","dB","bed ele","bmax","Surf DEM","remarks"]
@test typeof(r)==Matrix{String}
@test_throws ErrorException load("testtable1.org", tablenr=2)
@test_throws ErrorException load("testtable1.org", tablenr=3)

r,h=load("testtable2.org", T=Float64, fillval=NaN)
@test h==String["X","Y","r-dists contour","f","dB","bed ele","bmax","Surf DEM"]
@test typeof(r)==Matrix{Float64}
@test all([x[1]===x[2] for x in zip(r[1,:],[624479.0, 232874.0, 1.0, NaN, NaN, 376.0, NaN, 434.0])])
r,h=load("testtable2.org", tablenr=2)
@test h==String["X","Y","r-dists contour","f","dB","bed ele","bmax","Surf DEM"]
@test typeof(r)==Matrix{String}
@test all([x[1]==x[2] for x in zip(r[:,1],["a624479", "a619161"])])

r,h=load("testtable1.org", T=Float64, fillval=NaN, drop=[3,6,11])
@test h==String["X","Y","r-dists contour","f","dB","bed ele","bmax","Surf DEM"]
@test typeof(r)==Matrix{Float64}