text:
module SQLiteLoader

# https://github.com/JuliaDatabases/SQLite.jl
using SQLite # SQLite.jl 1.1.14
using .SQLite: Tables, DBInterface

using Octo: Repo, AdapterBase, DBMS, SQLElement, Structured
using .Repo: SQLKeyword, ExecuteResult, sql_startswith_insert_update_delete
using .AdapterBase: INSERT

# db_dbname
function db_dbname(nt::NamedTuple)::String
    (last ∘ splitdir)(nt.dbfile)
end

# db_connect
function db_connect(; kwargs...)
    dbfile = getindex(kwargs, :dbfile)
    SQLite.DB(dbfile)
end

# db_disconnect
function db_disconnect(db)
    DBInterface.close!(db)
end

# query
function get_query_result(q::SQLite.Query)
    _st = SQLite._stmt(q.stmt)
    cnt = SQLite.sqlite3_column_count(_st.handle)
    if iszero(cnt)
        NamedTuple[]
    else
        Tables.rowtable(q)
    end
end

function query(db, sql::String)
    q = DBInterface.execute(db, sql)
    get_query_result(q)
end

function query(db, prepared::String, vals::Vector)
    stmt = DBInterface.prepare(db, prepared)
    q = DBInterface.execute(stmt, vals)
    get_query_result(q)
end

# execute
function execute(db, sql::String)::ExecuteResult
    DBInterface.execute(db, sql)
    sql_startswith_insert_update_delete_then_get_num_affected_rows(sql, db)
end

function execute(db, prepared::String, vals::Vector)::ExecuteResult
    stmt = DBInterface.prepare(db, prepared)
    DBInterface.execute(stmt, vals)
    sql_startswith_insert_update_delete_then_get_num_affected_rows(prepared, db)
end

function execute(db, prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult
    stmt = DBInterface.prepare(db, prepared)
    num_affected_rows = 0
    for nt in nts
        DBInterface.execute(stmt, values(nt))
        num_affected_rows += get_num_affected_rows(db)
    end
    (num_affected_rows=num_affected_rows,)
end

function sql_startswith_insert_update_delete_then_get_num_affected_rows(sql::String, db)
    if sql_startswith_insert_update_delete(sql)
        num_affected_rows = get_num_affected_rows(db)
        (num_affected_rows=num_affected_rows,)
    else
        nothing
    end
end

function get_num_affected_rows(db)::Int
    SQLite.sqlite3_total_changes(db.handle)
end

# execute_result
function execute_result(db, command::SQLKeyword)::NamedTuple
    if INSERT === command
        last_insert_id = SQLite.last_insert_rowid(db)
        (id=last_insert_id,)
    else
        NamedTuple()
    end
end

function Base.show(io::IO, mime::MIME"text/plain", element::Union{E,Structured} where E<:SQLElement)
    dbms = DBMS.SQLite()
    AdapterBase._show(io, mime, dbms, element)
end

end # module Octo.Backends.SQLiteLoader


text:
export FunctionalTable, columns, ordering, rename

"""
$(TYPEDEF)

# Internal notes

- Use accessors `length`, `colums`, and `ordering` to access the fields, property accessors
  are forwarded to `columns`.

- The only inner constructor is the one where both the length and the ordering is trusted
  (and thus unchecked). Outer constructors should first wrap the ordering rule, then
  compute/verify length.
"""
struct FunctionalTable{C <: NamedTuple, O <: TableOrdering}
    len::Int
    columns::C
    ordering::O
    function FunctionalTable(trust_length::TrustLength, columns::C, ordering_rule::R
                             ) where {C <: NamedTuple, R <: TrustOrdering}
        @unpack ordering = ordering_rule
        @argcheck key_issubset(orderkey.(ordering), keys(columns))
        new{C, typeof(ordering)}(trust_length.len, columns, ordering)
    end
end

####
#### Outer constructors from FunctionalTable or NamedTuple
####

function FunctionalTable(ft::FunctionalTable,
                         ordering_rule::OrderingRule{K} = VerifyOrdering()) where K
    new_ordering = ordering_rule.ordering
    K ≡ :trust && return FunctionalTable(TrustLength(length(ft)), columns(ft), ordering_rule)
    rule = is_prefix(new_ordering, ordering(ft)) ? TrustOrdering(new_ordering) : ordering_rule
    FunctionalTable(TrustLength(length(ft)), columns(ft), rule)
end

function FunctionalTable(len::TrustLength, columns::NamedTuple,
                         ordering_rule::VerifyOrdering = VerifyOrdering(()))
    ft = FunctionalTable(len, columns, TrustOrdering(ordering_rule))
    @argcheck issorted(ft; lt = (a, b) -> isless_ordering(ordering(ft), a, b))
    ft
end

function FunctionalTable(len::TrustLength, columns::NamedTuple, ::TryOrdering)
    error("not implemented yet, maybe open an issue?")
end

function FunctionalTable(len::Integer, columns::NamedTuple, ordering_rule::OrderingRule)
    @argcheck all(column -> length(column) == len, values(columns))
    FunctionalTable(TrustLength(len), columns, ordering_rule)
end

function FunctionalTable(columns::NamedTuple, ordering_rule::OrderingRule)
    @argcheck !isempty(columns) "At least one column is needed to determine length."
    len = length(first(columns))
    @argcheck all(column -> length(column) == len, Base.tail(values(columns)))
    FunctionalTable(TrustLength(len), columns, ordering_rule)
end

FunctionalTable(len::Integer) =
    FunctionalTable(TrustLength(len), NamedTuple(), TrustOrdering())

####
#### accessors for fields and property overloading
####

"""
$(SIGNATURES)

Return the columns in a `NamedTuple`.

Each column is an iterable, but not necessarily an `<: AbstractVector`.

!!! note
    **Never mutate columns obtained by this method**, as that will violate invariants
    assumed by the implementation. Use `map(collect, columns(ft))` or similar to obtain
    mutable vectors.
"""
@inline columns(ft::FunctionalTable) = getfield(ft, :columns)

"""
$(SIGNATURES)

Return the ordering of the table, which is a tuple of `ColumnOrdering` objects.
"""
ordering(ft::FunctionalTable) = getfield(ft, :ordering)

@inline Base.propertynames(ft::FunctionalTable) = propertynames(columns(ft))

@inline Base.getproperty(ft::FunctionalTable, key::Symbol) = getproperty(columns(ft), key)

@inline Base.keys(ft::FunctionalTable) = keys(columns(ft))

@inline Base.pairs(ft::FunctionalTable) = pairs(columns(ft))

@inline Base.values(ft::FunctionalTable) = values(columns(ft))

####
#### Iteration interface and constructor
####

Base.IteratorSize(::FunctionalTable) = Base.HasLength()

@inline Base.length(ft::FunctionalTable) = getfield(ft, :len)

Base.IteratorEltype(::FunctionalTable) = Base.HasEltype()

@inline Base.eltype(::Type{<: FunctionalTable{C}}) where C =
    NamedTuple{fieldnames(C), Tuple{map(eltype, fieldtypes(C))...}}

"""
$(SIGNATURES)

Create a `FunctionalTable` from an iterable that returns `NamedTuple`s.

Returned values need to have the same names (but not necessarily types).

`ordering_rule` specifies sorting. The `VerifyOrdering` (default), `TrustOrdering`, and
`TryOrdering`  constructors take a tuple of a tuple of `:key` or `:key => reverse` elements.

`cfg` determines sink configuration for collecting elements of the columns, see
[`SinkConfig`](@ref).
"""
function FunctionalTable(itr, ordering_rule::OrderingRule = TrustOrdering();
                         cfg::SinkConfig = SINKCONFIG)
    FunctionalTable(collect_columns(cfg, itr, ordering_rule)...)
end

function Base.iterate(ft::FunctionalTable, states...)
    ys = map(iterate, columns(ft), states...)
    any(isequal(nothing), ys) && return nothing
    map(first, ys), map(last, ys)
end

"""
Shows this many values from each column in a `FunctionalTable`.
"""
const SHOWROWS = 5

function _showcolcontents(io::IO, itr)
    elts = collect(Iterators.take(itr, SHOWROWS + 1))
    print(io, eltype(itr), "[")
    for (i, elt) in enumerate(elts)
        i > 1 && print(io, ", ")
        i > SHOWROWS ? print(io, "…") : show(io, elt)
    end
    print(io, "]")
end

function Base.show(io::IO, ft::FunctionalTable)
    print(io, "FunctionalTable of $(length(ft)) rows, ")
    printstyled(io, ordering_repr(ordering(ft)); color = :green)
    ioc = IOContext(io, :compact => true)
    for (key, col) in pairs(columns(ft))
        println(ioc)
        print(ioc, "    ")
        printstyled(ioc, key; color = :blue)
        print(ioc, " = ")
        _showcolcontents(ioc, col)
    end
end

"""
$(SIGNATURES)

With a tuple of symbols returns `FunctionalTable` with a subset of the columns.

With a single symbol, return that column (an iterable).

`[drop = spec]` will keep *all but* the given columns, where `spec` is a `Tuple` of
`Symbol`s.

# Example

```julia
ft[(:a, :b)]
ft[:a]
ft[drop = (:a, :b)]
 ```
"""
function Base.getindex(ft::FunctionalTable, keep::Keys)
    @argcheck key_issubset(keep, keys(ft)) "Not all keys in $(keep) are in $(keys(ft))."
    FunctionalTable(TrustLength(length(ft)), NamedTuple{keep}(columns(ft)),
                    TrustOrdering(mask_ordering(ordering(ft), keep)))
end

Base.getindex(ft::FunctionalTable, key::Symbol) = columns(ft)[key]

function Base.getindex(ft::FunctionalTable; drop::Keys)
    @argcheck key_issubset(drop, keys(ft)) "Not all keys in $(drop) are in $(keys(ft))."
    keep = key_setdiff(keys(ft), drop)
    FunctionalTable(TrustLength(length(ft)), NamedTuple{keep}(columns(ft)),
                    TrustOrdering(mask_ordering(ordering(ft), drop, true)))
end

"""
$(SIGNATURES)

Rename the columns of a `FunctionalTable`. The second argument should be a `NamedTuple` of
`src = dest` pairs, where `dest` is a symbol.

# Example

```julia
rename(ft, (α = :a, β = :b)) # rename `α` to `a` and `β` to `b`
```
"""
function rename(ft::FunctionalTable, changes::NamedTuple{srckeys, <: Keys}) where srckeys
    destkeys = values(changes)
    @argcheck allunique(destkeys) "Non-unique destination keys."
    @argcheck srckeys ⊆ keys(columns(ft)) "Source keys not in table."
    change(key) = (key ∈ srckeys ? changes[key] : key)
    newkeys = map(change, keys(columns(ft)))
    newordering = map(o -> ColumnOrdering{change(orderkey(o)), orderrev(o)}(), ordering(ft))
    FunctionalTable(TrustLength(length(ft)),
                    NamedTuple{newkeys}(values(columns(ft))),
                    TrustOrdering(newordering))
end

"""
$(SIGNATURES)

Merge two `FunctionalTable`s.

When `replace == true`, columns in the first one are replaced by second one, otherwise an
error is thrown if column names overlap.

The second table can be specified as a `NamedTuple` of columns.
"""
function Base.merge(a::FunctionalTable, b::FunctionalTable; replace = false)
    @argcheck length(a) == length(b)
    if !replace
        dup = tuple((keys(a) ∩ keys(b))...)
        @argcheck isempty(dup) "Duplicate columns $(dup). Use `replace = true`."
    end
    FunctionalTable(TrustLength(length(a)), merge(columns(a), columns(b)),
                    TrustOrdering(mask_ordering(ordering(a), keys(b), true)))
end

"""
$(SIGNATURES)
"""
Base.map(f, ft::FunctionalTable; cfg = SINKCONFIG) = FunctionalTable(imap(f, ft); cfg = cfg)

"""
$(SIGNATURES)

Map `ft` using `f` by rows, then `merge` the two. See [`map(f, ::FunctionalTable)`](@ref).

`cfg` is passed to `map`, `replace` governs replacement of overlapping columns in `merge`.
"""
function Base.merge(f, ft::FunctionalTable; cfg = SINKCONFIG, replace = false)
    merge(ft, map(f, ft; cfg = cfg); replace = replace)
end

Base.merge(ft::FunctionalTable, newcolumns::NamedTuple; replace = false) =
    merge(ft, FunctionalTable(newcolumns); replace = replace)

"""
$(SIGNATURES)

Filter a FunctionalTable using a prediate that maps rows (as `NamedTuple`s) to `Bool`.
Preserves type (and thus ordering) of the table.
"""
function Base.filter(f, ft::FunctionalTable; cfg = SINKCONFIG)
    FunctionalTable(collect_columns(cfg, Iterators.filter(f, ft),
                                    TrustOrdering(ordering(ft)), eltype(ft))...)
end

"""
$(SIGNATURES)

A `FunctionalTable` of the first `n` rows.

Useful for previews and data exploration. Preserves type (and thus ordering) of the table.
"""
function Base.first(ft::FunctionalTable, n::Integer; cfg::SinkConfig = SINKVECTORS)
    i = 0
    filter(_ -> (i += 1; i ≤ n), ft; cfg = cfg)
end

ext: jl


text:
# Julia1_3.6.jl Modules
# defining a Module

module Gravity      # introduces new global scope
#using somelib       # presuming somelib contains something you need

#import Base.show # if you just want to use show you could use using

# import allows you to extend show with new methods
g=9.80665           # m/s^2
v(t) = g*t            # velocity as a function of time
s(t) = 0.5*g*t^2   # Distance as a function of time
export v, s, g      # make these functions and variables global
end

using Main.Gravity
for t=0.0:0.5:5.0
    print("t=$t v=$(v(t))\n      s=$(s(t)) \n")
end
#
###
#
module Astronomy
# constants
const c = 299792.458     # speed of light km/s
const dsec=24.0*3600.0   # seconds in a day
const ly = 365.25*dsec*c # light year in km
const au = 149597870.7   # astronomical unit in km (Radius of earths orbit)
const parsec = au*648000.0/π   # parsec in km (Distance to create parallax
#                              of 1 arcsecond of parallax from earths orbit)
export c
# functions
parsectoly(ps)= ps*parsec/ly  #  cnv parsecs to light years
stardis(parallax)=parsectoly(1.0/parallax)  # calc star dist
# An angle of 1 arcsecond of a degree  1/3600 of a degree
# from an object
export stardis
end

using Printf
using Main.Astronomy
nearstars = Dict("ProximaCentauri"=>.76850, "AlphaCentauri"=>.74723,
 "Bernard's Star"=>.54745, "Luhman16"=>.501398,
 "Wise0855-714"=>.449, "Wolf"=>.41516,
 "Lalande21185"=>.30264, "Sirius A"=>.37668)
@printf("       Star name   parallax     d.lightyears\n")
for star in nearstars
    name,plax=star
    @printf("%16s %10.6f %10.4f\n",name,plax,stardis(plax))
end

println(ly)  # error since not exported

println(Astronomy.ly) # But can still get ly using Module name
println(Astronomy.parsectoly(1.0))  # same with functions


text:
# Testing a pretrained normalizing flow for 2D Rosenbrock example
# Authors: Ali Siahkoohi, alisk@gatech.edu
#          Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
#          Philipp Witte, pwitte3@gatech.edu
# Date: September 2020
"""
Example:
julia train_unsupervised_hint_rosenbrock.jl
"""

using DrWatson
@quickactivate :FastApproximateInference

using InvertibleNetworks
using LinearAlgebra
using Random
using Statistics
using ArgParse
using BSON
using PyPlot
using Seaborn
using Printf
using Flux
using Flux.Optimise: update!

set_style("whitegrid")
rc("font", family="serif", size=13)
font_prop = matplotlib.font_manager.FontProperties(
    family="serif",
    style="normal",
    size=13
)
sfmt=matplotlib.ticker.ScalarFormatter(useMathText=true)
sfmt.set_powerlimits((0, 0))

# Random seed
Random.seed!(19)

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 25
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-3
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 1
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 64
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 64
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 8
    "--sigma"
        help = "noise std"
        arg_type = Float32
        default = 4f-1
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "unsupervised-hint-rosenbrock"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
sigma = parsed_args["sigma"]
sim_name = parsed_args["sim_name"]

# Define network
nx = 1
ny = 1
n_in = 2

CH_rev_precond = NetworkConditionalHINT(n_in, n_hidden, depth)
CH_rev = NetworkConditionalHINT(n_in, n_hidden, depth)
CH = NetworkConditionalHINT(n_in, n_hidden, depth)

# Loading the experiment—only network weights and training loss
Params, fval, exp_path = load_experiment(parsed_args; return_path=true)

Y_fixed = wload(exp_path)["Y_obs"]
Zy_fixed = wload(exp_path)["Zy_obs"]
X_fixed = wload(exp_path)["X_true"]
Zx = wload(exp_path)["Z_train"]
A = wload(exp_path)["A.A"]
warmstarted_args = wload(exp_path)["warmstarted_args"]

put_params!(CH_rev, Params)
CH_rev = reverse(CH_rev)

Params_precond, fval_precond, exp_path_precond = load_experiment(
    warmstarted_args;
    return_path=true
)
pretrained_args = wload(exp_path_precond)["pretrained_args"]

put_params!(CH_rev_precond, Params_precond)
CH_rev_precond = reverse(CH_rev_precond)

Params0, fval0 = load_experiment(pretrained_args)
put_params!(CH, Params0)

# Forward operator
A = SquareCSop(A, n_in)

# Testing data
test_size = size(Zx, 4)
RB_dist = RosenbrockDistribution(0f0, 5f-1, ones(Float32, 2, 1))
X_test = rand(RB_dist, 20*test_size)
Y_test = A*X_test + sigma*randn(Float32, nx, ny, n_in, 20*test_size)
Y_train = X_test + sigma*randn(Float32, nx, ny, n_in, 20*test_size)

# Drawing posterior samples given y* (from trained network w/ warmstart)
X_post_precond = CH_rev_precond.forward(
    randn(Float32, nx, ny, n_in, 20*test_size),
    repeat(Zy_fixed, 1, 1, 1, 20*test_size)
)[1]

X_post = CH_rev.forward(
    randn(Float32, nx, ny, n_in, 20*test_size),
    repeat(Zy_fixed, 1, 1, 1, 20*test_size)
)[1]

X_post0 = CH.inverse(
    randn(Float32, nx, ny, n_in, 20*test_size),
    repeat(Zy_fixed, 1, 1, 1, 20*test_size)
)[1]

# Create Flux parameter
x̂ = randn(Float32, size(X_fixed))
θ = Flux.Params([x̂])

# Optimizer
η = 0.08
opt = pSGLD(η)
max_itr = 40*test_size
X_sgld = zeros(Float32, 1, 1, 2, max_itr)

function objective(x_in, y_in, RB_dist)
    return (1f0/(2f0*sigma^2f0))*sum((A*x_in - y_in).^2f0) + sum(nl_pdf(x_in, RB_dist))
end

for j = 1:max_itr

    # Evaluate objective and gradients
    obj, back = Flux.pullback(θ) do
      objective(x̂, Y_fixed, RB_dist)
    end
    grads = back(1.0f0)
    update!(opt, x̂, grads[x̂])

    if j%100 == 0
      @printf("SGLD iterations: [%d/%d] | Objective: %4.2f\r", j, max_itr, obj)
    end
    X_sgld[:, :, :, j:j] = x̂
end

# Warm-up phase
X_sgld = X_sgld[:, :, :, fld(max_itr, 2):end]

# Thinning
X_sgld = X_sgld[:, :, :, 1:5:end]

save_dict = @strdict max_epoch lr lr_step batchsize n_hidden depth
save_path = plotsdir(sim_name, savename(save_dict))

Training loss
fig = figure("training logs", figsize=(5, 5))
plot(
    range(1, max_epoch, length=length(fval)), fval,
    label="Equation (2)", color="#addbc2"
)
plot(
    range(1, warmstarted_args["max_epoch"], length=length(fval_precond)), fval_precond,
    label="Equation (6)", color="#ff5900"
)
plot(
    range(0, max_epoch+1, length=length(fval)), 4.3f0*ones(Float32, length(fval)),
    color="k", linestyle="--", linewidth=.5
)
title("Training objective")
ylabel(L"$\mathbb{D}_{KL} + const.$")
xlabel("Epochs")
xlim([.6, max_epoch+1])
legend(loc="upper right", ncol=1)
wsave(joinpath(save_path, "training-obj.png"), fig)
close(fig)

# True samples from Rosenbrock distribuition
fig = figure("rosenbrock samples", figsize=(5, 5))
cmap = matplotlib.cm.Greys(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[13:end,1:end])
ax = kdeplot(
    X_test[1, 1, 1, :], X_test[1, 1, 2, :], shade=false, cut=0, levels=15,
    cmap=cmap, bw_adjust=1.5, alpha=1.0
)
scatter(
    X_fixed[1, 1, 1, :], X_fixed[1, 1, 2, :], s=30.0, color="#000000", marker="*",
    label="true model"
)
ax.set_xlim([-3, 3])
ax.set_ylim([-2.5, 7])
ax.set_ylabel(L"$x_2$")
ax.set_xlabel(L"$x_1$")
ax.legend(loc="upper right", ncol=2)
ax.set_title("Rosenbrock prior density")
wsave(joinpath(save_path, "true-samples.png"), fig)
close(fig)

# Data samples
cmap = matplotlib.cm.Greys_r(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[13:end,1:end])
fig = figure("data samples", figsize=(5, 5))
ax = kdeplot(
    Y_test[1, 1, 1, :], Y_test[1, 1, 2, :], shade=false, cut=0, levels=15,
    bw_adjust=1.7, alpha=1.0, cmap=cmap
)
label1 = matplotlib.patches.Patch(color="#4a4a4a", label=L"$\mathbf{y} = \mathbf{A} \mathbf{x} + \mathbf{\epsilon}$")
cmap = matplotlib.cm.Reds(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[20:end,1:end])
fig = figure("data samples", figsize=(5, 5))
ax = kdeplot(
    Y_train[1, 1, 1, :], Y_train[1, 1, 2, :], shade=false, cut=0, levels=15,
    bw_adjust=1.7, cmap=cmap, alpha=0.5
)
label2 = matplotlib.patches.Patch(color="#e05a5a", label=L"$\mathbf{y} = \mathbf{x} + \mathbf{\epsilon}$")
ax.legend(
    handles=[label1, label2], loc="upper left", ncol=2,
    markerscale=5
)
ax.set_xlim([-3, 3])
ax.set_ylim([-2.5, 7])
ax.set_ylabel(L"$y_2$")
ax.set_xlabel(L"$y_1$")
ax.set_title("Data space")
wsave(joinpath(save_path, "test-data-samples.png"), fig)
close(fig)

# Predicted posterior samples
fig = figure("rosenbrock samples", figsize=(5, 5))
cmap = matplotlib.cm.BuGn(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[15:end,1:end])
ax = kdeplot(
    X_post[1, 1, 1, :], X_post[1, 1, 2, :], shade=true, cut=0, levels=10,
    cmap=cmap, alpha=0.4, bw_adjust=1.7
)
label1 = matplotlib.patches.Patch(color="#addbc2", label="Equation (2)")

cmap = matplotlib.cm.Reds(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[13:end,1:end])
kdeplot(
    X_post_precond[1, 1, 1, :], X_post_precond[1, 1, 2, :], cmap=cmap, shade=false,
    alpha=1.0, cut=0, levels=10, linewidth=1, bw_adjust=1.7
)
label2 = matplotlib.patches.Patch(color="#ff5900", label="Equation (6)")

label3 = scatter(
    X_sgld[1, 1, 1, :], X_sgld[1, 1, 2, :], s=.5, color="#000000", alpha=.3,
    label="MCMC"
)
ax.set_xlim([0.75, 2.5])
ax.set_ylim([1.75, 5.0])
ax.set_ylabel(L"$x_2$")
ax.set_xlabel(L"$x_1$")
ax.set_title("Posterior density")
ax.legend(
    handles=[label1, label2, label3], loc="upper left", ncol=2,
    markerscale=5
)
wsave(joinpath(save_path, "unsupervised-posterior-samples.png"), fig)
close(fig)

fig = figure("rosenbrock samples", figsize=(5, 5))
cmap = matplotlib.cm.Greys(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[13:end,1:end])
ax = kdeplot(
    X_test[1, 1, 1, :], X_test[1, 1, 2, :], shade=false, cut=0, levels=15,
    cmap=cmap, bw_adjust=1.5, alpha=0.9
)
label1 = matplotlib.patches.Patch(color="#4a4a4a", label="Prior")

cmap = matplotlib.cm.Reds(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[13:end,1:end])
kdeplot(
    X_post_precond[1, 1, 1, :], X_post_precond[1, 1, 2, :], cmap=cmap, shade=false,
    alpha=0.8, cut=0, levels=10, linewidth=1, bw_adjust=1.7
)
label2 = matplotlib.patches.Patch(color="#ff5900", label="Equation (6)")

cmap = matplotlib.cm.Blues(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[15:end,1:end])
ax = kdeplot(
    X_post0[1, 1, 1, :], X_post0[1, 1, 2, :], shade=true, cut=0, levels=10,
    cmap=cmap, alpha=0.4, bw_adjust=1.7
)
label3 = matplotlib.patches.Patch(color="#65aadb", label="Equation (3)")
ax.set_xlim([-3, 3])
ax.set_ylim([-2.5, 7])
ax.set_ylabel(L"$x_2$")
ax.set_xlabel(L"$x_1$")
ax.legend(
    handles=[label1, label2, label3], loc="upper left", ncol=2,
    markerscale=5
)
ax.set_title("Prior vs. posterior density")
wsave(joinpath(save_path, "warmstarted-overlaid-posterior.png"), fig)
close(fig)

# Preeicted prior samples and posterior samples
fig = figure("rosenbrock samples", figsize=(5, 5))
cmap = matplotlib.cm.Blues(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[15:end,1:end])
ax = kdeplot(
    X_post0[1, 1, 1, :], X_post0[1, 1, 2, :], shade=true, cut=0, levels=10,
    cmap=cmap, alpha=0.4, bw_adjust=1.7
)
label1 = matplotlib.patches.Patch(color="#65aadb", label="Equation (3)")

cmap = matplotlib.cm.OrRd(range(0, 1, length=50))
cmap = matplotlib.colors.ListedColormap(cmap[13:end,1:end])
kdeplot(
    X_post_precond[1, 1, 1, :], X_post_precond[1, 1, 2, :], cmap=cmap, shade=false, alpha=1.0, cut=0,
    levels=10, linewidth=1, bw_adjust=1.7
)
label2 = matplotlib.patches.Patch(color="#ff5900", label="Equation (6)")
ax.set_xlim([0.75, 2.5])
ax.set_ylim([1.75, 5.0])
ax.set_ylabel(L"$x_2$")
ax.set_xlabel(L"$x_1$")
ax.set_title("Low- and high-fidelity Posterior")
ax.legend(
    handles=[label1, label2], loc="upper left", ncol=2,
    markerscale=5
)
wsave(joinpath(save_path, "prior-posterior-samples.png"), fig)
close(fig)


text:
export monitorpowersupplies

"""
    enablePositivePowerSupply!(dev::DWFDevice)

Enable positive power supply. Remember to AnalogIOEnable!(dev).
"""
function enablePositivePowerSupply!(dev::DWFDevice)
    checkopen(dev)
    FDwfAnalogIOChannelNodeSet(dev.hdwf[], 0, 0, 3)
    return nothing
end

"""
    setPositivePowerSupply!(dev::DWFDevice, voltage::T) where {T<:AbstractFloat}

Set the positive power supply voltage.
"""
function setPositivePowerSupply!(dev::DWFDevice, voltage::T) where {T<:Real}
    checkopen(dev)
    (voltage < 0) && error("Positive voltage should be positive")
    FDwfAnalogIOChannelNodeSet(dev.hdwf[], 0, 1, Cdouble(voltage))
    return nothing
end

"""
    enableNegativePowerSupply!(dev::DWFDevice)

Enable negative power supply. Remember to AnalogIOEnable!(dev).
"""
function enableNegativePowerSupply!(dev::DWFDevice)
    checkopen(dev)
    FDwfAnalogIOChannelNodeSet(dev.hdwf[], 1, 0, 3)
    return nothing
end

"""
    setNegativePowerSupply!(dev::DWFDevice, voltage::T) where {T<:AbstractFloat}

Set the negative power supply voltage.
"""
function setNegativePowerSupply!(dev::DWFDevice, voltage::T) where {T<:Real}
    checkopen(dev)
    (voltage > 0) && error("Negative voltage should be negative")
    FDwfAnalogIOChannelNodeSet(dev.hdwf[], 1, 1, Cdouble(voltage))
    return nothing
end

"""
    AnalogIOEnable!(dev::DWFDevice; state::Bool=true)

Enable/disable analog IO. (Enable by default)
"""
function AnalogIOEnable!(dev::DWFDevice; state::Bool=true)
    checkopen(dev)
    FDwfAnalogIOEnableSet(dev.hdwf[], state)
    return nothing
end

"""
    AnalogIOReset!(dev::DWFDevice)

Disable then enable analog IO. Useful for resetting after overcurrent protection.
"""
function AnalogIOReset!(dev::DWFDevice)
    checkopen(dev)
    FDwfAnalogIOEnableSet(dev.hdwf[], false)
    FDwfAnalogIOEnableSet(dev.hdwf[], true)
    return nothing
end

"""
    AnalogIOEnabled(dev::DWFDevice)

Return whether the analogIO is enabled.
"""
function AnalogIOEnabled(dev::DWFDevice)
    checkopen(dev)
    fOn = Ref(Int32(0))
    FDwfAnalogIOEnableStatus(hdwf[], fOn)
    return fOn[] == 1
end

"""
    PowerSupplyMonitor(dev::DWFDevice)

Monitor the device power draw via both USB and AUX power supplies.
"""
function PowerSupplyMonitor(dev::DWFDevice)
    checkopen(dev)
    vUSB = Ref(Cdouble(0))
    aUSB = Ref(Cdouble(0))
    vAUX = Ref(Cdouble(0))
    aAUX = Ref(Cdouble(0))

    # fetch analogIO status from device
    FDwfAnalogIOStatus(dev.hdwf[])

    # supply monitor
    FDwfAnalogIOChannelNodeStatus(dev.hdwf[], 2, 0, vUSB)
    FDwfAnalogIOChannelNodeStatus(dev.hdwf[], 2, 1, aUSB)
    FDwfAnalogIOChannelNodeStatus(dev.hdwf[], 3, 0, vAUX)
    FDwfAnalogIOChannelNodeStatus(dev.hdwf[], 3, 1, aAUX)

    return ("USB" => ("voltage" => vUSB[], "current" => aUSB[]),
    "AUX" => ("voltage" => vAUX[], "current" => aAUX[]))

end

function TemperatureMonitor(dev::DWFDevice)
    checkopen(dev)
    degDevice = Ref(Cdouble(0))

    # fetch analogIO status from device
    FDwfAnalogIOStatus(dev.hdwf[])

    FDwfAnalogIOChannelNodeStatus(dev.hdwf[], 2, 2, degDevice); #added from systemmonitor example
    return degDevice[]
end


text:
@testset "utility.jl" begin
    @test SphericalGeometry.normalize(-270.0, -180.0, 180.0) ≈ 90.0 atol = 0.1
    @test SphericalGeometry.normalize(181.0, -180.0, 180.0) ≈ -179.0 atol = 0.1

    # Source: https://gist.github.com/missinglink/d0a085188a8eab2ca66db385bb7c023a
    # Latitude positive
    @test SphericalGeometry.normalize(Point(55.555, 22.222)).ϕ ≈ 55.555 atol = 0.00001
    @test SphericalGeometry.normalize(Point(55.555, 22.222)).λ ≈ 22.222 atol = 0.00001

    @test SphericalGeometry.normalize(Point(1.0, 0.0)).ϕ ≈ 1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(1.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(91.0, 0.0)).ϕ ≈ 89.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(91.0, 0.0)).λ ≈ 180.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(181.0, 0.0)).ϕ ≈ -1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(181.0, 0.0)).λ ≈ 180.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(271.0, 0.0)).ϕ ≈ -89.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(271.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(361.0, 0.0)).ϕ ≈ 1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(361.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(631.0, 0.0)).ϕ ≈ -89.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(631.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(721.0, 0.0)).ϕ ≈ 1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(721.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    # Latitude negative
    @test SphericalGeometry.normalize(Point(-1.0, 0.0)).ϕ ≈ -1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-1.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(-91.0, 0.0)).ϕ ≈ -89.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-91.0, 0.0)).λ ≈ 180.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(-181.0, 0.0)).ϕ ≈ 1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-181.0, 0.0)).λ ≈ 180.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(-271.0, 0.0)).ϕ ≈ 89.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-271.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(-361.0, 0.0)).ϕ ≈ -1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-361.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(-631.0, 0.0)).ϕ ≈ 89.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-631.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(-721.0, 0.0)).ϕ ≈ -1.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(-721.0, 0.0)).λ ≈ 0.0 atol = 0.00001

    # Longitude positive
    @test SphericalGeometry.normalize(Point(0.0, 1.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, 1.0)).λ ≈ 1.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, 181.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, 181.0)).λ ≈ -179.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, 271.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, 271.0)).λ ≈ -89.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, 361.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, 361.0)).λ ≈ 1.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, 631.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, 631.0)).λ ≈ -89.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, 721.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, 721.0)).λ ≈ 1.0 atol = 0.00001

    # Longitude negative
    @test SphericalGeometry.normalize(Point(0.0, -1.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, -1.0)).λ ≈ -1.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, -181.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, -181.0)).λ ≈ 179.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, -271.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, -271.0)).λ ≈ 89.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, -361.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, -361.0)).λ ≈ -1.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, -631.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, -631.0)).λ ≈ 89.0 atol = 0.00001

    @test SphericalGeometry.normalize(Point(0.0, -721.0)).ϕ ≈ 0.0 atol = 0.00001
    @test SphericalGeometry.normalize(Point(0.0, -721.0)).λ ≈ -1.0 atol = 0.00001
end


text:
# parameters
mass_bin_name = ARGS[1]
tslice = "t1"
path_wavelist = "src"
path_to_working_folder = "data"

######################################################
using DelimitedFiles
push!(LOAD_PATH,"src")
using SDMHelper
using FittingPWALikelihood
using amplitudes_compass
using PWAHelper
using Random
rng = MersenneTwister()
serial= round(Int,rand(Random.seed!(rng))*100000)
print(serial)

# read wavelist throw waves below threshol
wavelist = get_wavelist(joinpath(path_wavelist,"wavelist_formated.txt");
     path_to_thresholds=joinpath(path_wavelist,"thresholds_formated.txt"),
     M3pi=Meta.parse(mass_bin_name[1:4])/1000)

const Nwaves = size(wavelist,1)

# Blocks in wavelist
const noϵ =  [1] # flat wave
const posϵ = [i for (i,ϵ) in enumerate(wavelist[:,6]) if ϵ=="+"]
const negϵ = [i for (i,ϵ) in enumerate(wavelist[:,6]) if ϵ=="-"]

# Model description
const ModelBlocks = [noϵ, posϵ, negϵ, negϵ]

# integrals from MC data
BmatMC = read_cmatrix(
    joinpath(path_to_working_folder,"integrmat_$(mass_bin_name)_$(tslice)_mc.txt"));

# normalize B
Bscale = [BmatMC[i,i] for i in 1:size(BmatMC,1)];
for i=1:Nwaves, j=1:Nwaves
    BmatMC[i,j] *= 1/sqrt(Bscale[i]*Bscale[j])
end

# load precalculated data array
const PsiRD = read_precalc_basis(
    joinpath(path_to_working_folder,"functions_$(mass_bin_name)_$(tslice)_rd.bin"));
# normalize Psi
for i in 1:size(PsiRD,2)
    PsiRD[:,i] .*= 1.0/sqrt(Bscale[i])
end

# get functions to calculate LLH, derivative and hessian
LLH, GRAD, LLH_and_GRAD!, HES = createLLHandGRAD(PsiRD, BmatMC, ModelBlocks);

# generate random parameters on a surface of intensity constraint
const Npar = get_npars(ModelBlocks)
const pars0 = rand(Npar);
pars0 .*= get_parameter_ranges(BmatMC, ModelBlocks)
normalize_pars!(pars0, BmatMC, ModelBlocks)

# fit
# start vast algorithm which goes directly to the minimum
@time minpars = minimize(LLH, LLH_and_GRAD!;
    algorithm = :LD_LBFGS, verbose=1, starting_pars=pars0,
    llhtolerance = 1e-4)

# scale parameters back
parscale = abs.(extnd(sqrt.(Bscale), get_parameter_map(ModelBlocks, Nwaves)))
final_pars = minpars ./ parscale

println("Done! Saving...")
writedlm(joinpath(path_to_working_folder,"llhfit_$(mass_bin_name)_$(tslice)_$(serial).txt"), final_pars)


text:
"""
ahs.jl

Types and functions for additive hard-sphere reference systems.

References:

* S. B. Yuste et al: J. Chem. Phys., Vol. 108, No.9, 1 (1998), 3683-3693.
* K. Hoshino. and W. H. Young: J. Phys. F: Metal Phys., Vol. 10 (1980), 1365-1374.

"""

immutable AHS <: IndependentReferenceSystem
  approx::Symbol # approximation
  σ::Vector{Float64} # hard-sphere diameters
  ρ::Vector{Float64} # number densities
  α::Float64 # α parameter for RFA approximation
end # type

function AHS(approx::Symbol, σ::Vector{Float64}, ρ::Vector{Float64})
  optimizealpha(AHS(approx, σ, ρ, 0.))
end

@doc doc"""
AHS(; kwargs...)

Initializes an additive hard-sphere reference system.
You may use one of possible combinations of keyword arguments.

One-component case:

* two out of σ, ρ, and η (all scalar)

Multi-component case:

* σ (vector) and ρ (vector)
* ρ (scalar), σ (vector), and c (vector)
* η (scalar), σ (vector), and c (vector)

You may use sigma, rho, or eta instead of σ, ρ, or η, respectively.
""" ->
function AHS(; kwargs...)
  keys, vals = expandkwargs(kwargs)

  # The type of approximation
  if in(:approx, keys)
    ind = findfirst(keys, :approx)
    approx = vals[ind]

    @assert in(approx, [:PY, :RFA]) "unsuportted type of approximation for AHS"

    deleteat!(keys, ind)
    deleteat!(vals, ind)
  else
    approx = :RFA # default
  end

  @assert 2 ≤ length(keys) ≤ 3 "wrong number of arguments for AHS"

  if all(x -> !(typeof(x) <: Vector), vals) # one-component case

    if greekin(:σ, keys) && greekin(:ρ, keys)
      σ = vals[greekfind(keys, :σ)]
      ρ = vals[greekfind(keys, :ρ)]

    elseif greekin(:ρ, keys) && greekin(:η, keys)
      ρ = vals[greekfind(keys, :ρ)]
      η = vals[greekfind(keys, :η)]
      σ = (6/π * η / ρ)^(1/3)

    elseif greekin(:η, keys) && greekin(:σ, keys)
      η = vals[greekfind(keys, :η)]
      σ = vals[greekfind(keys, :σ)]
      ρ = 6/π * η / σ^3

    else
      error("invalid arguments for AHS")
    end # if

    return AHS(approx, [σ], [ρ])
  else # Multi-component case

    if length(kwargs) == 2
      @assert all(x-> typeof(x) <: Vector, vals) "invalid arguments for AHS"
      σ = vals[greekfind(keys, :σ)]
      ρ = vals[greekfind(keys, :ρ)]

    else
      @assert greekin(:σ, keys) && greekin(:c, keys) "invalid arguments for AHS"
      σ = vals[greekfind(keys, :σ)]
      c = vals[findfirst(keys, :c)]

      if greekin(:ρ, keys)
        ρ = c .* vals[greekfind(keys, :ρ)]

      elseif greekin(:η, keys)
        η = vals[greekfind(keys, :η)]
        α = σ / maximum(σ)
        ρ = (6/π) ./ σ.^3 .* (c .* α.^3) ./ sum(c .* α.^3) * η
      end
    end

    @assert length(σ) == length(ρ) "invalid arguments for AHS"

    return AHS(approx, σ, ρ)
  end
end

ncomp(ahs::AHS)::Int = length(ahs.σ)
composition(ahs::AHS)::Vector{Float64} = ahs.ρ / sum(ahs.ρ)
numberdensity(ahs::AHS)::Vector{Float64} = ahs.ρ
totalnumberdensity(ahs::AHS)::Float64 = sum(ahs.ρ)
packingfraction(ahs::AHS)::Vector{Float64} = (π/6) * ahs.ρ .* ahs.σ.^3
totalpackingfraction(ahs::AHS)::Float64 = sum(packingfraction(ahs))
temperature(ahs::AHS)::Float64 = InvalTemp

function hsdiameter(sys::AHS)::Array{Float64,2}
  N = ncomp(sys)

  ret = Array{Float64,2}(N,N)

  for i in 1:N, j in 1:N
    if i == j
      ret[i,j] = sys.σ[i]
    else
      ret[i,j] = (sys.σ[i] + sys.σ[j]) / 2
    end
  end

  return ret
end

emptyradius(ahs::AHS)::Array{Float64,2} = hsdiameter(ahs)

function pairpotential(ahs::AHS)::Array{Function,2}
  N = ncomp(ahs)
  σ = hsdiameter(ahs)

  ret = Array{Function,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    function u(r)::Float64
      if r < σ[i,j]
        Inf
      else
        0
      end
    end

    ret[i,j] = ret[j,i] = u
  end

  return ret
end

function zeta(ahs::AHS, n::Int)::Float64
  N = ncomp(ahs)
  sum(ahs.ρ[i] * ahs.σ[i]^n for i = 1:N)
end

@doc doc"""
pressure(ahs::AHS)::Float64

Returns pressure of an additive hard-sphere system multiplied by inverse
temperature based on BMCSL equation of state.

""" ->
function pressure(ahs::AHS)::Float64
  ρ = totalnumberdensity(ahs)
  η = totalpackingfraction(ahs)
  ζ₁ = zeta(ahs, 1)
  ζ₂ = zeta(ahs, 2)

  ρ/(1-η) + π/2*ζ₁*ζ₂/(1-η)^2 + π^2/36*ζ₂^3*(3-η)/(1-η)^3
end

@doc doc"""
chemicalpotential(ahs::AHS)::Vector{Float64}

Returns a vector of excess (residual) chemical potential of each specie
multiplied by inverse temperature.

""" ->
function chemicalpotential(ahs::AHS)::Vector{Float64}
  N = ncomp(ahs)
  c = composition(ahs)
  σ = ahs.σ
  ρ = totalnumberdensity(ahs)

  R = σ/2
  S = π * σ.^2
  V = π/6 * σ.^3

  v = ρ * sum(c .* V)
  s = ρ * sum(c .* S)
  q = ρ * sum(c .* R.^2)
  r = ρ * sum(c .* R)

  βΔμ = Vector{Float64}(N)
  for i in 1:N
    βΔμ[i] = -log(1-v) + (R[i]*s + S[i]*r + V[i]*ρ)/(1-v) +
             (R[i]^2*s^2 + 2S[i]*q*s + 6V[i]*r*s)/6/(1-v)^2 +
             V[i]*q*s^2*(2 - v/3)/9/(1-v)^3
  end

  return βΔμ
end

function chemicalpotentialPYC(ahs::AHS)::Vector{Float64}
  N = ncomp(ahs)
  σ = hsdiameter(ahs)
  ρ = numberdensity(ahs)
  η = π/6*ρ
  ξ = sum( η[i]*σ[i,i]^3 for i in 1:N )
  X = sum( η[i]*σ[i,i]^2 for i in 1:N )
  Y = sum( η[i]*σ[i,i]^1 for i in 1:N )

  βp = ( sum(ρ)*(1 + ξ + ξ^2) -
         18/π * sum( η[i]*η[j]*(σ[i,i] - σ[j,j])^2*(2σ[i,j] + σ[i,i]*σ[j,j]*X)
                     for i in 1:N, j in 1:N if i < j )
       ) / (1-ξ)^3

  βΔμ = Vector{Float64}(N)
  for i in 1:N
    βΔμ[i] = log(1-ξ) + π/6*βp*σ[i,i]^3 + σ[i,i]^2/(1-ξ)^3 *
             ( 3Y*(1-2ξ) + 9/2*X^2*(1-ξ) + 3ξ^2*Y ) +
             σ[i,i]/(1-ξ)^3*(3X - 6ξ*X + 3ξ^2*X)
  end

  return βΔμ
end

# Ref: A. Meyer et al. Chemical Physics 49 (1980) 147-152
function contactvaluePY(ahs::AHS)::Matrix{Float64}
  N = ncomp(ahs)
  if N > 1
    error("calculation of contact values of multi-component system in PY approxmation is not implemented")
  end
  η = totalpackingfraction(ahs)
  ξ₀ = (1 + η/2) / (1 - η)^2
  g₀ = Matrix{Float64}(N,N)
  g₀[1,1] = ξ₀
  g₀
end

function contactvalueRFA(ahs::AHS)::Matrix{Float64}
  N::Int = ncomp(ahs)
  σ::Array{Float64,2} = hsdiameter(ahs)
  ρ::Vector{Float64} = numberdensity(ahs)

  # scalar constans
  ζ₂::Float64 = sum(ρ[i] * σ[i,i]^2 for i = 1:N)
  ζ₃::Float64 = sum(ρ[i] * σ[i,i]^3 for i = 1:N)
  η::Float64 = (π/6) * ζ₃
  λ::Float64 = 2π / (1-η)
  λ′::Float64 = π^2 * ζ₂ / (1-η)^2

  [ 1/2π * (λ + λ′/2 * σ[i,i]*σ[j,j]/σ[i,j] +
            λ′^2 / 18λ * σ[i,i]^2 * σ[j,j]^2 / σ[i,j]^2)
    for i = 1:N, j = 1:N ]
end

function contactvalue(ahs::AHS)::Matrix{Float64}
  if ahs.approx == :PY
    contactvaluePY(ahs)
  elseif ahs.approx == :RFA
    contactvalueRFA(ahs)
  end
end

function contactgradientPY(ahs::AHS)::Matrix{Float64}
  N = ncomp(ahs)
  if N > 1
    error("calculation of contact gradients of multi-component system in PY approxmation is not implemented")
  end
  η = totalpackingfraction(ahs)
  ξ₀ = (1 + η/2) / (1 - η)^2
  ξ₁ = 9/2 * η * (1 + η) / (1 - η)^3
  Y = Matrix{Float64}(N,N)
  Y[1,1] = ξ₀ / ξ₁
  Y
end

function contactgradientRFA(ahs::AHS)::Matrix{Float64}
  N::Int = ncomp(ahs)
  σ::Array{Float64,2} = hsdiameter(ahs)
  ρ::Vector{Float64} = numberdensity(ahs)
  α::Float64 = ahs.α
  g_σ::Array{Float64,2} = contactvalue(ahs)

  ζ₂::Float64 = sum(ρ[i] * σ[i,i]^2 for i = 1:N)
  ζ₃::Float64 = sum(ρ[i] * σ[i,i]^3 for i = 1:N)
  η::Float64 = (π/6) * ζ₃
  λ::Float64 = 2π / (1-η)
  λ′::Float64 = π^2 * ζ₂ / (1-η)^2

  L²::Array{Float64,2} =
    [ 2π * α * σ[i,j] * g_σ[i,j] for i = 1:N, j = 1:N ]

  L¹::Array{Float64,2} =
    [ ( λ*σ[i,j] + λ′/2*σ[i,i]*σ[j,j] + (λ + λ′*σ[i,i])*α -
        λ/2*σ[i,i]*sum(ρ[k]*σ[k,k]*L²[k,j] for k = 1:N) )
      for i = 1:N, j = 1:N ]

  g′::Array{Float64,2} =
    [ 1 / (2π*α*σ[i,j]) * (L¹[i,j] - L²[i,j]*(1/α + 1/σ[i,j]))
      for i in 1:N, j in 1:N ]
end

function contactgradient(ahs::AHS)::Matrix{Float64}
  if ahs.approx == :PY
    contactgradientPY(ahs)
  elseif ahs.approx == :RFA
    contactgradientRFA(ahs)
  end
end

#
# Pair correlation function in Laplace space
# ref: S. B. Yuste et al: J. Chem. Phys., Vol. 108, No.9, 1 (1998), 3683-3693.
#
function optimizealpha(ahs::AHS)
  ahs.approx == :PY && return ahs

  N::Int = ncomp(ahs)
  x::Vector{Float64} = composition(ahs)
  σ::Array{Float64,2} = hsdiameter(ahs)
  ρ::Vector{Float64} = numberdensity(ahs)

  # scalar constans
  ζ₁::Float64 = sum(ρ[i] * σ[i,i] for i = 1:N)
  ζ₂::Float64 = sum(ρ[i] * σ[i,i]^2 for i = 1:N)
  ζ₃::Float64 = sum(ρ[i] * σ[i,i]^3 for i = 1:N)
  η::Float64 = (π/6) * ζ₃
  λ::Float64 = 2π / (1-η)
  λ′::Float64 = π^2 * ζ₂ / (1-η)^2

  # auxiliary qunatities and functions
  I::Array{Float64,2} = δ::Array{Float64,2} = eye(N)
  ϕ₀(x::Complex{Float64}) = x^-1 * (1 - exp(-x))
  ϕ₁(x::Complex{Float64}) = x^-2 * (1 - x - exp(-x))
  ϕ₂(x::Complex{Float64}) = x^-3 * (1 - x + x^2/2 - exp(-x))

  # isothermal susceptibity
  χ::Float64 = sum(ρ) / (sum(ρ)/(1-η)^2 + π*ζ₁*ζ₂/(1-η)^3 +
                         π^2/36 * ζ₂^3 * (9-4η+η^2)/(1-η)^4)

  # the contact value of pair correlation function
  g_σ::Array{Float64,2} = contactvalue(ahs)

  L⁰ = Array{Float64,2}(N,N)
  L¹ = Array{Float64,2}(N,N)
  L² = Array{Float64,2}(N,N)
  A⁰ = Array{Float64,2}(N,N)
  A¹ = Array{Float64,2}(N,N)
  A² = Array{Float64,2}(N,N)
  A³ = Array{Float64,2}(N,N)
  B⁰ = Array{Float64,2}(N,N)
  B¹ = Array{Float64,2}(N,N)
  H⁰ = Array{Float64,2}(N,N)
  H¹ = Array{Float64,2}(N,N)

  function fopt(v::Vector{Float64}, g::Vector)::Float64
    α = v[1]

    # express L² in terms of α
    L² = [ 2π * α * σ[i,j] * g_σ[i,j] for i = 1:N, j = 1:N ]

    #
    # express H¹ in terms of L² and α
    #
    for i in 1:N, j in 1:N
      L⁰[i,j] = λ + λ′*σ[j,j] + 2λ′*α - λ*sum(ρ[k]*σ[k,k]*L²[k,j] for k = 1:N)
      L¹[i,j] = λ*σ[i,j] + λ′/2*σ[i,i]*σ[j,j] + (λ + λ′*σ[i,i])*α -
                λ/2*σ[i,i]*sum(ρ[k]*σ[k,k]*L²[k,j] for k = 1:N)

      Aᵢⱼ(n)::Float64 = (-1)^n * ρ[i] *
        ( σ[i,i]^(n+3) / factorial(n+3) * L⁰[i,j] -
          σ[i,i]^(n+2) / factorial(n+2) * L¹[i,j] +
          σ[i,i]^(n+1) / factorial(n+1) * L²[i,j] )

      A⁰[i,j] = Aᵢⱼ(0)
      A¹[i,j] = Aᵢⱼ(1)
      A²[i,j] = Aᵢⱼ(2)
      A³[i,j] = Aᵢⱼ(3)
    end

    for i in 1:N, j in 1:N
      B⁰[i,j] = 1/2π * L²[i,j] + sum(A²[k,j] for k = 1:N) -
        sum(σ[i,k]*(α*δ[k,j] - A¹[k,j]) for k = 1:N) -
        sum(1/2 * σ[i,k]^2 * (δ[k,j] - A⁰[k,j]) for k = 1:N)
    end

    H⁰ = B⁰ * inv(I - A⁰)

    for i in 1:N, j in 1:N
      B¹[i,j] = sum(A³[k,j] for k = 1:N) + sum(σ[i,k]*A²[k,j] for k = 1:N) -
        sum((1/2*σ[i,k]^2 + H⁰[i,k])*(α*δ[k,j] - A¹[k,j]) for k = 1:N) -
        sum((1/6*σ[i,k]^3 + σ[i,k]*H⁰[i,k])*(δ[k,j] - A⁰[k,j]) for k = 1:N)
    end

    H¹ = B¹ * inv(I - A⁰)

    #
    # calculate inversed isothermal susceptibity for given α
    #
    h0::Array{Float64,2} =
      [ -4π * √(ρ[i]*ρ[j]) * H¹[i,j] for i = 1:N, j = 1:N ]

    χ⁻¹ = sum(√(x[i]*x[j]) * inv(I + h0)[i,j] for i = 1:N, j = 1:N)

    abs(χ - 1/χ⁻¹)
  end # fopt

  # heuristic initial and boundary condition
  σ_m::Float64 = sum(x .* σ)
  α_init::Float64 = 0.01 * σ_m / η
  α_min::Float64 = 0.0 * α_init
  α_max::Float64 = 3.0 * α_init

  # opt = Optim.optimize(fopt, α_min, α_max, abs_tol=1e-6)
  # show(opt)
  # α::Float64 = Optim.minimizer(opt)

  opt = NLopt.Opt(:LN_BOBYQA, 1)
  NLopt.min_objective!(opt, fopt)
  NLopt.xtol_rel!(opt, 1e-4)
  NLopt.lower_bounds!(opt, [α_min])
  NLopt.upper_bounds!(opt, [α_max])

  (fmin, xmin, res) = NLopt.optimize(opt, [α_init])
  α::Float64 = xmin[1]

  if isapprox(α, α_min, rtol=1e-2) || isapprox(α, α_max, rtol=1e-2)
    warn("α = $α is on a bound of optimization (α_init = $α_init)")
  end

  AHS(ahs.approx, ahs.σ, ahs.ρ, α)
end # optimizealpha()

#
# Pair correlation function in Laplace space
# ref: S. B. Yuste et al: J. Chem. Phys., Vol. 108, No.9, 1 (1998), 3683-3693.
#
function paircorrelationlaplace(ahs::AHS)::Array{Function,2}
  N::Int = ncomp(ahs)
  x::Vector{Float64} = composition(ahs)
  ρ::Vector{Float64} = numberdensity(ahs)
  σ::Array{Float64,2} = hsdiameter(ahs)
  g_σ::Array{Float64,2} = contactvalue(ahs)

  # scalar constans
  ζ₂::Float64 = sum(ρ[i] * σ[i,i]^2 for i = 1:N)
  ζ₃::Float64 = sum(ρ[i] * σ[i,i]^3 for i = 1:N)
  η::Float64 = (π/6) * ζ₃
  λ::Float64 = 2π / (1-η)
  λ′::Float64 = π^2 * ζ₂ / (1-η)^2

  # auxiliary qunatities and functions
  I::Array{Float64,2} = δ::Array{Float64,2} = eye(N)
  ϕ₀(x::Complex{Float64}) = x^-1 * (1 - exp(-x))
  ϕ₁(x::Complex{Float64}) = x^-2 * (1 - x - exp(-x))
  ϕ₂(x::Complex{Float64}) = x^-3 * (1 - x + x^2/2 - exp(-x))

  α::Float64 = ahs.α

  L²::Array{Float64,2} =
    [ 2π * α * σ[i,j] * g_σ[i,j] for i = 1:N, j = 1:N ]

  L⁰::Array{Float64,2} =
    [ ( λ + λ′*σ[j,j] + 2λ′*α - λ*sum(ρ[k]*σ[k,k]*L²[k,j] for k = 1:N) )
      for i = 1:N, j = 1:N ]

  L¹::Array{Float64,2} =
    [ ( λ*σ[i,j] + λ′/2*σ[i,i]*σ[j,j] + (λ + λ′*σ[i,i])*α -
        λ/2*σ[i,i]*sum(ρ[k]*σ[k,k]*L²[k,j] for k = 1:N) )
      for i = 1:N, j = 1:N ]

  𝐋::Array{Function,2} =
    [ s::Complex{Float64} -> L⁰[i,j] + L¹[i,j]*s + L²[i,j]*s^2
      for i = 1:N, j = 1:N ]

  𝐀::Array{Function,2} =
    [ s::Complex{Float64} -> ρ[i] * ( ϕ₂(σ[i,i]*s)*σ[i,i]^3*L⁰[i,j] +
                                      ϕ₁(σ[i,i]*s)*σ[i,i]^2*L¹[i,j] +
                                      ϕ₀(σ[i,i]*s)*σ[i,i]^1*L²[i,j] )
      for i = 1:N, j = 1:N ]

  ret = Array{Function,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    function G(s::Complex{Float64})::Complex{Float64}
      L::Array{Complex{Float64},2} = [ 𝐋[i,j](s) for i = 1:N, j = 1:N ]
      A::Array{Complex{Float64},2} = [ 𝐀[i,j](s) for i = 1:N, j = 1:N ]

      exp(-σ[i,j]*s) / (2π*s^2) * (L * inv((1 + α*s)*I - A))[i,j]
    end

    ret[i,j] = ret[j,i] = G
  end

  return ret
end

@doc doc"""
structurefactor(ahs::AHS)

Returns a structure factor of additive hard-sphere system in wave-number space.
""" ->
function structurefactor(ahs::AHS)::Array{Function,2}
  N::Int = ncomp(ahs)
  ρ::Vector{Float64} = numberdensity(ahs)

  # g(r)*r in Laplace space
  G::Array{Function,2} = paircorrelationlaplace(ahs)

  I = eye(N) # identity matrix

  # matrix of structure factors
  ret = Array{Function,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    F(s::Complex{Float64})::Complex{Float64} = (s^2 * G[i,j](s) - 1) / s^3
    hᵢⱼ(q::Float64)::Float64 = -4π * real(F(im*q))

    S(q::Real)::Float64 = I[i,j] + √(ρ[i]*ρ[j]) * hᵢⱼ(convert(Float64, q))

    ret[i,j] = ret[j,i] = S
  end

  return ret
end

function paircorrelation(ahs::AHS)::Array{Function,2}
  N::Int = ncomp(ahs)
  σ::Array{Float64,2} = hsdiameter(ahs)
  ρ::Vector{Float64} = numberdensity(ahs)

  # g(r)*r in Laplace space
  G::Array{Function,2} = paircorrelationlaplace(ahs)

  # g(r) at r = σ
  g_σ::Array{Float64,2} = contactvalue(ahs)

  ret = Array{Function,2}(N,N)

  for i in 1:N, j in 1:N
    i > j && continue

    function g(r)::Float64
      if r < σ[i,j]
        0
      elseif r == σ[i,j]
        g_σ[i,j]
      else
        nterm = trunc(Int, 2048/sqrt(σ[i,j]))
        inverselaplace(G[i,j], r, nterm, 30) / r
      end
    end

    ret[i,j] = ret[j,i] = g
  end

  return ret
end

function cavityfunction(ahs::AHS)
  N = ncomp(ahs)
  δ = eye(N,N)
  σ = hsdiameter(ahs)
  ρ = numberdensity(ahs)
  g₀ = contactvalue(ahs)

  M = chemicalpotential(ahs)
  S = [ -π*δ[i,j] * sum(ρ[k]*σ[k,j]^2*g₀[k,j] for k in 1:N) for i in 1:N, j in 1:N ]
  Y = contactvalue(ahs)
  Y′ = contactgradient(ahs)
  a = [ abs(σ[i,i] - σ[j,j]) / 2 for i in 1:N, j in 1:N ]

  ret = Matrix{Function}(N,N)
  for i in 1:N, j in 1:N
    i > j && continue
    s = σ[i,j] - a[i,j]
    α = Vector{Float64}(4)
    α[1] = M[i]
    α[2] = S[i,j]
    α[3] = s^-2 * ( 3log(Y[i,j]) - s*Y′[i,j]/Y[i,j] - 3α[1] - 2s*α[2] )
    α[4] = s^-3 * ( -2log(Y[i,j]) + s*Y′[i,j]/Y[i,j] + 2α[1] + s*α[2] )
    function y(r::Real)::Float64
      if r < a[i,j]
        exp(M[i])
      elseif r < σ[i,j]
        exp(sum( α[k]*(r - a[i,j])^(k-1) for k in 1:4 ))
      else
        paircorrelation(ahs)[i,j](r)
      end
    end
    ret[i,j] = ret[j,i] = y
  end

  return ret
end

# Ref: K. Hoshino. and W. H. Young: J. Phys. F: Metal Phys. 10 (1980) 1365-1374.
function entropy(ahs::AHS)::Float64
  @attach(ahs, σ)

  N::Int = ncomp(ahs)
  c::Vector{Float64} = composition(ahs)
  η::Float64= totalpackingfraction(ahs)
  ζ::Float64 = 1 / (1 - η)

  y₁::Float64 = 0.
  y₂::Float64 = 0.

  for i in 1:N, j in 1:N
    i ≥ j && continue

    σ_sum = sum(c .* σ)

    y₁ += c[i]*c[j] * (σ[i] + σ[j]) * (σ[i] - σ[j])^2 / σ_sum^3
    y₂ += c[i]*c[j] * σ[i]*σ[j] * (σ[i] - σ[j])^2 / σ_sum^6 * sum(c .* σ.^2)
  end

  S_η = - (ζ - 1)*(ζ + 3)

  S_σ = (3/2*(ζ^2 - 1) - 1/2*(ζ - 1)*(ζ - 3) - log(ζ))*y₁ + (3/2*(ζ - 1)^2 - 1/2*(ζ - 1)*(ζ - 3) - log(ζ))*y₂

  S = S_η + S_σ
end

# Perturbation-independent part of internal energy
function internal(ahs::AHS)::Float64
  U = 0
end


text:
"""Julia program to implement Subset Sum problem.
Given a set, we have to find if there exist a subset with the given sum.
The given problem can be solved using Dynamic Programming.
"""

function subset_sum(set, n, S)

    # Initialize the entire 2D matrix with 0 as value, signifying false.
    dp = zeros(Int, n + 1, S + 1)

    # When the sum is zero, the answer is always true, by forming an empty subset.
    for i in 1:n + 1
        dp[i,1] = 1
    end

    # When the given set is empty and the required sum is not zero, then such a subset cannot e formed.
    for i in 1:S + 1
        dp[1,i] = 0
    end

    """Tabulating the dp array in Bottom-Up manner. In such a way that each cell of the dp[i, j]
       array contains whether a subset of sum j can be formed with elements in the given set till i"""
    for i in 2:n + 1
        for j in 2:S + 1
            if j < (set[i - 1] + 1)
                dp[i,j] = dp[(i - 1),j]
            else
                if (dp[(i - 1),j] == 1)
                    dp[i, j] = 1
                elseif (dp[(i - 1),(j - set[i - 1])] == 1)
                    dp[i, j] = 1
                else
                    dp[i, j] = 0
                end
            end
        end
    end
    return dp[n + 1,S + 1]
end

print("What is the size of the set? ")
n = readline()
n = parse(Int, n)
arr = Int[]
if (n > 0)
    print("Enter the numbers: ")
    arr = [parse(Int, num) for num in split(readline())]
else
    print("The set is Empty!!!")
    exit()
end
print("Enter the sum? ")
S = readline()
S = parse(Int, S)
res = subset_sum(arr, n, S)
if (res == 1)
    print("A subset with the given sum is present in the array.")
else
    print("A subset with the given sum is not present in the array.")
end

"""
Time Complexity - O(n * S), where 'n' is the size of the array and 'S is the Sum
Space Complexity - O(n * S)

SAMPLE INPUT AND OUTPUT

SAMPLE I

What is the size of the set? 3
Enter the numbers: 3 2 1
Enter the sum? 55
A subset with the given sum is not present in the array.

SAMPLE II

What is the size of the set? 4
Enter the numbers: 10 20 30 40
Enter the sum? 60
A subset with the given sum is present in the array.
"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
include("PUnits.jl")
include("Uncertainty.jl")
module Quantities
using PUnits, Uncertainty
# adding methods to:
import Base: promote_rule, convert, show, sqrt, +, *, -, /, ^, .*, ./, .^, ==, getindex, setindex!, size, ndims, endof, length, isapprox,
<, >, >=, <=, .!=, .==, .<, .>, .>=, .<=

typealias QValue  Union{Number, AbstractArray, Uncertain} # things that go inside a Quantity
# Quantity is where most of the action happens
# Quantity combines a value with some units
# Quantities can be reduced to base units via asbase, which uses the UnitSytem Dict
# all other conversions are call on asbase
type Quantity{T<:QValue}
    value::T
    unit::Unit
end
function Quantity_(x::QValue, y::Unit)
    if asbase(y) == Unitless
        return asbase(Quantity(x,y))
    end
    Quantity(x,y)
end
Quantity(x::Quantity, y::Unit) = error("Quantity{Quantity} not allowed")
# UnitSytem by example with SI units
# UnitSystem["J"] = a quantity with other units that is equal to a Joule
# UnitSystem["m"] = 0 which markes it as a non-prefixed base unit
UnitSystem = Dict{String, Union{Quantity, Int}}() # I should figure out how to give DefaultUnitSystem a type
reset_unit_system() = [pop!(UnitSystem, k) for (k,v) in UnitSystem]
function QUnit(x::AbstractString; prefix=0)
    return Quantity(1,Unit(x,prefix))
end
function BaseUnit(x::AbstractString; prefix=0, latex="", system=UnitSystem)
    system[x] = prefix # base units have their prefix as the dict value
    return Quantity(1,Unit(x,prefix))
end
function DerivedUnit(x::AbstractString, y::Quantity; latex="", system=UnitSystem)
    x=convert(String, x)
    system[x] = y
    return Quantity(1,Unit(x))
end
function asbase(x::Quantity, system=UnitSystem)
    x, x.unit
    prefactor, prefixless_unit = remove_prefix(x.unit)
    out=prefactor*x.value
    for (unitsymbol,n) in prefixless_unit.d
        next = system[unitsymbol.sym]
        if isa(next,Quantity)
            out *= asbase(next, system)^n
        elseif isa(next,Int) # next is the prefix of a base unit
            out *= (10.^-(next*n))*QUnit(unitsymbol.sym, prefix=next)^n
        else
            throw(ValueError())
        end
    end
    return out
end

function asbase(x::Unit, system=UnitSystem)
    out=Unitless
    for (unitsymbol,n) in x.d
        next = system[unitsymbol.sym]
        if isa(next,Quantity)
            out *=  asbase(next.unit, system)^n
        elseif isa(next,Int) # next is the prefix of a base unit
            out *= Unit(unitsymbol.sym,next)^n
        else
            throw(ValueError())
        end
    end
    return out
end
function as(from::Quantity, to::Quantity, system=UnitSystem) # warning, ignores the value of to, only looks at the units
    out = Quantity(1, to.unit)
    f,t = asbase(from, system), asbase(out, system)
    f.unit == t.unit ? out *= f.value./t.value : error("incompatible base units $(f.unit) and $(t.unit)")
    return out
end
# this is purposfully the only way to use a Prefix
*(x::Prefix, y::Quantity) = length(y.unit.d) == 1 ? Quantity(y.value, x*y.unit) : error("Prefix*Quantity only defined for Quantities with only one unit")
function *{T,S}(x::Quantity{T}, y::Quantity{S})
    xp,xu = remove_prefix(x.unit)
    yp,yu = remove_prefix(y.unit)
    Quantity_(xp*yp*x.value*y.value, xu*yu)
end
*(x::Quantity, y::QValue) = Quantity_(x.value*y, x.unit)
*(x::QValue, y::Quantity) = y*x
function .*{T,S}(x::Quantity{T}, y::Quantity{S})
    xp,xu = remove_prefix(x.unit)
    yp,yu = remove_prefix(y.unit)
    Quantity_(xp*yp*(x.value.*y.value), xu*yu)
end
.*(x::Quantity, y::QValue) = Quantity_(x.value.*y, x.unit)
.*(x::QValue, y::Quantity) = y.*x
function /{T,S}(x::Quantity{T}, y::Quantity{S})
    xp,xu = remove_prefix(x.unit)
    yp,yu = remove_prefix(y.unit)
    Quantity_((xp/yp)*(x.value/y.value), xu/yu)
end
/(x::Quantity, y::QValue) = Quantity_(x.value/y, x.unit)
/(x::QValue, y::Quantity) = Quantity_(x/y.value, y.unit^-1)
function ./{T,S}(x::Quantity{T}, y::Quantity{S})
    xp,xu = remove_prefix(x.unit)
    yp,yu = remove_prefix(y.unit)
    Quantity_((xp/yp)*(x.value./y.value), xu/yu)
end
./(x::Quantity, y::QValue) = Quantity_(x.value./y, x.unit)
./(x::QValue, y::Quantity) = Quantity_(x/y.value, y.unit^-1)
+{T,S}(x::Quantity{T}, y::Quantity{S}) = Quantity_(x.value + as(y,x).value, x.unit)
-{T,S}(x::Quantity{T}, y::Quantity{S}) = Quantity_(x.value - as(y,x).value, x.unit)
-{T}(x::Quantity{T}) = Quantity_(-x.value, x.unit)
^{T}(x::Quantity{T}, y::Rational) = Quantity_(x.value^convert(AbstractFloat,y), x.unit^convert(AbstractFloat,y))
^{T}(x::Quantity{T}, y::Integer) = Quantity_(x.value^convert(AbstractFloat,y), x.unit^convert(AbstractFloat,y))
^{T}(x::Quantity{T}, y::Number) = Quantity_(x.value^convert(AbstractFloat,y), x.unit^convert(AbstractFloat,y))
.^{T}(x::Quantity{T}, y::Rational) = Quantity_(x.value.^convert(AbstractFloat,y), x.unit.^convert(AbstractFloat,y))
.^{T}(x::Quantity{T}, y::Integer) = Quantity_(x.value.^convert(AbstractFloat,y), x.unit.^convert(AbstractFloat,y))
.^{T}(x::Quantity{T}, y::Number) = Quantity_(x.value.^convert(AbstractFloat,y), x.unit.^convert(AbstractFloat,y))
for f in (:(==), :<, :>, :>=, :<=, :!=, :(.==), :.<, :.>, :.>=, :.<=, :.!=, :isapprox)
    @eval begin function ($f)(x::Quantity, y::Quantity)
                    a = asbase(x)
                    b = asbase(y)
                    assert(a.unit == b.unit)
                    ($f)(a.value,b.value)
                end end
end
sqrt(x::Quantity) = Quantity_(sqrt(x.value), x.unit^.5)
getindex(x::Quantity, y...) = Quantity_(getindex(x.value, y...),x.unit)
setindex!(x::Quantity, y::Quantity, z...) = setindex!(x.value, as(y,x).value, z...)
setindex!(x::Quantity, y, z...) = error("x[z]=y requires same units, x.unit=$(x.unit), y has no units. use x[z] = y*same_units or x.value[z] = y instead")
size(x::Quantity) = size(x.value)
ndims(x::Quantity) = ndims(x.value)
endof(x::Quantity) = endof(x.value)
length(x::Quantity) = length(x.value)
Base.round(x::Quantity, y) = Quantity(round(x.value,y), x.unit)

function show{T}(io::IO, x::Quantity{T})
    show(io, x.value)
    print(io, " ")
    show(io, x.unit)
end

export BaseUnit, DerivedUnit, Prefix, asbase, as, Uncertain

end #end module


text:
#=
map_many.jl
2019-06-13, Jeff Fessler, University of Michigan
=#

export map_many

"""
    y = map_many(fun::Function, x::AbstractArray{<:Any}, idim::Dims)

apply a function `fun` to leading slices of input `x`;
cousin of `mapslices`

in
- `fun::Function` maps input of size `idim` to output of some size `odim`
- `x [idim ldim]`

out
- `y [odim ldim]`

Example: if `fun` maps array of size (1,2) to array of size (3,4,5)
and if input `x` has size (1,2,7,8)
then output `y` will have size (3,4,5,7,8)
where `y[:,:,:,i,j] = fun(x[:,:,i,j])`
"""
function map_many(fun::Function, x::AbstractArray{<:Any}, idim::Dims)
	xdim = size(x)
	D = length(idim)
	ndims(x) < D && throw(DimensionMismatch("ndims(x) < D=$D)"))
	xdim[1:D] != idim && throw(DimensionMismatch("size(x) vs idim"))

	if xdim == idim
		return fun(x)
	end

	ldim = xdim[(D+1):end]
	L = prod(ldim) # *ldim
	x = reshape(x, prod(idim), :) # [*idim L]
	tmp = fun(reshape((@view x[:,1]), idim)) # [odim]
	odim = size(tmp)
	T = eltype(tmp)
	out = Array{T}(undef, prod(odim), L) # [*odim L]
	out[:,1] = vec(tmp)
	for l=2:L
		out[:,l] = vec(fun(reshape((@view x[:,l]), idim)))
	end
	return reshape(out, odim..., ldim...) # [odim ldim]
end


text:
"This module contains functions that can only be applied to CSV files."
module CSV

using ..Constants
using ..Exceptions
using ..PDF

export _get_authors_from_CSV,
       _get_highlights_from_CSV,
       _get_titles_from_CSV,
       get_all,
       get_locations,
       get_notes,
       get_titles,
       import_highlights,
       initialize

# Load the functions

# Initialize a table
include("Functions/initialize.jl")

# Check if the structure of the table is correct
include("Functions/_check.jl")

# Import the highlights from the PDF
include("Functions/import_highlights.jl")

# Get the column values from the table
include("Functions/_get_authors_from_CSV.jl")
include("Functions/_get_highlights_from_CSV.jl")
include("Functions/_get_titles_from_CSV.jl")
include("Functions/get_all.jl")
include("Functions/get_notes.jl")
include("Functions/get_locations.jl")

end


text:
# # Binary Classification with Laplace approximation
#
# This example demonstrates how to carry out non-conjugate Gaussian process
# inference using the Laplace approximation.
#
# For a basic introduction to the functionality of this library, please refer
# to the [User Guide](@ref).
#
# ## Setup

using ApproximateGPs
using LinearAlgebra
using Distributions
using LogExpFunctions: logistic, softplus, invsoftplus
using Zygote
using Optim

using Plots
default(; legend=:outertopright, size=(700, 400))

using Random
Random.seed!(1);

# ## Generate training data
#
# We create a binary-labelled toy dataset:

Xgrid = -4:0.1:29  # for visualization
X = range(0, 23.5; length=48)  # training inputs
f(x) = 3 * sin(10 + 0.6x) + sin(0.1x) - 1  # latent function
fs = f.(X)  # latent function values at training inputs

lik = BernoulliLikelihood()  # has logistic invlink by default
## could use other invlink, e.g. normcdf(f) = cdf(Normal(), f)

invlink = lik.invlink  # logistic function
ps = invlink.(fs)  # probabilities at the training inputs
Y = [rand(Bernoulli(p)) for p in ps]  # observations at the training inputs
## could do this in one call as `Y = rand(lik(fs))`

function plot_data()
    plot(; xlims=extrema(Xgrid), xticks=0:6:24)
    plot!(Xgrid, invlink ∘ f; label="true probabilities")
    return scatter!(X, Y; label="observations", color=3)
end

plot_data()

# ## Create a latent GP
#
# Here we write a function that creates our latent GP prior, given the
# hyperparameter vector `theta`. Compared to a "vanilla" GP, the `LatentGP`
# requires a function or functor that maps from the latent GP `f` to the
# distribution of observations `y`. This functor is commonly called "the
# likelihood".

function build_latent_gp(theta)
    ## `theta` is unconstrained, but kernel variance and lengthscale must be positive:
    variance = softplus(theta[1])
    lengthscale = softplus(theta[2])

    kernel = variance * with_lengthscale(SqExponentialKernel(), lengthscale)

    dist_y_given_f = BernoulliLikelihood()  # has logistic invlink by default
    ## We could also be explicit and define it as a function:
    ## dist_y_given_f(f) = Bernoulli(invlink(f))

    jitter = 1e-8  # required for numeric stability
    return LatentGP(GP(kernel), dist_y_given_f, jitter)
end;

# We define a latent GP at our initial hyperparameter values, here with
# variance 1.0 and lengthscale 5.0:

theta0 = [invsoftplus(1.0), invsoftplus(5.0)]

lf = build_latent_gp(theta0)

lf.f.kernel

# We can now compute the Laplace approximation ``q(f)`` to the true posterior
# ``p(f | y)``:

f_post = posterior(LaplaceApproximation(), lf(X), Y)

# This finds a mode of the posterior (for the given values of the
# hyperparameters) using iterated Newton's method (i.e. solving an optimisation
# problem) and then constructs a Gaussian approximation to the posterior by
# matching the curvature at the mode.

# Let's plot samples from this approximate posterior:

function plot_samples!(Xgrid, fpost; samples=100, color=2)
    fx = fpost(Xgrid, 1e-8)
    fsamples = rand(fx, samples)
    plot!(Xgrid, invlink.(fsamples); color, alpha=0.2, label="")
    return plot!(Xgrid, invlink.(mean(fx)); color, lw=2, label="posterior fit")
end

p1 = plot_data()
plot_samples!(Xgrid, f_post)

# We can improve this fit by optimising the hyperparameters. For exact Gaussian
# process regression, the maximization objective is the marginal likelihood.
# Here, we can only optimise an _approximation_ to the marginal likelihood.

# ## Optimise the hyperparameters
#
# ApproximateGPs provides a convenience function `build_laplace_objective` that
# constructs an objective function for optimising the hyperparameters, based on
# the Laplace approximation to the log marginal likelihood.

objective = build_laplace_objective(build_latent_gp, X, Y);

# We pass this objective to Optim.jl's LBFGS optimiser:

training_results = Optim.optimize(
    objective, θ -> only(Zygote.gradient(objective, θ)), theta0, LBFGS(); inplace=false
)

# Now that we have (hopefully) better hyperparameter values, we need to construct a LatentGP prior with these values:

lf2 = build_latent_gp(training_results.minimizer)

lf2.f.kernel

# Finally, we need to construct again the (approximate) posterior given the
# observations for the latent GP with optimised hyperparameters:

f_post2 = posterior(LaplaceApproximation(; f_init=objective.f), lf2(X), Y)

# By passing `f_init=objective.f` we let the Laplace approximation "warm-start"
# at the last point of the inner-loop Newton optimisation; `objective.f` is a
# field on the `objective` closure.

# Let's plot samples from the approximate posterior for the optimised hyperparameters:

p2 = plot_data()
plot_samples!(Xgrid, f_post2)


text:
const Tuple3F64 = Tuple{Float64,Float64,Float64}
const Tuple3I64 = Tuple{Int64,Int64,Int64}


text:
function _permutedims{dim}(S::FourthOrderTensor{dim}, idx::NTuple{4,Int})
    sort([idx...]) == [1,2,3,4] || throw(ArgumentError("Missing index."))
    neworder = sortperm([idx...])
    f = (i,j,k,l) -> S[[i,j,k,l][neworder]...]
    return Tensor{4,dim}(f)
end

@testset "tensor operations" begin;
for T in (Float32, Float64), dim in (1,2,3)
AA = rand(Tensor{4, dim, T})
BB = rand(Tensor{4, dim, T})
A = rand(Tensor{2, dim, T})
B = rand(Tensor{2, dim, T})
a = rand(Tensor{1, dim, T})
b = rand(Tensor{1, dim, T})

AA_sym = rand(SymmetricTensor{4, dim, T})
BB_sym = rand(SymmetricTensor{4, dim, T})
A_sym = rand(SymmetricTensor{2, dim, T})
B_sym = rand(SymmetricTensor{2, dim, T})

i,j,k,l = rand(1:dim,4)

@testset "double contraction" begin
    # 4 - 4
    # Value tests
    @test vec(dcontract(AA, BB)) ≈ vec(reshape(vec(AA), (dim^2, dim^2)) * reshape(vec(BB), (dim^2, dim^2)))
    @test vec(dcontract(AA_sym, BB)) ≈ vec(reshape(vec(AA_sym), (dim^2, dim^2)) * reshape(vec(BB), (dim^2, dim^2)))
    @test vec(dcontract(AA, BB_sym)) ≈ vec(reshape(vec(AA), (dim^2, dim^2)) * reshape(vec(BB_sym), (dim^2, dim^2)))
    @test vec(dcontract(AA_sym, BB_sym)) ≈ vec(reshape(vec(AA_sym), (dim^2, dim^2)) * reshape(vec(BB_sym), (dim^2, dim^2)))
    @test dcontract(convert(Tensor, AA_sym), convert(Tensor, BB_sym)) ≈ dcontract(AA_sym, BB_sym)
    # Type tests
    @test isa(dcontract(AA, BB), Tensor{4, dim, T})
    @test isa(dcontract(AA_sym, BB), Tensor{4, dim, T})
    @test isa(dcontract(AA, BB_sym), Tensor{4, dim, T})
    @test isa(dcontract(AA_sym, BB_sym), SymmetricTensor{4, dim, T})

    # 2 - 4
    # Value tests
    @test dcontract(AA, A) ≈ reshape(reshape(vec(AA), (dim^2, dim^2)) * reshape(vec(A), (dim^2,)), dim, dim)
    @test dcontract(AA_sym, A) ≈ reshape(reshape(vec(AA_sym), (dim^2, dim^2)) * reshape(vec(A), (dim^2,)), dim, dim)
    @test dcontract(AA, A_sym) ≈ reshape(reshape(vec(AA), (dim^2, dim^2)) * reshape(vec(A_sym), (dim^2,)), dim, dim)
    @test dcontract(AA_sym, A_sym) ≈ reshape(reshape(vec(AA_sym), (dim^2, dim^2)) * reshape(vec(A_sym), (dim^2,)), dim, dim)
    @test dcontract(A, AA) ≈ reshape(reshape(vec(AA), (dim^2, dim^2))' * reshape(vec(A), (dim^2,)), dim, dim)
    @test dcontract(A_sym, AA) ≈ reshape(reshape(vec(AA), (dim^2, dim^2))' * reshape(vec(A_sym), (dim^2,)), dim, dim)
    @test dcontract(A, AA_sym) ≈ reshape(reshape(vec(AA_sym), (dim^2, dim^2))' * reshape(vec(A), (dim^2,)), dim, dim)
    @test dcontract(A_sym, AA_sym) ≈ reshape(reshape(vec(AA_sym), (dim^2, dim^2))' * reshape(vec(A_sym), (dim^2,)), dim, dim)
    @test dcontract(convert(Tensor, AA_sym), convert(Tensor, A_sym)) ≈ dcontract(AA_sym, A_sym)
    # Type tests
    @test isa(dcontract(AA, A), Tensor{2, dim, T})
    @test isa(dcontract(AA_sym, A), SymmetricTensor{2, dim, T})
    @test isa(dcontract(AA, A_sym), Tensor{2, dim, T})
    @test isa(dcontract(AA_sym, A_sym), SymmetricTensor{2, dim, T})
    @test isa(dcontract(A, AA), Tensor{2, dim, T})
    @test isa(dcontract(A_sym, AA), Tensor{2, dim, T})
    @test isa(dcontract(A, AA_sym), SymmetricTensor{2, dim, T})
    @test isa(dcontract(A_sym, AA_sym), SymmetricTensor{2, dim, T})

    # 2 - 2
    # Value tests
    @test dcontract(A, B) ≈ sum(vec(A) .* vec(B))
    @test dcontract(A_sym, B) ≈ sum(vec(A_sym) .* vec(B))
    @test dcontract(A, B_sym) ≈ sum(vec(A) .* vec(B_sym))
    @test dcontract(A_sym, B_sym) ≈ sum(vec(A_sym) .* vec(B_sym))
    # Type tests
    @test isa(dcontract(A, B), T)
    @test isa(dcontract(A_sym, B), T)
    @test isa(dcontract(A, B_sym), T)
    @test isa(dcontract(A_sym, B_sym), T)
end # of testset

@testset "outer product" begin
    # Value tests
    @test otimes(a, b) ≈ Array(a) * Array(b)'
    @test reshape(vec(otimes(A, B)), dim^2, dim^2) ≈ vec(A) * vec(B)'
    @test reshape(vec(otimes(A_sym, B)), dim^2, dim^2) ≈ vec(A_sym) * vec(B)'
    @test reshape(vec(otimes(A, B_sym)), dim^2, dim^2) ≈ vec(A) * vec(B_sym)'
    @test reshape(vec(otimes(A_sym, B_sym)), dim^2, dim^2) ≈ vec(A_sym) * vec(B_sym)'

    # Type tests
    @test isa(otimes(a, b), Tensor{2, dim, T})
    @test isa(otimes(A, B), Tensor{4, dim, T})
    @test isa(otimes(A_sym, B), Tensor{4, dim, T})
    @test isa(otimes(A, B_sym), Tensor{4, dim, T})
    @test isa(otimes(A_sym, B_sym), SymmetricTensor{4, dim, T})
end # of testset

@testset "dot products" begin
    # 1 - 2
    # Value tests
    @test dot(a, b) ≈ sum(Array(a) .* Array(b))
    @test dot(A, b) ≈ Array(A) * Array(b)
    @test dot(A_sym, b) ≈ Array(A_sym) * Array(b)
    @test dot(a, B) ≈ Array(B)' * Array(a)
    @test dot(a, B_sym) ≈ Array(B_sym)' * Array(a)

    # Type tests
    @test isa(dot(a, b), T)
    @test isa(dot(A, b), Vec{dim, T})
    @test isa(dot(A_sym, b), Vec{dim, T})
    @test isa(dot(b, A), Vec{dim, T})
    @test isa(dot(b, A_sym), Vec{dim, T})

    # 2 - 2
    # Value tests
    @test dot(A, B) ≈ Array(A) * Array(B)
    @test dot(A_sym, B) ≈ Array(A_sym) * Array(B)
    @test dot(A, B_sym) ≈ Array(A) * Array(B_sym)
    @test dot(A_sym, B_sym) ≈ Array(A_sym) * Array(B_sym)

    @test tdot(A) ≈ Array(A)' * Array(A)
    @test tdot(A_sym) ≈ Array(A_sym)' * Array(A_sym)
    @test dott(A) ≈ Array(A) * Array(A)'
    @test dott(A_sym) ≈ Array(A_sym) * Array(A_sym)'
    @test tdot(A) ≈ dott(transpose(A))
    @test tdot(transpose(A)) ≈ dott(A)
    @test tdot(A_sym) ≈ dott(transpose(A_sym))
    @test tdot(transpose(A_sym)) ≈ dott(A_sym)

    # Type tests
    @test isa(dot(A, B), Tensor{2, dim, T})
    @test isa(dot(A_sym, B), Tensor{2, dim, T})
    @test isa(dot(A, B_sym), Tensor{2, dim, T})
    @test isa(dot(A_sym, B_sym), Tensor{2, dim, T})

    @test isa(tdot(A), SymmetricTensor{2, dim, T})
    @test isa(tdot(A_sym), SymmetricTensor{2, dim, T})
    @test isa(dott(A), SymmetricTensor{2, dim, T})
    @test isa(dott(A_sym), SymmetricTensor{2, dim, T})
end # of testset

@testset "symmetric/skew-symmetric" begin
    if dim == 1 # non-symmetric tensors are symmetric
        @test issymmetric(A)
        @test issymmetric(AA)
    elseif dim != 1
        @test !issymmetric(A)
        @test !issymmetric(AA)
        @test !ismajorsymmetric(AA)
        @test !isminorsymmetric(AA)
        @test_throws InexactError convert(typeof(A_sym),A)
        @test_throws InexactError convert(typeof(AA_sym),AA)
    end
    @test issymmetric(A_sym)
    @test issymmetric(AA_sym)
    @test isminorsymmetric(AA_sym)
    @test issymmetric(symmetric(A))
    @test issymmetric(A + A.')

    @test symmetric(A) ≈ 0.5(A + A.')
    @test symmetric(A_sym) ≈ A_sym
    @test isa(symmetric(A), SymmetricTensor{2, dim, T})
    @test convert(typeof(A_sym),convert(Tensor,symmetric(A))) ≈ symmetric(A)

    @test symmetric(AA) ≈ minorsymmetric(AA)
    @test minorsymmetric(AA)[i,j,k,l] ≈ minorsymmetric(AA)[j,i,l,k]
    @test issymmetric(convert(Tensor,minorsymmetric(AA)))
    @test isminorsymmetric(convert(Tensor,minorsymmetric(AA)))
    @test symmetric(AA_sym) ≈ minorsymmetric(AA_sym)
    @test minorsymmetric(AA_sym)[i,j,k,l] ≈ minorsymmetric(AA_sym)[j,i,l,k]
    @test minorsymmetric(AA_sym) ≈ AA_sym
    @test issymmetric(convert(Tensor,minorsymmetric(AA_sym)))
    @test isminorsymmetric(convert(Tensor,minorsymmetric(AA_sym)))

    @test isa(minorsymmetric(AA), SymmetricTensor{4, dim, T})
    @test convert(typeof(AA_sym),convert(Tensor,minorsymmetric(AA))) ≈ minorsymmetric(AA)

    @test majorsymmetric(AA)[i,j,k,l] ≈ 0.5*(AA[i,j,k,l] + AA[k,l,i,j])
    @test majorsymmetric(AA)[i,j,k,l] ≈ majorsymmetric(AA)[k,l,i,j]
    @test ismajorsymmetric(majorsymmetric(AA))
    @test majorsymmetric(AA_sym)[i,j,k,l] ≈ 0.5*(AA_sym[i,j,k,l] + AA_sym[k,l,i,j])
    @test majorsymmetric(AA_sym)[i,j,k,l] ≈ majorsymmetric(AA_sym)[k,l,i,j]
    @test ismajorsymmetric(majorsymmetric(AA_sym))
    @test isa(majorsymmetric(AA), Tensor{4, dim, T})
    @test isa(majorsymmetric(AA_sym), Tensor{4, dim, T})

    @test skew(A) ≈ 0.5(A - A.')
    @test isa(skew(A), Tensor{2, dim, T})
    @test skew(A_sym) ≈ zero(A_sym)
    @test isa(skew(A_sym), Tensor{2, dim, T})

    # Identities
    @test A ≈ symmetric(A) + skew(A)
    @test skew(A) ≈ -skew(A).'
    @test trace(skew(A)) ≈ 0.0
    @test trace(symmetric(A)) ≈ trace(A)
end # of testset

@testset "transpose" begin
    @test transpose(a) ≈ a' ≈ a
    @test isa(transpose(a), Vec{dim, T})
    @test transpose(A) ≈ Array(A).'
    @test transpose(transpose(A)) ≈ A
    @test transpose(A_sym) ≈ Array(A_sym).'
    @test transpose(A_sym) ≈ A_sym
    @test transpose(transpose(A_sym)) ≈ A_sym

    @test transpose(AA) ≈ minortranspose(AA)
    @test AA[i,j,k,l] ≈ minortranspose(AA)[j,i,l,k]
    @test AA_sym[i,j,k,l] ≈ minortranspose(AA_sym)[j,i,l,k]
    @test minortranspose(AA_sym) ≈ AA_sym
    @test minortranspose(minortranspose(AA)) ≈ AA
    @test minortranspose(minortranspose(AA_sym)) ≈ AA_sym

    @test majortranspose(majortranspose(AA)) ≈ AA
    @test majortranspose(majortranspose(AA_sym)) ≈ AA_sym
    @test AA[i,j,k,l] ≈ majortranspose(AA)[k,l,i,j]
    @test AA_sym[i,j,k,l] ≈ majortranspose(AA_sym)[k,l,i,j]
    @test typeof(majortranspose(AA_sym)) <: Tensor{4,dim}

    @test minortranspose(AA) ≈ _permutedims(AA,(2,1,4,3))
    @test minortranspose(AA_sym) ≈ _permutedims(AA_sym,(2,1,4,3))
    @test majortranspose(AA) ≈ _permutedims(AA,(3,4,1,2))
    @test majortranspose(AA_sym) ≈ _permutedims(AA_sym,(3,4,1,2))
end # of testset

@testset "cross product" begin
    @test a × a ≈ Vec{3, T}((0.0,0.0,0.0))
    @test a × b ≈ -b × a
    if dim == 2
        ad = Vec{2, T}((1.0,0.0))
        ad2 = Vec{2, T}((0.0,1.0))
        @test ad × ad2 ≈ Vec{3, T}((0.0, 0.0, 1.0))
    end
    if dim == 3
        ad = Vec{3, T}((1.0,0.0,0.0))
        ad2 = Vec{3, T}((0.0,1.0,0.0))
        @test ad × ad2 ≈ Vec{3, T}((0.0, 0.0, 1.0))
    end
end # of testset

@testset "special" begin
    AAT = Tensor{4, dim, T}((i,j,k,l) -> AA_sym[i,l,k,j])
    @test AAT ⊡ (b ⊗ a) ≈ dotdot(a, AA_sym, b)
end # of testset

@testset "rotation" begin
    x = eᵢ(Vec{3}, 1)
    y = eᵢ(Vec{3}, 2)
    z = eᵢ(Vec{3}, 3)

    @test rotate(z, z, rand()) ≈ z
    @test rotate(2*z, y, π/2) ≈ 2*x
    @test rotate(3*z, y, π) ≈ -3*z
    @test rotate(x+y+z, z, π/4) ≈ Vec{3}((0.0,√2,1.0))

    a = rand(Vec{3})
    b = rand(Vec{3})
    @test rotate(a, b, 0) ≈ a
    @test rotate(a, b, π) ≈ rotate(a, b, -π)
    @test rotate(a, b, π/2) ≈ rotate(a, -b, -π/2)
end # of testset

end
end # of testset


text:
################################################################################
# Variance estimate via OLS and QR decomposition.
################################################################################
function initvar(y::Vector, X::Matrix{T}) where T
    qrx  = qr(X)
    β    = inv(qrx.R) * qrx.Q' * y
    r    = y .- X * β
    sum(x -> x * x, r)/(length(r) - size(X, 2)), β
end
################################################################################
function nterms(lmm::LMM)
    nterms(lmm.mf)
end
function nterms(mf::ModelFrame)
    mf.schema.schema.count
end
function nterms(rhs::Union{Tuple{Vararg{AbstractTerm}}, Nothing, AbstractTerm})
    if isa(rhs, Term)
        p = 1
    elseif isa(rhs, Tuple)
        p = length(rhs)
    else
        p = 0
    end
    p
end
################################################################################
#                        VAR LINK
################################################################################

function vlink(σ::T) where T <: Real
    if σ < -21.0 return one(T)*7.582560427911907e-10 end #Experimental
    exp(σ)
end
function vlinkr(σ::T) where T <: Real
    log(σ)
end

function vlinksq(σ::T) where T <: Real
    σ*σ
end
function vlinksqr(σ::T) where T <: Real
    sqrt(σ)
end

function rholinkpsigmoid(ρ::T) where T <: Real
    return 1.0/(1.0 + exp(-ρ * 0.5))
end
function rholinkpsigmoidr(ρ::T) where T <: Real
    return - log(1.0/ρ - 1.0) / 0.5
end

function rholinksigmoid(ρ::T) where T <: Real
    return 1.0/(1.0 + exp(- ρ * 0.1)) * 2.0 - 1.0
end
function rholinksigmoidr(ρ::T) where T <: Real
    return - log(1.0/(ρ+1.0)*2.0 - 1.0) / 0.1
end

function rholinksqsigmoid(ρ::T) where T <: Real
    return ρ/sqrt(1.0 + (ρ)^2)
end
function rholinksqsigmoidr(ρ::T) where T <: Real
    return sign(ρ)*sqrt(ρ^2/(1.0 - ρ^2))
end

function rholinksigmoidatan(ρ::T) where T <: Real
    return atan(ρ)/pi*2.0
end
function rholinksigmoidatanr(ρ::T) where T <: Real
    return tan(ρ*pi/2.0)
end

################################################################################
################################################################################
function varlinkvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm)
    @inbounds for i = 1:length(v)
        if p[i] == :var
            if varlinkf == :exp
                v[i] = vlink(v[i])
            elseif varlinkf == :sq
                v[i] = vlinksq(v[i])
            elseif varlinkf == :identity
                v[i] = identity(v[i])
            end
        elseif p[i] == :rho
            if rholinkf == :sigm
                v[i] = rholinksigmoid(v[i])
            elseif rholinkf == :atan
                v[i] = rholinksigmoidatan(v[i])
            elseif rholinkf == :sqsigm
                v[i] = rholinksqsigmoid(v[i])
            elseif rholinkf == :psigm
                v[i] = rholinkpsigmoid(v[i])
            end
        end
    end
    v
end
function varlinkrvecapply!(v, p; varlinkf = :exp, rholinkf = :sigm)
    @inbounds for i = 1:length(v)
        if p[i] == :var
            if varlinkf == :exp
                v[i] = vlinkr(v[i])
            elseif varlinkf == :sq
                v[i] = vlinksqr(v[i])
            elseif varlinkf == :identity
                v[i] = identity(v[i])
            end
        elseif p[i] == :rho
            if rholinkf == :sigm
                v[i] = rholinksigmoidr(v[i])
            elseif rholinkf == :atan
                v[i] = rholinksigmoidatanr(v[i])
            elseif rholinkf == :sqsigm
                v[i] = rholinksqsigmoidr(v[i])
            elseif rholinkf == :psigm
                v[i] = rholinkpsigmoidr(v[i])
            end
        end
    end
    v
end
function varlinkvecapply(v, p; varlinkf = :exp, rholinkf = :sigm)
    s = similar(v)
    @inbounds for i = 1:length(v)
        if p[i] == :var
            if varlinkf == :exp
                s[i] = vlink(v[i])
            elseif varlinkf == :sq
                s[i] = vlinksq(v[i])
            elseif varlinkf == :identity
                s[i] = identity(v[i])
            end
        elseif p[i] == :rho
            if rholinkf == :sigm
                s[i] = rholinksigmoid(v[i])
            elseif rholinkf == :atan
                s[i] = rholinksigmoidatan(v[i])
            elseif rholinkf == :sqsigm
                s[i] = rholinksqsigmoid(v[i])
            elseif rholinkf == :psigm
                s[i] = rholinkpsigmoid(v[i])
            end
        else
            s[i] = v[i]
        end
    end
    s
end
################################################################################
function m2logreml(lmm)
    lmm.result.reml
end
function logreml(lmm)
    -m2logreml(lmm)/2
end
function m2logreml(lmm, theta)
    reml_sweep_β(lmm, LMMDataViews(lmm), theta)[1]
end
################################################################################

function optim_callback(os)
    false
end
################################################################################
"""
    gmatrix(lmm::LMM{T}, r::Int) where T
"""
function gmatrix(lmm::LMM{T}, r::Int) where T
    if isnothing(lmm.result.theta) error("No results or model not fitted!") end
    if r > length(lmm.covstr.random) error("Invalid random effect number: $(r)!") end
    G = zeros(T, lmm.covstr.q[r], lmm.covstr.q[r])
    gmat_switch!(G, lmm.result.theta, lmm.covstr, r)
    Symmetric(G)
end
"""
    rmatrix(lmm::LMM{T}, i::Int) where T
"""
function rmatrix(lmm::LMM{T}, i::Int) where T
    if !lmm.result.fit error("Model not fitted!") end
    if i > length(lmm.covstr.vcovblock) error("Invalid block number: $(i)!") end
    q    = length(lmm.covstr.vcovblock[i])
    R    = zeros(T, q, q)
    rmat_base_inc!(R, lmm.result.theta[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
    Symmetric(R)
end
"""
    vmatrix!(V, θ, lmm, i)

Update variance-covariance matrix V for i bolock. Upper triangular updated.
"""
function vmatrix!(V, θ, lmm, i)
    zgz_base_inc!(V, θ, lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
    rmat_base_inc!(V, θ[lmm.covstr.tr[end]], lmm.covstr, lmm.covstr.vcovblock[i], lmm.covstr.sblock[i])
end
function vmatrix(θ, lmm, i)
    V = zeros(length(lmm.covstr.vcovblock[i]), length(lmm.covstr.vcovblock[i]))
    vmatrix!(V, θ, lmm, i)
    Symmetric(V)
end
function vmatrix(lmm, i)
    vmatrix(lmm.result.theta, lmm, i)
end

function nblocks(lmm)
    return length(lmm.covstr.vcovblock)
end
"""
    hessian(lmm, theta)

Calculate Hessian matrix of REML for theta.
"""
function hessian(lmm, theta)
    #if !lmm.result.fit error("Model not fitted!") end
    vloptf(x) = reml_sweep_β(lmm, x, lmm.result.beta)[1]
    chunk  = ForwardDiff.Chunk{1}()
    hcfg   = ForwardDiff.HessianConfig(vloptf, theta, chunk)
    ForwardDiff.hessian(vloptf, theta, hcfg, Val{false}())
end
function hessian(lmm)
    if !lmm.result.fit error("Model not fitted!") end
    hessian(lmm, lmm.result.theta)
end
################################################################################

function get_symb(t::T; v = Vector{Symbol}(undef, 0)) where T <: Union{ConstantTerm, InterceptTerm, FunctionTerm}
    v
end
function get_symb(t::T; v = Vector{Symbol}(undef, 0)) where T <: Union{Term, CategoricalTerm}
    push!(v, t.sym)
    v
end
function get_symb(t::T; v = Vector{Symbol}(undef, 0)) where T <: InteractionTerm
    for i in t.terms
        get_symb(i; v = v)
    end
    v
end
function get_symb(t::T; v = Vector{Symbol}(undef, 0)) where T <: Tuple{Vararg{AbstractTerm}}
    for i in t
        get_symb(i; v = v)
    end
    v
end

################################################################################
# logdet with check
#=
function logdet_(C::Cholesky, noerror)
    dd = zero(real(eltype(C)))
    @inbounds for i in 1:size(C.factors,1)
        v = real(C.factors[i,i])
        if v > 0
            dd += log(v)
        else
            C.factors[i,i] = LOGDETCORR
            dd += log(real(C.factors[i,i]))
            noerror = false
        end
    end
    dd + dd, noerror
end
=#

function StatsModels.termvars(ve::VarEffect)
    termvars(ve.formula)
end
function StatsModels.termvars(ve::Vector{VarEffect})
    union(termvars.(ve)...)
end


text:
using Atom

# calling Atom's serve()

# using Pkg
# Pkg.add("Atom")
# Pkg.activate("..")

# Atom.handle("ping") do
#   "pong"
# end

Atom.serve(5601)
# This should connect to the elisp side server
# Atom.connect(5602)

# using MacroTools

# @destruct [word, path || nothing] = [Dict("word"=>"hello"), Dict("path"=>nothing)]

max_stars_repo_name: lihebi/julia-repl

max_issues_repo_name: lihebi/julia-repl

max_forks_repo_name: lihebi/julia-repl


text:
# Linear Radon transform modeling with a synthetic example.

using PyPlot, Seismic

close("all")

# 1- Make a 2D clean data with linear events
nt = 80
dt = 0.004;
nx = 32
dx = 10
h = collect(0:nx-1)*dx
d = zeros(Float32, nt, nx)
d[40, :] = 1.f0
d[25:nt, :] = d[25:nt, :] +  eye(nt-24, nx)
d[10:2:nt, :] = d[10:2:nt, :] + eye(40-4, nx)

# 2- Define linear moveouts p
pmin = 0
pmax = 2.5*dt/dx
dp = 2*dt/maximum(h)
p = collect(pmin:dp:pmax)

# 3- Transform to Radon domain
m = SeisRadonFreqInv(d, order="linear", dt=dt, h=h, p=p, flow=0, fhigh=125,
                     mu=1.e-5 )

# 4- Transform back to data domain
d2 = SeisRadonFreqFor(m, nt, order="linear", dt=dt, h=h, p=p, flow=0, fhigh=125)

# 5- Comput and print modeling error (%)
e2 = 100*maximum(abs(d-d2 ))/maximum(abs(d2))
println("Relative error = ", e2, " %")

# 6- Plot input data, Radon gather, modeled data and modeling Error
figure(1, figsize=(12,6));
subplot(141)
SeisPlot(d, title="Data", xlabel="Offset [m]", ylabel="Time [s]", vmin=-0.8,
         vmax = 0.8, ox=h[1], dx=dx, dy=dt, fignum=1)
subplot(142)
SeisPlot(m[1:nt, :], title="Radon gather", xlabel="Residual moveout [s]",
         ox=p[1], dx=p[2]-p[1], dy=dt, fignum=1)
subplot(143)
SeisPlot(d2, title="Modeled data", xlabel="Offset [m]", vmin = -0.8, vmax = 0.8,
         ox=h[1], dx=dx, dy=dt, fignum=1)
subplot(144)
SeisPlot(d-d2, title="Error", xlabel="Offset [m]", ox=h[1], dx=dx, dy=dt,
         fignum=1)
colorbar()
tight_layout()


text:
# Convenience functions - make a Random Number Generator object
mk_rng(rng::AbstractRNG) = rng
mk_rng(seed::Int) = MersenneTwister(seed)

function _split_mse{T<:Float64, U<:Real}(labels::Vector{T}, features::Matrix{U}, nsubfeatures::Int, rng)
    nr, nf = size(features)

    best = NO_BEST
    best_val = -Inf

    if nsubfeatures > 0
        r = randperm(rng, nf)
        inds = r[1:nsubfeatures]
    else
        inds = 1:nf
    end

    for i in inds
        # Sorting used to be performed only when nr <= 100, but doing it
        # unconditionally improved fitting performance by 20%. It's a bit of a
        # puzzle. Either it improved type-stability, or perhaps branch
        # prediction is much better on a sorted sequence.
        ord = sortperm(features[:,i])
        features_i = features[ord,i]
        labels_i = labels[ord]
        if nr > 100
            if VERSION >= v"0.4.0-dev"
                domain_i = quantile(features_i, linspace(0.01, 0.99, 99);
                                    sorted=true)
            else  # sorted=true isn't supported on StatsBase's Julia 0.3 version
                domain_i = quantile(features_i, linspace(0.01, 0.99, 99))
            end
        else
            domain_i = features_i
        end
        value, thresh = _best_mse_loss(labels_i, features_i, domain_i)
        if value > best_val
            best_val = value
            best = (i, thresh)
        end
    end

    return best
end

""" Finds the threshold to split `features` with that minimizes the
mean-squared-error loss over `labels`.

Returns (best_val, best_thresh), where `best_val` is -MSE """
function _best_mse_loss{T<:Float64, U<:Real}(labels::Vector{T}, features::Vector{U}, domain)
    # True, but costly assert. However, see
    # https://github.com/JuliaStats/StatsBase.jl/issues/164
    # @assert issorted(features) && issorted(domain)
    best_val = -Inf
    best_thresh = 0.0
    s_l = s2_l = zero(T)
    su = sum(labels)::T
    su2 = zero(T); for l in labels su2 += l*l end  # sum of squares
    nl = 0
    n = length(labels)
    i = 1
    # Because the `features` are sorted, below is an O(N) algorithm for finding
    # the optimal threshold amongst `domain`. We simply iterate through the
    # array and update s_l and s_r (= sum(labels) - s_l) as we go. - @cstjean
    @inbounds for thresh in domain
        while i <= length(labels) && features[i] < thresh
            l = labels[i]

            s_l += l
            s2_l += l*l
            nl += 1

            i += 1
        end
        s_r = su - s_l
        s2_r = su2 - s2_l
        nr = n - nl
        # This check is necessary I think because in theory all labels could
        # be the same, then either nl or nr would be 0. - @cstjean
        if nr > 0 && nl > 0
            loss = s2_l - s_l^2/nl + s2_r - s_r^2/nr
            if -loss > best_val
                best_val = -loss
                best_thresh = thresh
            end
        end
    end
    return best_val, best_thresh
end

function build_stump{T<:Float64, U<:Real}(labels::Vector{T}, features::Matrix{U}; rng=Base.GLOBAL_RNG)
    S = _split_mse(labels, features, 0, rng)
    if S == NO_BEST
        return Leaf(mean(labels), labels)
    end
    id, thresh = S
    split = features[:,id] .< thresh
    return Node(id, thresh,
                Leaf(mean(labels[split]), labels[split]),
                Leaf(mean(labels[neg(split)]), labels[neg(split)]))
end

function build_tree{T<:Float64, U<:Real}(labels::Vector{T}, features::Matrix{U}, maxlabels=5, nsubfeatures=0, maxdepth=-1; rng=Base.GLOBAL_RNG)
    if maxdepth < -1
        error("Unexpected value for maxdepth: $(maxdepth) (expected: maxdepth >= 0, or maxdepth = -1 for infinite depth)")
    end
    if length(labels) <= maxlabels || maxdepth==0
        return Leaf(mean(labels), labels)
    end
    S = _split_mse(labels, features, nsubfeatures, rng)
    if S == NO_BEST
        return Leaf(mean(labels), labels)
    end
    id, thresh = S
    split = features[:,id] .< thresh
    return Node(id, thresh,
                build_tree(labels[split], features[split,:], maxlabels, nsubfeatures, max(maxdepth-1, -1); rng=rng),
                build_tree(labels[neg(split)], features[neg(split),:], maxlabels, nsubfeatures, max(maxdepth-1, -1); rng=rng))
end

function build_forest{T<:Float64, U<:Real}(labels::Vector{T}, features::Matrix{U}, nsubfeatures::Integer, ntrees::Integer, maxlabels=5, partialsampling=0.7, maxdepth=-1; rng=Base.GLOBAL_RNG)
    rng = mk_rng(rng)::AbstractRNG
    partialsampling = partialsampling > 1.0 ? 1.0 : partialsampling
    Nlabels = length(labels)
    Nsamples = _int(partialsampling * Nlabels)
    forest = @parallel (vcat) for i in 1:ntrees
        inds = rand(rng, 1:Nlabels, Nsamples)
        build_tree(labels[inds], features[inds,:], maxlabels, nsubfeatures, maxdepth; rng=rng)
    end
    return Ensemble([forest;])
end


text:
## Exercise 6-3
## Write a function isbetween(x, y, z) that returns true if x ≤ y ≤ z or false otherwise.
println("Ans: ")

function isBetween(x, y, z)
    return x ≤ y ≤ z
end

println(isBetween(1, 2, 3))

println("End.")


text:
using BenchmarkTools
using Printf

"""
    evaluatefunctions(N)

Evaluates trigonometric functions on N values spaced
evenly between -1500.0 and 1500.0.
"""
function evaluatefunctions(N)
    #x = linspace(-1500.0, 1500.0, N)
    x = range(-1500.0, stop=1500.0, length=N)
    M = 10000
    for i in 1:M
        y = sin.(x)
        x = asin.(y)
        y = cos.(x)
        x = acos.(y)
        y = tan.(x)
        x = atan.(y)
    end
end

# Get the number of values from the command line.
N = parse(Int, ARGS[1])

println("-----------------------------------------")
println("Evaluate functions: ", N)
println("-----------------------------------------")
println(" ")

@btime evaluatefunctions(N)