text:
struct SumKernel{K1<:Kernel, K2<:Kernel} <: PairKernel{K1, K2}
    kleft::K1
    kright::K2
end
leftkern(sumkern::SumKernel) = sumkern.kleft
rightkern(sumkern::SumKernel) = sumkern.kright

get_param_names(sumkern::SumKernel) = composite_param_names(components(sumkern), :sk)

function Statistics.cov(sk::SumKernel, x::VecF64, y::VecF64)
    cov(sk.kleft, x, y) + cov(sk.kright, x, y)
end

@inline cov_ij(k::K, X::MatF64, i::Int, j::Int, dim::Int) where {K<:SumKernel} = cov_ij(k.kleft, X, i, j, dim) + cov_ij(k.kright, X, i, j, dim)
@inline cov_ij(k::K, X::MatF64, data::PairData, i::Int, j::Int, dim::Int) where {K<:SumKernel} = cov_ij(k.kleft, X, data.data1, i, j, dim) + cov_ij(k.kright, X, data.data2, i, j, dim)

@inline function dKij_dθp(sumkern::SumKernel, X::MatF64, i::Int, j::Int, p::Int, dim::Int)
    np = num_params(sumkern.kleft)
    if p<=np
        return dKij_dθp(sumkern.kleft, X, i, j, p, dim)
    else
        return dKij_dθp(sumkern.kright, X, i, j, p-np, dim)
    end
end
@inline function dKij_dθp(sumkern::SumKernel, X::MatF64, data::PairData, i::Int, j::Int, p::Int, dim::Int)
    np = num_params(sumkern.kleft)
    if p<=np
        return dKij_dθp(sumkern.kleft, X, data.data1, i, j, p, dim)
    else
        return dKij_dθp(sumkern.kright, X, data.data2, i, j, p-np, dim)
    end
end
@inline @inbounds function dKij_dθ!(dK::VecF64, sumkern::SumKernel, X::MatF64, data::PairData, i::Int, j::Int, dim::Int, npars::Int)
    npright = num_params(sumkern.kright)
    npleft = num_params(sumkern.kleft)
    dKij_dθ!(dK, sumkern.kright, X, data.data2, i, j, dim, npars-npleft)
    for ipar in npright:-1:1
        dK[npleft+ipar] = dK[ipar]
    end
    dKij_dθ!(dK,  sumkern.kleft,  X, data.data1, i, j, dim, npleft)
end

function grad_slice!(dK::MatF64, sumkern::SumKernel, X::MatF64, data::PairData, p::Int)
    np = num_params(sumkern.kleft)
    if p<=np
        return grad_slice!(dK, sumkern.kleft, X, data.data1, p)
    else
        return grad_slice!(dK, sumkern.kright, X, data.data2, p-np)
    end
end

# Addition operators
Base.:+(kleft::Kernel, kright::Kernel) = SumKernel(kleft,kright)


text:
## ------------------------------------------------------
# globals
const UNDONE_SIM_STATUS = :UNDONE
const DEAD_SIM_STATUS = :DEAD
const EXPLODED_SIM_STATUS = :EXPLODED
const NITERS_SIM_STATUS = :NITERS
const STST_SIM_STATUS = :STST
const STATUS_FILE_EXT = ".status.jls"
const BATCH_FILE_EXT = ".batch.jls"
const SIM_FILE_EXT = ".simdat.jls"
const BFILES_FILE_EXT = ".bfiles.jls"
const ME_FILE_EXT = ".maxent.jls"
const ME_ALGTEST_FILE_EXT = ".algtest.jls"

## ------------------------------------------------------
# utils
status_file(simid, simparams) = procdir(Dynamic, simid, simparams, STATUS_FILE_EXT)
set_status(status, simid, simparams) = (sdat(Dynamic, status, status_file(simid, simparams)); status)
get_status(simid, simparams) = ldat(() -> UNDONE_SIM_STATUS, status_file(simid, simparams))

## ------------------------------------------------------
# batch files
bfiles_file(simid, simparams) = procdir(Dynamic, simid, simparams, BFILES_FILE_EXT)
save_bfiles(bfiles, simid, simparams) = sdat(bfiles, bfiles_file(simid, simparams))
load_bfiles(simid, simparams) = ldat(bfiles_file(simid, simparams))

## ------------------------------------------------------
# batch
batch_file(simid, simparams, it) = procdir(Dynamic, simid, simparams, (;it), BATCH_FILE_EXT)
save_batch(dat, simid, simparams, it) =
   (bfile = batch_file(simid, simparams, it); sdat(dat, bfile); bfile)
load_batch(simid, simparams, it) = ldat(batch_file(simid, simparams, it))

## ------------------------------------------------------
# simdat
simdat_file(simid, simparams) = procdir(Dynamic, simid, simparams, SIM_FILE_EXT)
save_simdat(S, simid, simparams) = sdat(S, simdat_file(simid, simparams))
load_simdat(simid, simparams) = ldat(simdat_file(simid, simparams))

## ------------------------------------------------------
# maxent
maxent_file(simid, MEmode, simparams) =
   procdir(Dynamic, simid, simparams, (;MEmode = string(MEmode)), ME_FILE_EXT)
save_maxent(dat, simid, MEmode, simparams) =
   sdat(dat, maxent_file(simid, MEmode, simparams))
load_maxent(simid, MEmode, simparams) =
   ldat(maxent_file(simid, MEmode, simparams))

## ------------------------------------------------------
# maxent probe
algtest_file(simid, simparams...) =
   procdir(Dynamic, simid, simparams..., ME_ALGTEST_FILE_EXT)
save_algtest(dat, simid, simparams...) =
   sdat(dat, algtest_file(simid, simparams...))
load_algtest(simid, simparams...) =
   ldat(algtest_file(simid, simparams...))

## ------------------------------------------------------
function inparams(nparams::Dict; qparams...)
   qparams = Dict(qparams...)
   isempty(qparams) && return false
   isempty(nparams) && return false
   for (k, qparam) in qparams
      (get(nparams, string(k), nothing) != qparam) && return false
   end
   return true
end

## ------------------------------------------------------
function collect_ts(tsid, simid, qparams)
   dat_pool = Dict()
   bfile_params = load_bfiles(simid, qparams)
   t0 = 1
   for (simid, simparams, it) in bfile_params
      # batch
      fname = batch_file(simid, simparams, it)
      batch = ldat(fname)
      bs = getfield(batch, tsid)
      @assert (bs isa Vector)
      blen = length(bs)

      bt = (t0:(t0 + blen - 1)) .* batch.push_frec
      ts = get!(dat_pool, :ts, Float64[])
      push!(ts, bt...)
      bdat = identity.(bs)
      ds = get!(dat_pool, :ds, eltype(bdat)[])
      push!(ds, bdat...)
      t0 += length(bs) + 1
   end
   ts, ds = get(dat_pool, :ts, []), get(dat_pool, :ds, [])
   return (ts, ds)
end

## ------------------------------------------------------
function is_feasible_stst(net::MetNet, D, cgD_X; tol = 0.0)
   L, _ = fixxing(net, :z, D; tol = abs(D * tol)) do
      fva(net, :ug)
   end
   L < cgD_X
end

ext: jl


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 74207314-8c47-11eb-0989-f3dab6bcffd2
using PlutoUI

# ╔═╡ 2f26e06e-8c08-11eb-18c0-2f8e5b0cf47a
using CUDA

# ╔═╡ 590fdac2-8c08-11eb-1c42-93ff500817c1
using GPUArrays: allowscalar

# ╔═╡ 4068759c-8c08-11eb-1d2f-d79efe98c2b0
using Flux

# ╔═╡ 8297be82-8c08-11eb-39fb-eb37f75d09e4
using Flux: update!, reset!, loadparams!, onehot

# ╔═╡ 91e66b2c-8c08-11eb-1a3e-891e6401eb8b
using Transformers

# ╔═╡ 92baa6f8-8c08-11eb-3321-233e7bf14afe
using Transformers.Basic

# ╔═╡ 95f72c2e-8c08-11eb-2e18-0b166f648de6
using Transformers.Pretrain

# ╔═╡ 98c4351e-8c08-11eb-14fe-7df792e4b7cf
using BSON: @save, @load

# ╔═╡ 44df9b3a-8c16-11eb-17a7-3d6038911fc5
using Dates

# ╔═╡ 2ed19a74-8c18-11eb-0c06-95c930bbc8c1
using Statistics

# ╔═╡ 30edfac8-8c18-11eb-1d59-a5042a177b20
using Plots

# ╔═╡ 06f8595a-8c42-11eb-3b1c-011d0706850e
using PreSumm

# ╔═╡ 998b0178-8c60-11eb-2804-190b910dc7d8
begin
	# Workaround for issue:
	# https://github.com/fonsp/Pluto.jl/issues/301
	using Zygote, BSON
	Core.eval(Main, :(using Zygote))
end

# ╔═╡ 702d3820-84d9-11eb-1895-1d00242e5363
md"""
# Reproducing Text Summarization with Pretrained Encoders

This is a reproducability study on the ["Text Summarization with Pretrained Encoders"](https://doi.org/10.18653/v1/D19-1387) paper by Yang Liu and Mirella Lapata.
"""

# ╔═╡ fc420da6-8cfc-11eb-1e05-b71010313c62
TableOfContents()

# ╔═╡ 176ccb48-8c08-11eb-3068-3b684ff378b5
md"""
## Setup
"""

# ╔═╡ 400c92b8-8c17-11eb-170d-1f6b3fe06a43
md"""
Should we train the model before evaluating?

Yes $(@bind TRAIN CheckBox(default=(get(ENV, "TRAIN", "false") == "true")))

Keep in mind that training requires a lot of hardware resources and takes a long while! See the _Training loop_ section for more details.
"""

# ╔═╡ 19fec164-8c43-11eb-3e55-1dc0c30bf8c1
md"""
Should we evaluate the trained model?

Yes $(@bind EVALUATE CheckBox(default=(get(ENV, "EVALUATE", "true") == "true")))
"""

# ╔═╡ 6847039e-8c17-11eb-15b9-e70a268a56fb
md"""
Should we use a very simple model for less powerfull machines?

Yes $(@bind DEBUG CheckBox(default=(get(ENV, "DEBUG", "true") == "true")))

If this checkbox is ticked, we'll drastically over-simplify the model and data so that you should be able to train a tiny variant of the actual model on your own machine. See the _Model_ section for more details.
"""

# ╔═╡ 3351b1ca-8c08-11eb-14c3-8f61900721f4
md"""
### Load packages
"""

# ╔═╡ 1be31330-8c08-11eb-296f-6f5df8bf6e41
md"""
### Setup GPU with CUDA
"""

# ╔═╡ 2646eafe-8c08-11eb-2a25-c338673a3d2a
if CUDA.functional(true)
    # CUDA.device!(1)
    allowscalar(false)
    Flux.use_cuda[] = true
    enable_gpu(true)
else
    Flux.use_cuda[] = false
    enable_gpu(false)
end;

# ╔═╡ dbc654f4-8c18-11eb-059e-f91749357883
if CUDA.functional(true)
    md"Alright, you're running in GPU-accellerated mode using CUDA."
else
    md"""
	**⚠ You're running the notebook without GPU support.**

	It is likely that you won't be able to train our model.
	At least try setting `DEBUG = true`, better yet skip training with `TRAIN = false` and use pretrained model snapshots
	"""
end

# ╔═╡ 39e1ae40-8c18-11eb-3f12-b566f7f516e2
md"""
### Setup plots
"""

# ╔═╡ 48c4e332-8c18-11eb-05c8-f15576c14fcb
gr()

# ╔═╡ 493f6f90-8c26-11eb-0627-fb59af397184
theme(:vibrant)

# ╔═╡ 22e29d1a-84e2-11eb-3e35-f7050cfd6cc4
md"""
## Datasets
We're loading datasets with DataDeps.jl.
This way each dataset is cached and we only have to download once.
If you haven't previously downloaded the datasets, _this process will take a while_.
"""

# ╔═╡ 0d86a904-84e3-11eb-1c0c-9d37637c8420
ENV["DATADEPS_ALWAYS_ACCEPT"] = true; # Don't ask for downloading.

# ╔═╡ 29348670-84e4-11eb-076b-b78a92e94642
ENV["DATADEPS_PROGRESS_UPDATE_PERIOD"] = 5; # Log process to the console, though.

# ╔═╡ 98474e70-8c06-11eb-28c4-35bd54c8a558
md"""
### Data loading
"""

# ╔═╡ 44848686-8c0b-11eb-064c-7f56c697e69c
md"""
#### Training data
Load article-summary pairs from the CNN / Daily Mail dataset's train split.
We can get a fresh `Channel` (iterator) with training data each time we call `cnndm_train()`.
"""

# ╔═╡ 9109e4e8-8c46-11eb-296d-154cfbd4db73
data(split) = cnn_dm(split)

# ╔═╡ 32af5d7e-8c0a-11eb-1784-17984db6e6fc
data_train() = data("train")

# ╔═╡ 0d56f78e-8c0b-11eb-3d74-218c28817970
first(data_train())

# ╔═╡ d4b77c36-8c0b-11eb-3a32-b3fc1b7e80af
md"""
#### Validation data
Get a `Channel` of article-summary pairs from the dataset's development/validation split.
"""

# ╔═╡ 2a4a03e0-8c0b-11eb-0bea-2d96d67712fc
data_dev() = data("valid")

# ╔═╡ 3a4f4ce4-8c0b-11eb-3518-0907e3d2a62d
first(data_dev())

# ╔═╡ 95b154c8-8c0b-11eb-293f-f15c59a529e8
md"""
#### Test data
Get a `Channel` of article-summary pairs from the dataset's test split.
"""

# ╔═╡ 1c697a62-8c0b-11eb-0814-0d02d802f6fc
data_test() = data("test")

# ╔═╡ 395e445e-8c0b-11eb-0193-17ea21763de6
first(data_test())

# ╔═╡ 6462ab2e-84e3-11eb-33dd-13e5000c78da
md"Load preprocessed CNN/Dailymail data."

# ╔═╡ 58544bfa-84d9-11eb-3426-b3b5cc2a19a8
md"""
### Preprocessing
"""

# ╔═╡ a4f8ec0a-8c06-11eb-042e-290b3405e5e6
md"""
## Pretrained BERT model
Load the pretrained BERT model from Google using Transformers.jl.
We load the uncased BERT-base variant with 12 layers of hidden size 768.
The model is also loaded internally with DataDeps.jl.
"""

# ╔═╡ 0303fb1e-8c0c-11eb-28d1-1d9e631cb8af
bert_model, wordpiece, tokenizer = pretrain"bert-uncased_L-12_H-768_A-12"

# ╔═╡ ff141dee-8c0c-11eb-35e3-c7a7fe623dc3
md"""
### Vocabulary
When debugging, to keep the model trainable on smaller machines, limit the vocabulary to 100 words.
"""

# ╔═╡ 58ae6888-8c0c-11eb-359c-ed1a2cafbfde
function load_vocabulary()::Vocabulary
	if !DEBUG
		Vocabulary(wordpiece) |> gpu
	else
		Vocabulary(wordpiece.vocab[1:100], wordpiece.vocab[100]) |> gpu
	end
end

# ╔═╡ b1ac4c9e-8c0d-11eb-035a-4b075ffeb8f5
vocabulary = load_vocabulary()

# ╔═╡ fc930278-8c0d-11eb-1c70-2176689ce05c
md"""
## Model
"""

# ╔═╡ aab3df80-8c0e-11eb-37bd-bd40e5a498c7
md"""
When debugging on smaller machines, use a tiny transformer encoder-decoder variant instead of the large BERT-based model.
"""

# ╔═╡ 0cb47e48-8c0e-11eb-0a4d-8f1bd42be2b8
function load_model()::Translator
	if !DEBUG
		BertAbs(bert_model, length(vocabulary)) |> gpu
	else
		TransformerAbsTiny(length(vocabulary)) |> gpu
	end
end

# ╔═╡ e6db6ba4-8c0e-11eb-1bdb-3bc0b02e2b4f
model = load_model()

# ╔═╡ 3f4853dc-8c06-11eb-3a59-cdb9fc854879
md"""
## Training
"""

# ╔═╡ 29f82622-8c0f-11eb-211b-1d275e80ebbf
md"""
### Preprocessing
Preprocess a `String` to a sequence of tokens (not necessarily words), tokenized by the BERT tokenizer.
"""

# ╔═╡ 0326189c-8c0f-11eb-38d1-595283980478
function preprocess(text::String)::AbstractVector{String}
    tokens = text |> tokenizer |> wordpiece
    max_length = min(4096, length(tokens)) # Truncate to 4096 tokens
    tokens = tokens[1:max_length]
    return ["[CLS]"; tokens; "[SEP]"]
end

# ╔═╡ e1a94230-8c10-11eb-3a01-e765b58c2c29
md"""
### Prediction & loss calculation
"""

# ╔═╡ 5144a378-8c0f-11eb-3b32-6f7661781a0a
md"""
Label smoothing factor $\alpha$ to apply to the one-hot ground-truth labels when calculating the model loss.
"""

# ╔═╡ 619a2256-8c08-11eb-3544-2b20fba8a71d
# FIXME Label smoothing doesn't work yet.
label_smoothing_α = 0.0

# ╔═╡ c8cf0d90-8c0f-11eb-3a11-4bded8f5be4f
md"""
Calculate the model loss for given input, output, and ground truth token probabilities.
If $\alpha > 0$, return Kullback-Leibler divergence between the model's predictions and smoothed ground truth labels. Otherwise cross entropy between predictions and unsmoothed ground truth.
"""

# ╔═╡ 811777ce-8c0f-11eb-2ac4-b51b85cd6bc9
function loss(
    inputs::AbstractVector{String},
    outputs::AbstractVector{String},
    ground_truth::AbstractMatrix{<:Number},
    translator::Translator
)::AbstractFloat
    prediction = translator.transformers(vocabulary, inputs, outputs)
    loss = logtranslationloss(prediction, ground_truth, α=label_smoothing_α)
    return loss
end

# ╔═╡ a46bcb12-8c11-11eb-0493-1505139bb2bf
md"""
### Model parameters
Parameters are separated for the encoder and the rest of the model (i.e., embeddings, decoder, and generator).
This accounts for the former being pretrained while the latter have to be learned from scratch.
"""

# ╔═╡ 390e1432-8c12-11eb-2697-13d8f3c0e2bf
parameters_encoder = params(model.transformers.encoder)

# ╔═╡ 3ddf05d4-8c12-11eb-0a80-e50ef061dc33
parameters_decoder = params(
    model.transformers.embed,
    model.transformers.decoder,
    model.transformers.generator
)

# ╔═╡ ca3a1fe4-8c3c-11eb-243e-b58cc96b2d1d
md"""
We're going to tune
$(params_count(parameters_encoder))
parameters from the pretrained encoder and train
$(params_count(parameters_decoder))
parameters for the non-pretrained model parts.
(That is
$(params_count(params(model.transformers.embed)))
parameters from embeddings,
$(params_count(params(model.transformers.decoder)))
from the decoder, and
$(params_count(params(model.transformers.generator)))
from the generator.)
"""

# ╔═╡ bf061ff0-8c10-11eb-108d-cf008ae94342
md"""
### Optimizers
We use ADAM optimizers with a custom warmup schedule for the learning rate $\eta$.
The model's encoder is optimized slower, because it has already bee pretrained.
"""

# ╔═╡ 486cc03c-8c11-11eb-03cb-bb640ac822f5
md"""
Encoder schedule:
$\eta = 2\mathcal{e}^{-3} \cdot \min(\text{step}^{-0.5},\ \text{step} \cdot 20\,000)$
"""

# ╔═╡ 2c0590d8-8c11-11eb-1ac0-41a24f84ed58
optimizer_encoder = WarmupADAM(η=2ℯ^(-3), w=20_000, β=(0.9, 0.99)) |> gpu

# ╔═╡ 940c2c58-8c11-11eb-39b6-cf4dc61bd550
md"""
Decoder schedule:
$\eta = 0.1 \cdot \min(\text{step}^{-0.5},\ \text{step} \cdot 10\,000)$
"""

# ╔═╡ 1cd3487e-8c11-11eb-173a-af13bddc850f
optimizer_decoder = WarmupADAM(η=0.1, w=10_000, β=(0.9, 0.99)) |> gpu

# ╔═╡ 03a78c44-8c13-11eb-27da-e11097968bc7
max_steps = !DEBUG ? 200_000 : 200

# ╔═╡ d32d3994-8c12-11eb-1848-d9bd515709a5
md"""
### Training loop
Here we define the training loop that is iterated over for at most $max_steps steps.
"""

# ╔═╡ 07e53e8e-8c13-11eb-1188-5309e353c56a
snapshot_steps = !DEBUG ? 2500 : 10

# ╔═╡ 79b63cd0-8c1c-11eb-3b31-1395f64e979d
start_time = now()

# ╔═╡ ef68c91e-8c3c-11eb-393a-5145b8cff911
md"""
The actual training loop extracts tokens for each row's source and target text (original article vs. short summary).
Then the loss between the model's predicted token probabilities and the ground truth probability is compared.
The model is then updated with gradients for the encoder and gradients for decoder, embeddings, and generator.

Losses are saved from every step and snapshots of the model weights and losses/optimizers for both parameter sets are saved in BSON format to the output folder:
$(out_dir())
"""

# ╔═╡ ae61b104-8c57-11eb-3074-d5445b61f010
train_loss(inputs, outputs, ground_truth) = loss(inputs, outputs, ground_truth, model)

# ╔═╡ b4dbc124-8c15-11eb-009d-7b3e0a26ba4d
function save_snapshot(
		time::DateTime,
		step::Int,
		model::Translator,
		optimizer_encoder::WarmupADAM,
		optimizer_decoder::WarmupADAM,
		losses_encoder::AbstractVector{<:AbstractFloat},
		losses_decoder::AbstractVector{<:AbstractFloat})
	@info "Save model snapshot."
	weights = params(model)
	snapfile(time, step, name) = snapshot_file(time, step, name)
	@save snapfile(time, step, "model.bson") model
	@save snapfile(time, step, "weights.bson") weights
	@save snapfile(time, step, "optimizer-encoder.bson") optimizer_encoder
	@save snapfile(time, step, "optimizer-decoder.bson") optimizer_decoder
	@save snapfile(time, step, "losses-encoder.bson") losses_encoder
	@save snapfile(time, step, "losses-decoder.bson") losses_decoder
end

# ╔═╡ 3f117812-8c13-11eb-3b34-3988290737a7
function train!(model::Translator)
	losses_encoder = []
	losses_decoder = []

	data = data_train()
	for (step, summary) ∈ zip(1:max_steps, data)
   		@info "Training step $step/$max_steps."
		inputs = summary.source |> preprocess |> gpu
		outputs = summary.target |> preprocess |> gpu
		ground_truth = onehot(vocabulary, outputs) |> gpu

		@info "Calculate encoder gradients."
		local loss_encoder
		gradients_encoder = gradient(parameters_encoder) do
			loss_encoder = train_loss(inputs, outputs, ground_truth)
			return loss_encoder
		end
		push!(losses_encoder, loss_encoder)

		@info "Calculate decoder, embeddings, and generator gradients."
		local loss_decoder
		gradients_decoder = gradient(parameters_decoder) do
			loss_decoder = train_loss(inputs, outputs, ground_truth)
			return loss_decoder
		end
		push!(losses_decoder, loss_decoder)

		@info "Updating encoder parameters." loss_encoder
		update!(optimizer_encoder, parameters_encoder, gradients_encoder)

		@info "Updating decoder, embeddings, and generator parameters." loss_decoder
		update!(optimizer_decoder, parameters_decoder, gradients_decoder)

		if step % snapshot_steps == 0
			save_snapshot(
				start_time,
				step,
				model |> cpu,
				optimizer_encoder |> cpu,
				optimizer_decoder |> cpu,
				losses_encoder |> cpu,
				losses_decoder |> cpu,
			)
		end
	end
end

# ╔═╡ be5c7b3c-8c16-11eb-1a3f-01eadf6bc1bb
md"""
Now train the model. Lean back and take a sip of coffee, as this takes a long while. ☕

_Tip: If you've already trained the model, untick the training checkbox at the beginning of the notebook and the evaluation will continue without training the model first._
"""

# ╔═╡ ded97e7a-8c16-11eb-2b04-23a5dec4405e
if TRAIN
	train!(model)
end

# ╔═╡ 4e3ead50-8c06-11eb-39ff-a3834ba819c2
md"""
## Evaluation
"""

# ╔═╡ 566a15e4-8c44-11eb-3d27-a1b28fc3057d
if EVALUATE
	md"""
	Evaluate training loss, select best model by development/validation loss,
	and measure summary quality.
	"""
else
	md"""
	**⚠ Skipping evaluation is not yet fully supported.**
	"""
end

# ╔═╡ 35878a0a-8c1c-11eb-0e18-d1964db60a67
md"""
### Configuration
Define the training run we want to evaluate.
If training is activated, this is automatically set to the current training run.
Otherwise, you can select a file manually.
"""

# ╔═╡ 9bcff242-8c1d-11eb-007f-6b3ae5b2ef6a
md"""
When did the chosen model start training?
"""

# ╔═╡ 546c963e-8c1c-11eb-0762-cbecc179fdf0
# custom_start_time = DateTime(2021, 03, 24, 06, 06)
# custom_start_time = DateTime(2021, 03, 22, 17, 50)
custom_start_time = DateTime(2021, 03, 24, 12, 05)

# ╔═╡ 655f933a-8c1c-11eb-1fb0-636b021637b6
eval_start_time = TRAIN ? start_time : custom_start_time

# ╔═╡ 16441198-8c1e-11eb-05af-99619660cdd8
md"""
What was the rate at which snapshots were saved for that model?
"""

# ╔═╡ 87ceb29c-8c1d-11eb-2534-39649a087b4f
# custom_snapshot_steps = 10
custom_snapshot_steps = 2500

# ╔═╡ bd5cde20-8c1d-11eb-25f4-c1d10590778a
eval_snapshot_steps = TRAIN ? snapshot_steps : custom_snapshot_steps

# ╔═╡ 110108b4-8c3d-11eb-230e-f1f8ef3985ff
md"""
### Model loading
At this point we've either trained the model from the previous sections or copied pretrained training snapshots to the output folder:
$(out_dir())
"""

# ╔═╡ dc238b06-8c1d-11eb-3cf2-495db70ac90a
weights_snapshots, max_step = EVALUATE ? snapshot_files(
	eval_start_time,
	eval_snapshot_steps,
	"weights.bson"
) : ([],0)

# ╔═╡ c44ed108-8c49-11eb-30d2-73e41b46aed7
has_weights_snapshot = !isempty(weights_snapshots)

# ╔═╡ 873bb7ca-8c1e-11eb-08f7-755484003009
if !has_weights_snapshot
	md"""
	**⚠ Could not find any weights snapshot for training
	started at $(Dates.format(eval_start_time, "yyyy-mm-dd (at HH:MM)"))
	with snapshotting every $eval_snapshot_steps steps.**
	"""
else
	md"""
	Found $(length(weights_snapshots)) snapshots for training
	started at $(Dates.format(eval_start_time, "yyyy-mm-dd (at HH:MM)"))
	with snapshotting every $eval_snapshot_steps steps.
	The last snapshot is from step $max_step.
	"""
end

# ╔═╡ 9fcbdd2e-8c21-11eb-0e59-63649715d380
eval_snapshot_file(name) = snapshot_file(eval_start_time, max_step, name)

# ╔═╡ 216e718c-8c1c-11eb-066f-f705805cb70e
md"""
### Training loss
Show loss on the training set at each iteration
"""

# ╔═╡ 7e4da548-8c22-11eb-1b3c-81269c4f5f5a
function plot_train_loss()
	if !has_weights_snapshot
		return
	end

	@load eval_snapshot_file("losses-encoder.bson") losses_encoder
	@load eval_snapshot_file("losses-decoder.bson") losses_decoder
	loss_plot = plot(
		title = "Generator classification loss during training",
		xlabel = "steps",
		ylabel = "cross entropy loss",
    	# xticks = 0:eval_snapshot_steps:max_step,
	)
	plot!(loss_plot, losses_encoder[1:100], label="encoder")
	plot!(loss_plot, losses_decoder[1:100], label="decoder, embeddings, generator")
	savefig(loss_plot, joinpath(out_dir(), "training-loss.pdf"))
	savefig(loss_plot, joinpath(out_dir(), "training-loss.png"))
	loss_plot
end;

# ╔═╡ a39ba53e-8c22-11eb-04af-95d7d112f67d
plot_train_loss()

# ╔═╡ 16ccaf1a-8c2a-11eb-274a-5d8831732be1
md"""
### Best model selection
Next, we'll select the best model by classification loss on the development/validation set.
"""

# ╔═╡ b3692f86-8c1b-11eb-3425-35027febcec0
md"""
#### Development data
For evaluating loss on the development/validation set, we load its summary pairs.
Loading and caching to the GPU might take a short moment.
"""

# ╔═╡ 97423bde-8c18-11eb-2b01-8526b32161d4
dev_summaries = collect(data_dev()) |> gpu

# ╔═╡ fcd40322-8c6c-11eb-0259-1d2abf1c7b03
md"""
Select the best model by mean loss on validation data.
"""

# ╔═╡ 40f5d3de-8c2a-11eb-2e1b-d39b430a6550
function dev_loss(translator::Translator; agg=mean)::AbstractFloat
    losses = []
	summaries = !DEBUG ? dev_summaries : dev_summaries[1:10]
    for summary ∈ summaries
        inputs = summary.source |> preprocess |> gpu
        outputs = summary.target |> preprocess |> gpu
        ground_truth = onehot(vocabulary, outputs) |> gpu
        current_loss = loss(inputs, outputs, ground_truth, translator)
        push!(losses, current_loss)
    end
    return agg(losses)
end

# ╔═╡ 1fc901de-8c6d-11eb-0c4e-394aa014313f
md"""
Load weights of each snapshot, load weights into model and check if the model performs better on the dev set.
"""

# ╔═╡ 7cb9dfd2-8c2a-11eb-27ac-811c01ba9368
function find_best_model(;agg=mean)::Translator
    best_loss = Inf
	best_weights = nothing
    model = load_model()
    for snapshot ∈ weights_snapshots
        @info "Load model snapshot $snapshot."
		@load snapshot weights
		loadparams!(model, weights)
        loss = dev_loss(model, agg=agg)
        if loss < best_loss
			@info "New best model."
            best_loss = loss
			best_weights = weights
        end
    end
	if best_weights == nothing
		throw(ErrorException("Could not find model."))
	end
	loadparams!(model, best_weights)
    return model
end

# ╔═╡ c33a3092-8c2a-11eb-03c3-d74e98a3867c
best_model = find_best_model()

# ╔═╡ e2449708-8c06-11eb-0d09-991e4917b9d3
md"""
### Automatic evaluation
"""

# ╔═╡ 26d83e30-8c28-11eb-3107-a3df1788bc2f
md"""
#### Test data
For evaluating quality on the test set, we load its summary pairs.
Loading and caching to the GPU might take a short moment.
"""

# ╔═╡ 3e25c08a-8c28-11eb-143b-b996590c22d4
test_summaries = collect(data_test()) |> gpu

# ╔═╡ a5adebd0-8c6c-11eb-14af-019f37d78469
test_source = test_summaries[1].source |> preprocess |> gpu

# ╔═╡ aae5f0f0-8d02-11eb-2397-93f70ee24eb1
model(test_source, vocabulary)

# ╔═╡ ab6b0af8-8c62-11eb-02ad-d94033352713
best_model(test_source, vocabulary)

# ╔═╡ a2ef84f2-8c07-11eb-0a94-21ebe1ed471a
md"""
#### ROUGE benchmark
Measure ROUGE-1, ROUGE-2, and ROUGE-L on the summaries generated for the test dataset.
"""

# ╔═╡ f67d71ea-8c06-11eb-06a9-550a8d54c139
md"""
### Manual evaluation
We generate summaries for the first 100 articles from the CNN / Daily Mail dataset's test split.
These summaries are maually examined and scored on the following scale:
- 0 = unreadable summary or unrelated to original article
- 1 = readable and related to original article, some redundancy
- 2 = captures exactly the article's main points, no redundancy
"""

# ╔═╡ Cell order:
# ╟─702d3820-84d9-11eb-1895-1d00242e5363
# ╠═fc420da6-8cfc-11eb-1e05-b71010313c62
# ╟─176ccb48-8c08-11eb-3068-3b684ff378b5
# ╟─400c92b8-8c17-11eb-170d-1f6b3fe06a43
# ╟─19fec164-8c43-11eb-3e55-1dc0c30bf8c1
# ╟─6847039e-8c17-11eb-15b9-e70a268a56fb
# ╟─3351b1ca-8c08-11eb-14c3-8f61900721f4
# ╠═74207314-8c47-11eb-0989-f3dab6bcffd2
# ╠═2f26e06e-8c08-11eb-18c0-2f8e5b0cf47a
# ╠═590fdac2-8c08-11eb-1c42-93ff500817c1
# ╠═4068759c-8c08-11eb-1d2f-d79efe98c2b0
# ╠═8297be82-8c08-11eb-39fb-eb37f75d09e4
# ╠═91e66b2c-8c08-11eb-1a3e-891e6401eb8b
# ╠═92baa6f8-8c08-11eb-3321-233e7bf14afe
# ╠═95f72c2e-8c08-11eb-2e18-0b166f648de6
# ╠═98c4351e-8c08-11eb-14fe-7df792e4b7cf
# ╠═44df9b3a-8c16-11eb-17a7-3d6038911fc5
# ╠═2ed19a74-8c18-11eb-0c06-95c930bbc8c1
# ╠═30edfac8-8c18-11eb-1d59-a5042a177b20
# ╠═06f8595a-8c42-11eb-3b1c-011d0706850e
# ╠═998b0178-8c60-11eb-2804-190b910dc7d8
# ╟─1be31330-8c08-11eb-296f-6f5df8bf6e41
# ╠═2646eafe-8c08-11eb-2a25-c338673a3d2a
# ╟─dbc654f4-8c18-11eb-059e-f91749357883
# ╟─39e1ae40-8c18-11eb-3f12-b566f7f516e2
# ╠═48c4e332-8c18-11eb-05c8-f15576c14fcb
# ╠═493f6f90-8c26-11eb-0627-fb59af397184
# ╟─22e29d1a-84e2-11eb-3e35-f7050cfd6cc4
# ╠═0d86a904-84e3-11eb-1c0c-9d37637c8420
# ╠═29348670-84e4-11eb-076b-b78a92e94642
# ╟─98474e70-8c06-11eb-28c4-35bd54c8a558
# ╟─44848686-8c0b-11eb-064c-7f56c697e69c
# ╠═9109e4e8-8c46-11eb-296d-154cfbd4db73
# ╠═32af5d7e-8c0a-11eb-1784-17984db6e6fc
# ╠═0d56f78e-8c0b-11eb-3d74-218c28817970
# ╟─d4b77c36-8c0b-11eb-3a32-b3fc1b7e80af
# ╠═2a4a03e0-8c0b-11eb-0bea-2d96d67712fc
# ╠═3a4f4ce4-8c0b-11eb-3518-0907e3d2a62d
# ╟─95b154c8-8c0b-11eb-293f-f15c59a529e8
# ╠═1c697a62-8c0b-11eb-0814-0d02d802f6fc
# ╠═395e445e-8c0b-11eb-0193-17ea21763de6
# ╟─6462ab2e-84e3-11eb-33dd-13e5000c78da
# ╟─58544bfa-84d9-11eb-3426-b3b5cc2a19a8
# ╟─a4f8ec0a-8c06-11eb-042e-290b3405e5e6
# ╠═0303fb1e-8c0c-11eb-28d1-1d9e631cb8af
# ╟─ff141dee-8c0c-11eb-35e3-c7a7fe623dc3
# ╠═58ae6888-8c0c-11eb-359c-ed1a2cafbfde
# ╠═b1ac4c9e-8c0d-11eb-035a-4b075ffeb8f5
# ╟─fc930278-8c0d-11eb-1c70-2176689ce05c
# ╟─aab3df80-8c0e-11eb-37bd-bd40e5a498c7
# ╠═0cb47e48-8c0e-11eb-0a4d-8f1bd42be2b8
# ╠═e6db6ba4-8c0e-11eb-1bdb-3bc0b02e2b4f
# ╟─3f4853dc-8c06-11eb-3a59-cdb9fc854879
# ╟─29f82622-8c0f-11eb-211b-1d275e80ebbf
# ╠═0326189c-8c0f-11eb-38d1-595283980478
# ╟─e1a94230-8c10-11eb-3a01-e765b58c2c29
# ╟─5144a378-8c0f-11eb-3b32-6f7661781a0a
# ╠═619a2256-8c08-11eb-3544-2b20fba8a71d
# ╟─c8cf0d90-8c0f-11eb-3a11-4bded8f5be4f
# ╠═811777ce-8c0f-11eb-2ac4-b51b85cd6bc9
# ╟─a46bcb12-8c11-11eb-0493-1505139bb2bf
# ╠═390e1432-8c12-11eb-2697-13d8f3c0e2bf
# ╠═3ddf05d4-8c12-11eb-0a80-e50ef061dc33
# ╟─ca3a1fe4-8c3c-11eb-243e-b58cc96b2d1d
# ╟─bf061ff0-8c10-11eb-108d-cf008ae94342
# ╟─486cc03c-8c11-11eb-03cb-bb640ac822f5
# ╠═2c0590d8-8c11-11eb-1ac0-41a24f84ed58
# ╟─940c2c58-8c11-11eb-39b6-cf4dc61bd550
# ╠═1cd3487e-8c11-11eb-173a-af13bddc850f
# ╟─d32d3994-8c12-11eb-1848-d9bd515709a5
# ╠═03a78c44-8c13-11eb-27da-e11097968bc7
# ╠═07e53e8e-8c13-11eb-1188-5309e353c56a
# ╠═79b63cd0-8c1c-11eb-3b31-1395f64e979d
# ╟─ef68c91e-8c3c-11eb-393a-5145b8cff911
# ╠═ae61b104-8c57-11eb-3074-d5445b61f010
# ╠═3f117812-8c13-11eb-3b34-3988290737a7
# ╠═b4dbc124-8c15-11eb-009d-7b3e0a26ba4d
# ╟─be5c7b3c-8c16-11eb-1a3f-01eadf6bc1bb
# ╠═ded97e7a-8c16-11eb-2b04-23a5dec4405e
# ╟─4e3ead50-8c06-11eb-39ff-a3834ba819c2
# ╟─566a15e4-8c44-11eb-3d27-a1b28fc3057d
# ╟─35878a0a-8c1c-11eb-0e18-d1964db60a67
# ╟─9bcff242-8c1d-11eb-007f-6b3ae5b2ef6a
# ╠═546c963e-8c1c-11eb-0762-cbecc179fdf0
# ╠═655f933a-8c1c-11eb-1fb0-636b021637b6
# ╟─16441198-8c1e-11eb-05af-99619660cdd8
# ╠═87ceb29c-8c1d-11eb-2534-39649a087b4f
# ╠═bd5cde20-8c1d-11eb-25f4-c1d10590778a
# ╟─110108b4-8c3d-11eb-230e-f1f8ef3985ff
# ╠═dc238b06-8c1d-11eb-3cf2-495db70ac90a
# ╠═c44ed108-8c49-11eb-30d2-73e41b46aed7
# ╟─873bb7ca-8c1e-11eb-08f7-755484003009
# ╠═9fcbdd2e-8c21-11eb-0e59-63649715d380
# ╟─216e718c-8c1c-11eb-066f-f705805cb70e
# ╠═7e4da548-8c22-11eb-1b3c-81269c4f5f5a
# ╠═a39ba53e-8c22-11eb-04af-95d7d112f67d
# ╟─16ccaf1a-8c2a-11eb-274a-5d8831732be1
# ╟─b3692f86-8c1b-11eb-3425-35027febcec0
# ╠═97423bde-8c18-11eb-2b01-8526b32161d4
# ╟─fcd40322-8c6c-11eb-0259-1d2abf1c7b03
# ╠═40f5d3de-8c2a-11eb-2e1b-d39b430a6550
# ╟─1fc901de-8c6d-11eb-0c4e-394aa014313f
# ╠═7cb9dfd2-8c2a-11eb-27ac-811c01ba9368
# ╠═c33a3092-8c2a-11eb-03c3-d74e98a3867c
# ╟─e2449708-8c06-11eb-0d09-991e4917b9d3
# ╟─26d83e30-8c28-11eb-3107-a3df1788bc2f
# ╠═3e25c08a-8c28-11eb-143b-b996590c22d4
# ╠═a5adebd0-8c6c-11eb-14af-019f37d78469
# ╠═aae5f0f0-8d02-11eb-2397-93f70ee24eb1
# ╠═ab6b0af8-8c62-11eb-02ad-d94033352713
# ╟─a2ef84f2-8c07-11eb-0a94-21ebe1ed471a
# ╟─f67d71ea-8c06-11eb-06a9-550a8d54c139


text:
## Bubble pressure solver
function x0_bubble_pressure(model::EoSModel,T,x)
    #TODO
    #on sufficiently large temps,
    #the joule-thompson inversion occurs
    #making the virial coeff positive
    #on those cases, use an strategy that supposes pure gas on that side
    #Pbi = inf
    #xi = 0

    #check each T with T_scale, if treshold is over, replace Pi with inf
    pure = split_model(model)
    crit = crit_pure.(pure)

    T_c = [tup[1] for tup in crit]
    V_c = [tup[3] for tup in crit]
    _0 = zero(T+first(x))
    nan = _0/_0
    sat_nan = (nan,nan,nan)
    replaceP = ifelse.(T_c .< T,true,false)

    eachx = eachcol(Diagonal(ones(eltype(x),length(x))))
#     Bi = second_virial_coefficient.(model,T,eachx)
    #using P_B(2B) as a sat aproximation
    #z = 1 + B/v
    #P_B = RT/v(1+B/v)
    #P_B(2B) = -RT/2B(1-B/2B)
    #P_B(2B) = -0.25*RT/B
    sat = [if !replaceP[i] sat_pure(pure[i],T) else sat_nan end for i in 1:length(pure)]

    P_sat = [tup[1] for tup in sat]
    V_l_sat = [tup[2] for tup in sat]
    V_v_sat = [tup[3] for tup in sat]
#     P_Bi = @. -0.25*R̄*T/Bi
    #=xP0 = yP
    #dot(x,P0) = P
    P = dot(x,P0)
    =#
    P = zero(T)
    V0_l = zero(T)
    V0_v = zero(T)
    Pi   = zero(x)
    for i in 1:length(x)
        if !replaceP[i]
            Pi[i] = P_sat[i][1]
            P+=x[i]*Pi[i]
            V0_l += x[i]*V_l_sat[i]
        else
            Pi[i] = pressure(pure[i],V_c[i],T)
            P+=x[i]*Pi[i]
            V0_l += x[i]*V_c[i]
        end
    end
    #@show P_Bi
    #P = dot(x,P_Bi)
    y = @. x*Pi/P
    ysum = 1/∑(y)
    y    = y.*ysum

    for i in 1:length(x)
        if !replaceP[i]
            V0_v += y[i]*V_v_sat[i]
        else
            V0_v += y[i]*V_c[i]*1.2
        end
    end

    prepend!(y,log10.([V0_l,V0_v]))
    return y
end

function bubble_pressure(model::EoSModel, T, x; v0 =nothing)
    TYPE = promote_type(eltype(T),eltype(x))
#     lb_v = lb_volume(model,x)
    ts = T_scales(model,x)
    pmix = p_scale(model,x)
    if v0 === nothing
        v0 = x0_bubble_pressure(model,T,x)
    end
    len = length(v0[1:end-1])
    #xcache = zeros(eltype(x0),len)
    Fcache = zeros(eltype(v0[1:end-1]),len)
    f! = (F,z) -> Obj_bubble_pressure(model, F, T, exp10(z[1]), exp10(z[2]), x,z[3:end],ts,pmix)
    r  =Solvers.nlsolve(f!,v0[1:end-1],LineSearch(Newton()))
    sol = Solvers.x_sol(r)
    v_l = exp10(sol[1])
    v_v = exp10(sol[2])
    y = FractionVector(sol[3:end])
    P_sat = pressure(model,v_l,T,x)
    return (P_sat, v_l, v_v, y)
end

function Obj_bubble_pressure(model::EoSModel, F, T, v_l, v_v, x, y,ts,ps)
    y   = FractionVector(y) #julia magic, check misc.jl
    μ_l = VT_chemical_potential(model,v_l,T,x)
    μ_v = VT_chemical_potential(model,v_v,T,y)
    p_l = pressure(model,v_l,T,x)
    p_v = pressure(model,v_v,T,y)
    for i in 1:length(x)
        F[i] = (μ_l[i]-μ_v[i])/(R̄*ts[i])
    end
    F[end] = (p_l-p_v)/ps
    return F
end

    #j! = (J,z) -> Jac_bubble_pressure(model, J, T, exp10(z[1]), exp10(z[2]), x[i,:], z[3:end])
    #=
    _y0 = collect(FractionVector(v0[3:end]))
    @show _y0
    @show _pl0 = pressure(model,exp10(v0[1]),T,x)
    @show _pv0 = pressure(model,exp10(v0[2]),T,_y0)
    _P = 0.5*(_pl0+_pv0)
    for _ in 1:3
        _y0,_P = rr_bubble_pressure_refine(model,x,_y0,_P,T)
    end
    =#
# function bubble_pressure_rr(model, T, x; P = 40000)
#     sol0 = x0_bubble_pressure(model,T,x)
#     vl0 = exp10(sol0[1])
#     vv0 = exp10(sol0[2])

#     @show y0 = collect(FractionVector(sol0[3:end]))
#     @show vl0 = volume(model,P,T,x,phase=:l)
#     @show vv0 = volume(model,P,T,y0,phase=:v)
#     @show μ_l = vt_chemical_potential(model,vl0,T,x)
#     @show μ_v = vt_chemical_potential(model,vv0,T,y0)
#     y1 =  μ_l ./ μ_v .* x
#     @show y1 = y1 ./ sum(y1)
#     @show pl0 = pressure(model,vl0,T,x)
#     @show pv0 = pressure(model,vv0,T,y1)

#     @show P = (pl0 - pv0)/(log(pl0) - log(pv0))

#     @show vl = volume(model,P,T,x,phase=:l)
#     @show vv = volume(model,P,T,y1,phase=:v)
#     #=
#     μ_l = vt_chemical_potential(model,vl,T,x)
#     μ_v = vt_chemical_potential(model,vv,T,y)
#     K = log.(μ_v) ./ log.(μ_l)
#     y = K .* x
#     y = y./sum(y)
#     pl = pressure(model,vl,T,x)
#     pv = pressure(model,vv,T,y)
#     P = (pl+pv)/2
#     @show vl,vv,y
#     return y,P
#     =#
# end


text:
using Test

module X
    using MongocUtils, Mongoc

    module E
        using MongocUtils, Mongoc

        struct N
            x::Int64
        end

        function a()
            n = N(1)
            return Mongoc.BSON(n)
        end
    end

    function b()
        n = E.N(2)
        return Mongoc.BSON(n)
    end
end

n = X.E.a()
n_st = X.E.as_struct(X.E.N, n)
@test n_st.x == 1

n = X.b()
n_st = X.as_struct(X.E.N, n)
@test n_st.x == 2
##


text:
module BOMBS

# Load dependencies
using BayesianOptimization
using BlackBoxOptim
using CSV
using Calculus
using CmdStan
using DataFrames
using Dates
using DiffEqBase
using DifferentialEquations
using Distributed
using Distributions
using GaussianProcesses
using GaussianMixtures
using JLD
using LinearAlgebra
using ODEInterfaceDiffEq
using OrdinaryDiffEq
using Plots
using Random
using Roots
using ScikitLearn.GridSearch
using ScikitLearnBase
using SharedArrays
# using Statistics
using StatsBase
using StatsPlots
using Sundials

# laod logo
include("logoBOMBS.jl")
export printLogo
export versionBOMBS

# load info function
include("InfoStructs.jl")
export infoAll

# Load section to generate the model scripts
include("ModelGen.jl")
export defModStruct
export checkStruct
export GenerateModel

# Load section to simulate the model
include("ModelSim.jl")
export simulateODEs
export defSimulStruct
export checkStructSimul
export fileStructInfo
export defSimulStructFiles
export extractSimulCSV
export plotSimsODE

# Load the section to generate pseudo-data
include("ModelPDat.jl")
export GenPseudoDat
export plotPseudoDatODE
export defPseudoDatStruct
export checkStructPseudoDat
# export fileStructInfo
export defPseudoDatStructFiles
export extractPseudoDatCSV
export PDatCSVGen

include("ModelMLE.jl")
export defMLEStruct
export SimToMle
export checkStructMLE
export selectObsSim_te
export restructInputs_te
export UVloglike
export MVloglike
export plotMLEResults
export defCrossValMLEStruct
export checkStructCrossValMLE
export plotCrossValMLEResults
export CrossValMLE
export finishMLEres
export MLEtheta

include("ModelStanInfer.jl")
export defBayInfStruct
export defBayInfDataStruct
export defBayInfDataFromFilesStruct
export defBasicStanSettingsStruct
export convertBoundTo2
export convertBoundToReal
export fitPriorSamps
export fitPriorSampsMultiNorm
export checkStructBayInf
export checkStructBayInfData
export checkStructBayInfDataFiles
export checkStructBayInfStanSettings
export genStanInitDict
export reparamDictStan
export genStanModel
export restructureDataInference
export getStanInferenceElements
export saveStanResults
export runStanInference
export plotStanResults
export StanInfer

include("ModelEntropyTheta.jl")
export genSamplesPrior
export computeH
export computeHgain

include("ModelOEDSelection.jl")
export defODEModelSelectStruct
export checkStructOEDMS
export BhattacharyyaDist
export EuclideanDist
export genOptimMSFuncts
export plotOEDMSResults
export settingsBayesOpt
export mainOEDMS

include("ModelOEDCalibration.jl")
export defODEModelCalibrStruct
export checkStructOEDMC
export genOptimMCFuncts
export plotOEDMCResults
export settingsBayesOptMC
export mainOEDMC

end

max_stars_repo_name: DavidGomezC/BOMBSbeta

max_issues_repo_name: DavidGomezC/BOMBSbeta

max_forks_repo_name: DavidGomezC/BOMBSbeta


text:
module PackageStates

using Pkg
using Dates
using PrettyTables

include("utils.jl")
include("core.jl")
include("API.jl")

function __init__()
    push!(Base.package_callbacks, on_load_package)
end

end # module


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(AbstractPlotting.backend_display, (GLBackend, Scene))

    # These are awful and will go stale as gensyms change (check by putting `@assert` in front of each one).
    # It would be far better to fix the inference problems.
    isdefined(GLMakie, Symbol("#89#90")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#89#90")),Text{Tuple{String}}})   # time: 1.7054044
    isdefined(GLMakie, Symbol("#72#78")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#72#78")),SMatrix{4, 4, Float32, 16},SMatrix{4, 4, Float32, 16}})   # time: 0.18958463
    Base.precompile(Tuple{typeof(to_gl_text),String,Vector{Point{3, Float32}},Vector{Float32},Vector{FTFont},Vec{2, Float32},Vector{Quaternionf0},SMatrix{4, 4, Float32, 16},Float64,Float64})   # time: 0.11106048
    isdefined(GLMakie, Symbol("#44#46")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#44#46")),GLFW.Window})   # time: 0.098638594
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Symbol})   # time: 0.07212781
    Base.precompile(Tuple{typeof(draw_atomic),Screen,Scene,Union{Scatter{ArgType} where ArgType, MeshScatter{ArgType} where ArgType}})   # time: 0.06695828
    isdefined(GLMakie, Symbol("#89#90")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#89#90")),LineSegments{Tuple{Vector{Point{2, Float32}}}}})   # time: 0.050924074
    isdefined(GLMakie, Symbol("#104#109")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#104#109")),String,Vector{Point{3, Float32}},Vector{Float32},Vector{FTFont},Vec{2, Float32},Vector{Quaternionf0},SMatrix{4, 4, Float32, 16},Float64,Float64})   # time: 0.038516257
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Bool})   # time: 0.030008739
    isdefined(GLMakie, Symbol("#101#102")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#101#102")),Int64,Point{3, Float32},Float32,FTFont,Vec{2, Float32}})   # time: 0.029477166
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Symbol})   # time: 0.019217245
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Billboard})   # time: 0.016408404
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Type})   # time: 0.006553374
    Base.precompile(Tuple{typeof(renderloop),Screen})   # time: 0.005448615
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),RGBA{N0f8}})   # time: 0.004190753
    Base.precompile(Tuple{typeof(push!),Screen,Scene,RenderObject{GLMakie.GLAbstraction.StandardPrerender}})   # time: 0.002534885
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Any})   # time: 0.002144091
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Any})   # time: 0.001724354
    isdefined(GLMakie, Symbol("#12#20")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#12#20"))})   # time: 0.00144086
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Vector{RGBA{Float32}}})   # time: 0.001202468
    Base.precompile(Tuple{typeof(draw_atomic),Screen,Scene,Lines{ArgType} where ArgType})   # time: 0.001200726
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Tuple{Symbol, Float64}})   # time: 0.001081687
    isdefined(GLMakie, Symbol("#81#82")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#81#82")),Vector{Float32}})   # time: 0.001065034
    isdefined(GLMakie, Symbol("#89#90")) && Base.precompile(Tuple{getfield(GLMakie, Symbol("#89#90")),Annotations{Tuple{Vector{Tuple{String, Point{2, Float32}}}}}})   # time: 0.001031033
end


text:
const SimpleDiGraphEdge = SimpleEdge

"""
    SimpleDiGraph{T}

A type representing a directed graph.
"""
mutable struct SimpleDiGraph{T<:Integer} <: AbstractSimpleGraph
    ne::Int
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
    badjlist::Vector{Vector{T}} # [dst]: (src, src, src)
end

eltype(x::SimpleDiGraph{T}) where T = T

# DiGraph{UInt8}(6), DiGraph{Int16}(7), DiGraph{Int8}()
function (::Type{SimpleDiGraph{T}})(n::Integer = 0) where T<:Integer
    fadjlist = Vector{Vector{T}}()
    badjlist = Vector{Vector{T}}()
    for _ = one(T):n
        push!(badjlist, Vector{T}())
        push!(fadjlist, Vector{T}())
    end
    return SimpleDiGraph(0, fadjlist, badjlist)
end

# DiGraph()
SimpleDiGraph() = SimpleDiGraph{Int}()

# DiGraph(6), DiGraph(0x5)
SimpleDiGraph(n::T) where T<:Integer = SimpleDiGraph{T}(n)

# SimpleDiGraph(UInt8)
SimpleDiGraph(::Type{T}) where T<:Integer = SimpleDiGraph{T}(zero(T))

# sparse adjacency matrix constructor: DiGraph(adjmx)
function (::Type{SimpleDiGraph{T}})(adjmx::SparseMatrixCSC{U}) where T<:Integer where U<:Real
    dima, dimb = size(adjmx)
    isequal(dima,dimb) || error("Adjacency / distance matrices must be square")

    g = SimpleDiGraph(T(dima))
    maxc = length(adjmx.colptr)
    for c = 1:(maxc-1)
        for rind = adjmx.colptr[c]:adjmx.colptr[c+1]-1
            isnz = (adjmx.nzval[rind] != zero(U))
            if isnz
                r = adjmx.rowval[rind]
                add_edge!(g,r,c)
            end
        end
    end
    return g
end

# dense adjacency matrix constructor: DiGraph{UInt8}(adjmx)
function (::Type{SimpleDiGraph{T}})(adjmx::AbstractMatrix) where T<:Integer
    dima,dimb = size(adjmx)
    isequal(dima,dimb) || error("Adjacency / distance matrices must be square")

    g = SimpleDiGraph(T(dima))
    for i in find(adjmx)
        ind = ind2sub((dima,dimb),i)
        add_edge!(g,ind...)
    end
    return g
end

# DiGraph(adjmx)
SimpleDiGraph(adjmx::AbstractMatrix) = SimpleDiGraph{Int}(adjmx)

# converts DiGraph{Int} to DiGraph{Int32}
function (::Type{SimpleDiGraph{T}})(g::SimpleDiGraph) where T<:Integer
    h_fadj = [Vector{T}(x) for x in fadj(g)]
    h_badj = [Vector{T}(x) for x in badj(g)]
    return SimpleDiGraph(ne(g), h_fadj, h_badj)
end

# constructor from abstract graph: DiGraph(graph)
function SimpleDiGraph(g::AbstractSimpleGraph)
    h = SimpleDiGraph(nv(g))
    h.ne = ne(g) * 2 - num_self_loops(g)
    h.fadjlist = deepcopy(fadj(g))
    h.badjlist = deepcopy(badj(g))
    return h
end

edgetype(::SimpleDiGraph{T}) where T<: Integer = SimpleGraphEdge{T}

badj(g::SimpleDiGraph) = g.badjlist
badj(g::SimpleDiGraph, v::Integer) = badj(g)[v]

copy(g::SimpleDiGraph{T}) where T<:Integer =
SimpleDiGraph{T}(g.ne, deepcopy(g.fadjlist), deepcopy(g.badjlist))

==(g::SimpleDiGraph, h::SimpleDiGraph) =
vertices(g) == vertices(h) &&
ne(g) == ne(h) &&
fadj(g) == fadj(h) &&
badj(g) == badj(h)

is_directed(g::SimpleDiGraph) = true
is_directed(::Type{SimpleDiGraph}) = true
is_directed(::Type{SimpleDiGraph{T}}) where T = true

function add_edge!(g::SimpleDiGraph, e::SimpleDiGraphEdge)
    T = eltype(g)
    s, d = T.(Tuple(e))
    (s in vertices(g) && d in vertices(g)) || return false
    inserted = _insert_and_dedup!(g.fadjlist[s], d)
    if inserted
        g.ne += 1
    end
    return inserted && _insert_and_dedup!(g.badjlist[d], s)
end

function rem_edge!(g::SimpleDiGraph, e::SimpleDiGraphEdge)
    has_edge(g,e) || return false
    i = searchsorted(g.fadjlist[src(e)], dst(e))[1]
    deleteat!(g.fadjlist[src(e)], i)
    i = searchsorted(g.badjlist[dst(e)], src(e))[1]
    deleteat!(g.badjlist[dst(e)], i)
    g.ne -= 1
    return true
end

function add_vertex!(g::SimpleDiGraph)
    T = eltype(g)
    (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
    push!(g.badjlist, Vector{T}())
    push!(g.fadjlist, Vector{T}())

    return true
end

function has_edge(g::SimpleDiGraph, e::SimpleDiGraphEdge)
    u, v = Tuple(e)
    u > nv(g) || v > nv(g) && return false
    if degree(g,u) < degree(g,v)
        return length(searchsorted(fadj(g,u), v)) > 0
    else
        return length(searchsorted(badj(g,v), u)) > 0
    end
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# intrinsic functions

# For curmod_*
include("testenv.jl")

# bits types
@test isa((() -> Core.Intrinsics.bitcast(Ptr{Int8}, 0))(), Ptr{Int8})
@test isa(convert(Char, 65), Char)

# runtime intrinsics
@testset "runtime intrinsics" begin
    @test Core.Intrinsics.add_int(1, 1) == 2
    @test Core.Intrinsics.sub_int(1, 1) == 0
    @test_throws ErrorException("fpext: output bitsize must be >= input bitsize")    Core.Intrinsics.fpext(Int32, 0x0000_0000_0000_0000)
    @test_throws ErrorException("fptrunc: output bitsize must be < input bitsize")  Core.Intrinsics.fptrunc(Int32, 0x0000_0000)
    @test_throws ErrorException("fptrunc: output bitsize must be < input bitsize")  Core.Intrinsics.fptrunc(Int64, 0x0000_0000)
    @test_throws ErrorException("ZExt: output bitsize must be > input bitsize")     Core.Intrinsics.zext_int(Int8, 0x00)
    @test_throws ErrorException("SExt: output bitsize must be > input bitsize")     Core.Intrinsics.sext_int(Int8, 0x00)
    @test_throws ErrorException("ZExt: output bitsize must be > input bitsize")     Core.Intrinsics.zext_int(Int8, 0x0000)
    @test_throws ErrorException("SExt: output bitsize must be > input bitsize")     Core.Intrinsics.sext_int(Int8, 0x0000)
    @test_throws ErrorException("Trunc: output bitsize must be < input bitsize")    Core.Intrinsics.trunc_int(Int8, 0x00)
    @test_throws ErrorException("Trunc: output bitsize must be < input bitsize")    Core.Intrinsics.trunc_int(Int16, 0x00)
end

# issue #4581
primitive type Date4581{T} 64 end
let
    x = Core.Intrinsics.bitcast(Date4581{Int}, Int64(1234))
    xs = Date4581[x]
    ys = copy(xs)
    @test ys !== xs
    @test ys == xs
end

# issue #6591
function f6591(d)
    Core.Intrinsics.bitcast(Int64, d)
    return (f -> f(d))(identity)
end
let d = Core.Intrinsics.bitcast(Date4581{Int}, Int64(1))
    @test isa(f6591(d), Date4581)
end

# test functionality of non-power-of-2 primitive type constants
primitive type Int24 24 end
Int24(x::Int) = Core.Intrinsics.trunc_int(Int24, x)
Int(x::Int24) = Core.Intrinsics.zext_int(Int, x)
let x, y, f
    x = Int24(Int(0x12345678)) # create something (via truncation)
    @test Int(0x345678) === Int(x)
    f() = Int24(Int(0x02468ace))
    y = f() # invoke llvm constant folding
    @test Int(0x468ace) === Int(y)
    @test x !== y
    @test string(y) == "$(curmod_prefix)Int24(0x468ace)"
end

# test nonsensical valid conversions and errors

compiled_addi(x, y) = Core.Intrinsics.add_int(x, y)
@test compiled_addi(C_NULL, C_NULL) === C_NULL
@test_throws ErrorException compiled_addi(C_NULL, 1)
@test_throws ErrorException ((x)->compiled_addi(Float64(x), x))(1)
@test ((x)->compiled_addi(Float64(x), Float64(x)))(2) === -0.0
@test compiled_addi(0.5, 5.0e-323) === 0.5000000000000011
@test_throws ErrorException compiled_addi(Int8(1), UInt8(1))
@test compiled_addi(UInt8(1), UInt8(2)) === UInt8(3)
@test_throws ErrorException compiled_addi(UInt8(1), UInt16(2))
@test compiled_addi(Float32(.125), Float32(10)) === 2.1267648f38
@test compiled_addi(true, true) === false

compiled_addf(x, y) = Core.Intrinsics.add_float(x, y)
@test compiled_addf(C_NULL, C_NULL) === C_NULL
@test_throws ErrorException compiled_addf(C_NULL, 1)
@test compiled_addf(0.5, 5.0e-323) === 0.5
@test_throws ErrorException compiled_addf(im, im)
@test_throws ErrorException compiled_addf(true, true)

function compiled_conv(::Type{T}, x) where T
    t = Core.Intrinsics.trunc_int(T, x)
    z = Core.Intrinsics.zext_int(typeof(x), t)
    s = Core.Intrinsics.sext_int(typeof(x), t)
    fpt = Core.Intrinsics.fptrunc(T, x)
    fpe = Core.Intrinsics.fpext(typeof(x), fpt)
    return (t, z, s, fpt, fpe)
end
@test compiled_conv(UInt32, Int64(0x8000_0000)) ==
    (0x80000000, Int64(0x80000000), -Int64(0x80000000), 0x00000000, 0)
@test compiled_conv(UInt32, UInt64(0xC000_BA98_8765_4321)) ==
    (0x87654321, 0x0000000087654321, 0xffffffff87654321, 0xc005d4c4, 0xc000ba9880000000)
@test_throws ErrorException compiled_conv(Bool, im)

let f = Core.Intrinsics.ashr_int
    @test f(Int8(-17), 1) == -9
    @test f(Int32(-1), 33) == -1
    @test f(Int32(-1), -1) == -1
    @test f(Int32(-1), -10) == -1
    @test f(Int32(2), -1) == 0
end

# issue #29929
@test unsafe_store!(Ptr{Nothing}(C_NULL), nothing) === Ptr{Nothing}(0)
@test unsafe_load(Ptr{Nothing}(0)) === nothing
struct GhostStruct end
@test unsafe_load(Ptr{GhostStruct}(rand(Int))) === GhostStruct()

# macro to verify and compare the compiled output of an intrinsic with its runtime version
macro test_intrinsic(intr, args...)
    output = args[end]
    inputs = args[1:end-1]
    quote
        function f()
            $intr($(inputs...))
        end
        @test f() === Base.invokelatest($intr, $(inputs...))
        @test f() == $output
    end
end

@testset "Float16 intrinsics" begin
    # unary
    @test_intrinsic Core.Intrinsics.neg_float Float16(3.3) Float16(-3.3)
    @test_intrinsic Core.Intrinsics.fpext Float32 Float16(3.3) 3.3007812f0
    @test_intrinsic Core.Intrinsics.fpext Float64 Float16(3.3) 3.30078125
    @test_intrinsic Core.Intrinsics.fptrunc Float16 Float32(3.3) Float16(3.3)
    @test_intrinsic Core.Intrinsics.fptrunc Float16 Float64(3.3) Float16(3.3)

    # binary
    @test_intrinsic Core.Intrinsics.add_float Float16(3.3) Float16(2) Float16(5.3)
    @test_intrinsic Core.Intrinsics.sub_float Float16(3.3) Float16(2) Float16(1.301)
    @test_intrinsic Core.Intrinsics.mul_float Float16(3.3) Float16(2) Float16(6.6)
    @test_intrinsic Core.Intrinsics.div_float Float16(3.3) Float16(2) Float16(1.65)
    @test_intrinsic Core.Intrinsics.rem_float Float16(3.3) Float16(2) Float16(1.301)

    # ternary
    @test_intrinsic Core.Intrinsics.fma_float Float16(3.3) Float16(4.4) Float16(5.5) Float16(20.02)
    @test_intrinsic Core.Intrinsics.muladd_float Float16(3.3) Float16(4.4) Float16(5.5) Float16(20.02)

    # boolean
    @test_intrinsic Core.Intrinsics.eq_float Float16(3.3) Float16(3.3) true
    @test_intrinsic Core.Intrinsics.eq_float Float16(3.3) Float16(2) false
    @test_intrinsic Core.Intrinsics.ne_float Float16(3.3) Float16(3.3) false
    @test_intrinsic Core.Intrinsics.ne_float Float16(3.3) Float16(2) true
    @test_intrinsic Core.Intrinsics.le_float Float16(3.3) Float16(3.3) true
    @test_intrinsic Core.Intrinsics.le_float Float16(3.3) Float16(2) false

    # conversions
    @test_intrinsic Core.Intrinsics.sitofp Float16 3 Float16(3f0)
    @test_intrinsic Core.Intrinsics.uitofp Float16 UInt(3) Float16(3f0)
    @test_intrinsic Core.Intrinsics.fptosi Int Float16(3.3) 3
    @test_intrinsic Core.Intrinsics.fptoui UInt Float16(3.3) UInt(3)
end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

"""
Calculate field values to nodal points from Gauss points using least-squares fitting.
"""
function calc_nodal_values!(elements::Vector, field_name, field_dim, time;
                            F=nothing, nz=nothing, b=nothing, return_F_and_nz=false)

    if F == nothing
        A = SparseMatrixCOO()
        for element in elements
            gdofs = get_connectivity(element)
            for ip in get_integration_points(element)
                detJ = element(ip, time, Val{:detJ})
                w = ip.weight*detJ
                N = element(ip, time)
                add!(A, gdofs, gdofs, w*kron(N', N))
            end
        end
        A = sparse(A)
        nz = get_nonzero_rows(A)
        A = 1/2*(A + A')
        F = ldlt(A[nz,nz])
    end

    if b == nothing
        b = SparseMatrixCOO()
        for element in elements
            gdofs = get_connectivity(element)
            for ip in get_integration_points(element)
                if !haskey(ip, field_name)
                    @warn("integration point does not have field $field_name")
                    continue
                end
                detJ = element(ip, time, Val{:detJ})
                w = ip.weight*detJ
                f = ip(field_name, time)
                N = element(ip, time)
                for dim=1:field_dim
                    add!(b, gdofs, [dim], w*f[dim]*N')
                end
            end
        end
        b = sparse(b)
    end

    x = zeros(size(b)...)
    x[nz, :] = F \ b[nz, :]
    nodal_values = Dict()
    for i=1:size(x,1)
        nodal_values[i] = vec(x[i,:])
    end
    update!(elements, field_name, time => nodal_values)
    if return_F_and_nz
        return F, nz
    end
end

"""
Return node ids + vector of values
"""
function get_nodal_vector(elements::Vector, field_name::AbstractString, time::Float64)
    f = Dict()
    for element in elements
        for (c, v) in zip(get_connectivity(element), element(field_name, time))
            if haskey(f, c)
                @assert isapprox(f[c], v)
            end
            f[c] = v
        end
    end
    node_ids = sort(collect(keys(f)))
    field = [f[nid] for nid in node_ids]
    return node_ids, field
end

"""
    problem(field_name, X, time)

Interpolate field from a set of elements defined in problem. Here, `X` is the
location inside domain described by elements.

Internally, function loops through all the elements, finding the one containing
the point `X`. After that, using inverse isoparametric mapping, first find
dimensionless coordinates (ξ,η,ζ) of that element corresponding to the location
of point `X` and after that interpolate the values of field under investigation.
Algorithm can be expected to be somewhat slow for big models, but for tests
models the performance is good.

# Examples

Having a problem called `body`, one can query the field `displacement` at
position `X = (1.0, 2.0, 3.0)` and time `t = 1.0`, with the command
```julia
X = (1.0, 2.0, 3.0)
time = 1.0
u = body("displacement", X, time)
```
"""
function (problem::Problem)(field_name, X, time; fillna=NaN)
    for element in get_elements(problem)
        if inside(element, X, time)
            xi = get_local_coordinates(element, X, time)
            return element(field_name, xi, time)
        end
    end
    return fillna
end

function (problem::Problem)(field_name, X, time, ::Type{Val{:Grad}}; fillna=NaN)
    for element in get_elements(problem)
        if inside(element, X, time)
            xi = get_local_coordinates(element, X, time)
            return element(field_name, xi, time, Val{:Grad})
        end
    end
    return fillna
end

function (solver::Solver)(field_name::AbstractString, X::Vector, time::Float64; fillna=NaN)
    for problem in get_problems(solver)
        for element in get_elements(problem)
            if inside(element, X, time)
                xi = get_local_coordinates(element, X, time)
                return element(field_name, xi, time)
            end
        end
    end
    return fillna
end

""" Calculate area of cross-section. """
function calculate_area(problem::Problem, X=[0.0, 0.0], time=0.0)
    A = 0.0
    for element in get_elements(problem)
        elsize = size(element)
        elsize[1] == 2 || error("wrong dimension of problem for area calculation, element size = $elsize")
        for ip in get_integration_points(element)
            w = ip.weight*element(ip, time, Val{:detJ})
            A += w
        end
    end
    return A
end

""" Calculate center of mass of body with respect to X.
https://en.wikipedia.org/wiki/Center_of_mass
"""
function calculate_center_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
    M = 0.0
    Xc = zero(X)
    for element in get_elements(problem)
        for ip in get_integration_points(element)
            w = ip.weight*element(ip, time, Val{:detJ})
            M += w
            rho = haskey(element, "density") ? element("density", ip, time) : 1.0
            Xp = element("geometry", ip, time)
            Xc += w*rho*(Xp-X)
        end
    end
    return 1.0/M * Xc
end

""" Calculate second moment of mass with respect to X.
https://en.wikipedia.org/wiki/Second_moment_of_area
"""
function calculate_second_moment_of_mass(problem::Problem, X=[0.0, 0.0, 0.0], time=0.0)
    n = length(X)
    I = zeros(n, n)
    for element in get_elements(problem)
        for ip in get_integration_points(element)
            w = ip.weight*element(ip, time, Val{:detJ})
            rho = haskey(element, "density") ? element("density", ip, time) : 1.0
            Xp = element("geometry", ip, time) - X
            I += w*rho*Xp*Xp'
        end
    end
    return I
end


text:
#
# oidata.jl --
#
# Implement OI-FITS data types.
#
#------------------------------------------------------------------------------
#
# This file is part of OIFITS.jl which is licensed under the MIT "Expat"
# License:
#
# Copyright (C) 2015-2019: Éric Thiébaut.
#
#------------------------------------------------------------------------------

###############
# BASIC TYPES #
###############

# Conversion to real type (no-op. if already of the correct type).
to_real(x::Cdouble) = x
to_real(x::Array{Cdouble}) = x
to_real(x::Array{T}) where {T<:Real} = convert(Array{Cdouble}, x)
to_real(x::Real) = convert(Cdouble, x)

to_integer(x::Int) = x
to_integer(x::Array{Int}) = x
to_integer(x::Array{T}) where {T<:Integer} = convert(Array{Int}, x)
to_integer(x::Integer) = convert(Int, x)

to_complex(x::Complex{Cdouble}) = x
to_complex(x::Array{Complex{Cdouble}}) = x
to_complex(x::AbstractArray{T}) where {T<:Union{Real,Complex}} =
    convert(Array{Complex{Cdouble}}, x)
to_complex(x::Union{Real,Complex}) = convert(Complex{Cdouble}, x)

# OI-FITS files stores the following 4 different data types:
const _DTYPE_LOGICAL = 1 # for format letter 'L'
const _DTYPE_INTEGER = 2 # for format letters 'I' or 'J'
const _DTYPE_REAL    = 3 # for format letters 'D' or 'E'
const _DTYPE_COMPLEX = 4 # for format letter 'C'
const _DTYPE_STRING  = 5 # for format letter 'A'

# The following dictionary is used for quick conversion of FITS format
# letter to data type.
const _DATATYPES = Dict('l' =>  _DTYPE_LOGICAL,
                        'L' =>  _DTYPE_LOGICAL,
                        'i' =>  _DTYPE_INTEGER,
                        'I' =>  _DTYPE_INTEGER,
                        'j' =>  _DTYPE_INTEGER,
                        'J' =>  _DTYPE_INTEGER,
                        'e' =>  _DTYPE_REAL,
                        'E' =>  _DTYPE_REAL,
                        'd' =>  _DTYPE_REAL,
                        'D' =>  _DTYPE_REAL,
                        'c' =>  _DTYPE_COMPLEX,
                        'C' =>  _DTYPE_COMPLEX,
                        'a' =>  _DTYPE_STRING,
                        'A' =>  _DTYPE_STRING)

is_logical(::Any) = false
is_logical(::Bool) = true
is_logical(::Array{Bool}) = true

is_string(::Any) = false
is_string(::AbstractString) = true
is_string(::Array{T}) where {T<:AbstractString} = true

is_integer(::Any) = false
is_integer(::Integer) = true
is_integer(::Array{T}) where {T<:Integer} = true

is_real(::Any) = false
is_real(::Real) = true
is_real(::Array{T}) where {T<:Real} = true

is_complex(::Any) = false
is_complex(::Complex) = true
is_complex(::Array{T}) where {T<:Complex} = true

#################
# OI-FITS TYPES #
#################

const OIContents = Dict{Symbol,Any}

@compat abstract type OIDataBlock end

mutable struct OITarget <: OIDataBlock
    owner
    contents::OIContents
    OITarget(contents::OIContents) = new(nothing, contents)
end

mutable struct OIArray <: OIDataBlock
    owner
    contents::OIContents
    OIArray(contents::OIContents) = new(nothing, contents)
end

mutable struct OIWavelength <: OIDataBlock
    owner
    contents::OIContents
    OIWavelength(contents::OIContents) = new(nothing, contents)
end

mutable struct OICorrelation <: OIDataBlock
    owner
    contents::OIContents
    OICorrelation(contents::OIContents) = new(nothing, contents)
end

mutable struct OIPolarization <: OIDataBlock
    owner
    arr::Union{OIArray,Nothing}
    contents::OIContents
    OIPolarization(contents::OIContents) = new(nothing, nothing, contents)
end

mutable struct OIVis <: OIDataBlock
    owner
    arr::Union{OIArray,Nothing}
    ins::Union{OIWavelength,Nothing}
    corr::Union{OICorrelation,Nothing}
    contents::OIContents
    OIVis(contents::OIContents) = new(nothing, nothing, nothing, nothing, contents)
end

mutable struct OIVis2 <: OIDataBlock
    owner
    arr::Union{OIArray,Nothing}
    ins::Union{OIWavelength,Nothing}
    corr::Union{OICorrelation,Nothing}
    contents::OIContents
    OIVis2(contents::OIContents) = new(nothing, nothing, nothing, nothing, contents)
end

mutable struct OIT3 <: OIDataBlock
    owner
    arr::Union{OIArray,Nothing}
    ins::Union{OIWavelength,Nothing}
    corr::Union{OICorrelation,Nothing}
    contents::OIContents
    OIT3(contents::OIContents) = new(nothing, nothing, nothing, nothing, contents)
end

mutable struct OISpectrum <: OIDataBlock
    owner
    arr::Union{OIArray,Nothing}
    ins::Union{OIWavelength,Nothing}
    corr::Union{OICorrelation,Nothing}
    contents::OIContents
    OISpectrum(contents::OIContents) = new(nothing, nothing, nothing, nothing, contents)
end

# Correspondance between OI-FITS data-block names and Julia types.
const _DATABLOCKS = Dict("OI_TARGET"     => OITarget,
                         "OI_WAVELENGTH" => OIWavelength,
                         "OI_ARRAY"      => OIArray,
                         "OI_VIS"        => OIVis,
                         "OI_VIS2"       => OIVis2,
                         "OI_T3"         => OIT3,
                         "OI_SPECTRUM"   => OISpectrum,
                         "OI_CORR"       => OICorrelation,
                         "OI_INSPOL"     => OIPolarization)

const _EXTNAMES = Dict{valtype(_DATABLOCKS), keytype(_DATABLOCKS)}()

function __init__()
    # Dictionaries whose keys are types may not be precompiled (at least for
    # Julia versions anterior to 0.5).
    for (key, val) in _DATABLOCKS
        haskey(_EXTNAMES, val) && error("values must be unique")
        _EXTNAMES[val] = key
    end
end

get_dbname(db::OIDataBlock) = _EXTNAMES[typeof(db)]

# OIData is any OI-FITS data-block which contains interferometric data.
const OIData = Union{OIVis,OIVis2,OIT3}

# OIDataBlock can be indexed by the name (either as a string or as a
# symbol) of the field.
getindex(db::OIDataBlock, key::Symbol) = get(db.contents, key, nothing)
getindex(db::OIDataBlock, key::AbstractString) = getindex(db, Symbol(key))
haskey(db::OIDataBlock, key::Symbol) = haskey(db.contents, key)
haskey(db::OIDataBlock, key::AbstractString) = haskey(db.contents, Symbol(key))
keys(db::OIDataBlock) = keys(db.contents)

# OIMaster stores the contents of an OI-FITS file.  Data-blocks containing
# measurements (OI_VIS, OI_VIS2 and OI_T3) are stored into a vector and
# thus indexed by an integer.  Named data-blocks (OI_ARRAY and
# OI_WAVELENGTH) are indexed by their names (converted to upper case
# letters, with leading and trailing spaces stripped, multiple spaces
# replaced by a single ordinary space).
mutable struct OIMaster
    all::Vector{OIDataBlock}            # All data-blocks
    update_pending::Bool                # Update is needed?
    tgt::Union{OITarget,Nothing}
    arr::Dict{Name,OIArray}
    ins::Dict{Name,OIWavelength}
    corr::Dict{Name,OICorrelation}
    function OIMaster()
        new(Array{OIDataBlock}(undef, 0),
            false,
            nothing,
            Dict{Name,OIArray}(),
            Dict{Name,OIWavelength}(),
            Dict{Name,OICorrelation}())
    end
end

function show(io::IO, db::OITarget)
    print(io, "OI_TARGET: ")
    if haskey(db, :target)
        tgt = db[:target]
        n = length(tgt)
        if n > 0
            print(io, "target=")
            for i in 1:n
                print(io, (i == 1 ? "[\"" : ", \""),
                      tgt[i], (i == n ? "\"]" : "\""))
            end
        end
    else
        n = 0
    end
    n == 0 && print(io, "<empty>")
end

function show(io::IO, db::OIArray)
    print(io, "OI_ARRAY: ")
    ntels = (haskey(db, :sta_index) ? length(db[:sta_index]) : 0)
    if haskey(db, :arrname)
        print(io, "arrname=\"", db[:arrname], "\", ", ntels, " telescope(s)")
    else
        print(io, ntels, " telescope(s)")
    end
end

function show(io::IO, db::OIWavelength)
    print(io, "OI_WAVELENGTH: ")
    if ! haskey(db, :insname)
        print(io, "<empty>")
    else
        print(io, "insname=\"", db[:insname], "\"")
        if haskey(db, :eff_wave)
            wave = db[:eff_wave]
            n = length(db[:eff_wave])
            if n < 1
                print(io, " with 0 spectral channels")
            elseif n == 1
                print(io, " with 1 spectral channel at ",
                      round(wave[1]*1e6, digits=3), " µm")
            else
                print(io, " with ", n, " spectral channels from ",
                      round(minimum(wave)*1e6, digits=3), " µm to ",
                      round(maximum(wave)*1e6, digits=3), " µm")
            end
        end
    end
end

function show(io::IO, db::OISpectrum)
    print(io, "OI_SPECTRUM")
end

function show(io::IO, db::OIDataBlock,
              name::AbstractString, nwaves::Integer, ntimes::Integer)
    print(io, name, ": ")
    if nwaves > 0 && ntimes > 0
        print(io, nwaves*ntimes, " measurements in ", nwaves,
              " spectral channel(s) and ", ntimes, " exposure(s)")
    else
        print(io, name, "<empty>")
    end
end

function show(io::IO, db::OIVis)
    if haskey(db, :visamp)
        dims = size(db[:visamp])
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    show(io, db, "OI_VIS", nwaves, ntimes)
end

function show(io::IO, db::OIVis2)
    if haskey(db, :vis2data)
        dims = size(db[:vis2data])
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    show(io, db, "OI_VIS2", nwaves, ntimes)
end

function show(io::IO, db::OIT3)
    if haskey(db, :t3amp)
        dims = size(db[:t3amp])
        nwaves = dims[1]
        ntimes = dims[2]
    else
        nwaves = 0
        ntimes = 0
    end
    show(io, db, "OI_T3", nwaves, ntimes)
end

function show(io::IO, master::OIMaster)
    print(io, "OI_MASTER: (", length(master.all), " data-block(s))")
    for db in master.all
        print(io, "\n    ")
        show(io, db)
    end
end

####################################
# PARSING OF DATABLOCK DEFINITIONS #
####################################

# OIFieldDef is used to store the definition of a keyword/column field.
mutable struct OIFieldDef
    name::Name      # Keyword/column name as a string.
    symb::Symbol    # Keyword/column symbolic name.
    keyword::Bool   # Is keyword? (otherwise column)
    optional::Bool  # Optional field?
    multiplier::Int # Multiplier: 1 for keywords, number of cells for columns
                    # (a negative number -N means an array of N dimensions each
                    # equal to the number of spectral channels.
    dtype::Int      # Data type.
    units::Name     # Units.
    descr::Name     # Description.
end

# OIDataBlockDef is used to store the definition of data-block.
mutable struct OIDataBlockDef
    dbname::Name
    fields::Vector{Symbol}        # ordered field symbolic names
    spec::Dict{Symbol,OIFieldDef} # dictionary of field specifications
    function OIDataBlockDef(dbname::AbstractString, vect::Vector{OIFieldDef})
        spec = Dict{Symbol,OIFieldDef}()
        fields = Array{Symbol}(undef, length(vect))
        for j in 1:length(vect)
            entry = vect[j]
            fields[j] = entry.symb
            spec[entry.symb] = entry
        end
        new(dbname, fields, spec)
    end
end

# OIFormatDef is used to store all the data-block definitions for a given
# revision number.
const OIFormatDef = Dict{Name,OIDataBlockDef}

# _FORMATS table is indexed by the datablock name and then by the revision
# number of the corresponding OI-FITS table.
const _FORMATS = Dict{Name,Vector{Union{OIDataBlockDef,Nothing}}}()

# The default format version number.
default_revision() = 2

# _FIELDS is a dictionary indexed by the data-block name (e,g., "OI_VIS"), each
# entry stores a set of its fields.
const _FIELDS = Dict{Name,Set{Symbol}}()

"""
    get_def(db, revn = default_revision())

yields the `OIDataBlockDef` for datablock of type `db` (e.g., "OI_TARGET") in
revison `revn` of OI-FITS standard.

"""
function get_def(dbname::AbstractString, revn::Integer = default_revision())
    if ! haskey(_FORMATS, dbname)
        error("unknown data-block: $db")
    end
    v = _FORMATS[dbname]
    if revn < 1 || revn > length(v) || v[revn] === nothing
        error("unsupported revision number: $revn")
    end
    return v[revn]
end

function get_def(db::OIDataBlock)
    get_def(get_dbname(db), get_revn(db))
end

function get_descr(db::OIDataBlock)
    name = get_dbname(db)
    revn = get_revn(db)
    return (name, revn, get_def(name, revn))
end

function add_def(dbname::AbstractString,
                 revn::Integer,
                 tbl::Vector{S}) where {S<:AbstractString}
    if (! startswith(dbname, "OI_") || dbname != uppercase(dbname)
        || occursin(" ", dbname))
        error("invalid data-block name: \"$dbname\"")
    end
    if revn < 1
        error("invalid revision number: $revn")
    end
    if haskey(_FORMATS, dbname) && revn <= length(_FORMATS[dbname]) &&
        _FORMATS[dbname][revn] != nothing
        error("data-block \"$dbname\" version $revn already defined")
    end

    fields = get(_FIELDS, dbname, Set{Symbol}())
    def = Array{OIFieldDef}(undef, 0)
    keyword = true
    for j in 1:length(tbl)
        row = strip(tbl[j])
        m = match(r"^([^ ]+) +([^ ]+) +(.*)$", row)
        if m === nothing
            if match(r"^-+$", row) === nothing
                error("syntax error in OI_FITS definition: \"$row\"")
            end
            keyword = false
            continue
        end
        name = uppercase(m.captures[1])
        symb = symbolicname(name)
        format = m.captures[2]
        descr = m.captures[3]
        optional = (format[1] == '?')
        i = (optional ? 2 : 1)
        dtype = get(_DATATYPES, format[i], nothing)
        if dtype === nothing
            error("invalid format type in OI_FITS definition: \"$row\"")
        end
        if keyword
            if length(format) != i
                error("invalid keyword format in OI_FITS definition: \"$row\"")
            end
            multiplier = 1
        else
            # Very naive code to parse the dimension list of the column format.
            if length(format) < i + 3 || format[i+1] != '(' ||
                format[end] != ')'
                error("missing column dimensions in OI_FITS definition: \"$row\"")
            end
            format = uppercase(format[i+2:end-1])
            if format == "W"
                multiplier = -1
            elseif format == "W,W"
                multiplier = -2
            else
                multiplier = parse(Int, format)
                if multiplier <= 0
                    error("invalid multiplier in OI_FITS definition: \"$row\"")
                end
            end
        end
        mp = match(r"^(.*[^ ]) +\[([^\]])\]$", descr)
        if mp === nothing
            units = ""
        else
            descr = mp.captures[1]
            units = mp.captures[2]
        end
        push!(def, OIFieldDef(name, symb, keyword, optional, multiplier,
                              dtype, units, descr))
        push!(fields, symb)
    end

    # Insert the data-block definition in the global table.
    if haskey(_FORMATS, dbname)
        v = _FORMATS[dbname]
    else
        v = Array{Union{OIDataBlockDef,Nothing}}(undef, 0)
        _FORMATS[dbname] = v
    end
    while revn > length(v)
        push!(v, nothing)
    end
    v[revn] = OIDataBlockDef(dbname, def)
    _FIELDS[dbname] = fields
end

# Convert the name of an OI-FITS keyword/column into a valid symbol.
function symbolicname(name::AbstractString)
    key = lowercase(name)
    if key == "oi_revn"
        return :revn
    else
        return Symbol(replace(key, r"[^a-z0-9_]" => '_'))
    end
end

#########################
# BUILDING OF DATABLOCK #
#########################

# Build constructors for the various data-blocks types.
for (func, name) in ((:new_target,     "OI_TARGET"),
                     (:new_array,      "OI_ARRAY"),
                     (:new_wavelength, "OI_WAVELENGTH"),
                     (:new_vis,        "OI_VIS"),
                     (:new_vis2,       "OI_VIS2"),
                     (:new_t3,         "OI_T3"),
                     (:new_spectrum,   "OI_SPECTRUM"),
                     (:new_corr,       "OI_CORR"),
                     (:new_inspol,     "OI_INSPOL"))
    @eval begin
        function $func(master::Union{OIMaster, Nothing}=nothing;
                       revn::Integer = default_revision(), kwds...)
            db = build_datablock($name, revn, kwds)
            if master != nothing
                attach!(master, db)
            end
            return db
        end
    end
end

function build_datablock(dbname::AbstractString, revn::Integer, kwds)
    def = get_def(dbname, revn)
    contents = OIContents()
    contents[:revn] = to_integer(revn)
    rows = -1         # number of rows in the OI-FITS table
    channels = -1     # number of spectral channels
    for (field, value) in kwds
        # Check whether this field exists.
        spec = get(def.spec, field, nothing)
        if spec === nothing
            error("data-block $dbname has no field \"$field\"")
        end

        # Check value type.
        if spec.dtype == _DTYPE_LOGICAL
            if ! is_logical(value)
                error("expecting boolean value for field \"$field\" in $dbname")
            end
        elseif spec.dtype == _DTYPE_INTEGER
            if ! is_integer(value)
                error("expecting integer value for field \"$field\" in $dbname")
            end
            value = to_integer(value)
        elseif spec.dtype == _DTYPE_REAL
            if ! is_real(value)
                error("expecting real value for field \"$field\" in $dbname")
            end
            value = to_real(value)
        elseif spec.dtype == _DTYPE_COMPLEX
            if ! is_complex(value)
                error("expecting complex value for field \"$field\" in $dbname")
            end
            value = to_complex(value)
        elseif spec.dtype == _DTYPE_STRING
            if ! is_string(value)
                error("expecting string value for field \"$field\" in $dbname")
            end
        else
            error("*** CORRUPTED WORKSPACE ***")
        end

        # Check value dimensions.
        dims = (isa(value, Array) ? size(value) : ())
        rank = length(dims)
        if spec.keyword
            if rank != 0
                error("expecting a scalar value for field \"$field\" in $dbname")
            end
        else
            # Check array rank.
            mult = spec.multiplier
            maxrank = (spec.dtype == _DTYPE_STRING || mult == 1 ? 1 :
                       mult == -2 ? 3 :
                       2)
            if rank > maxrank
                error("bad number of dimensions for field \"$field\" in $dbname")
            end

            # Fields may have up to 3 dimensions.  For now, the last dimension
            # is the number of rows.  FIXME: The ordering of dimensions must
            # be changed: the number of rows should be the frist dimension.
            dim0 = (rank >= 1 ? dims[end] : 1)
            dim1 = (rank >= 2 ? dims[1] : 1)
            dim2 = (rank >= 3 ? dims[2] : 1)

            if rows == -1
                rows = dim0
            elseif rows != dim0
                error("incompatible number of rows for field \"$field\" in $dbname")
            end

            if mult < 0
                # Expecting an N-by-W or N-by-W-by-W array (N is the number of
                # rows and W the number of channels).
                if mult == -2 && dim1 != dim2
                    error("bad dimensions for field \"$field\" in $dbname")
                end
                if channels == -1
                    channels = dim1
                elseif channels != dim1
                    error("incompatible number of spectral channels for field \"$field\" in $dbname")
                end
            elseif spec.dtype != _DTYPE_STRING && dim1 != mult
                error("bad dimensions for field \"$field\" in $dbname")
            end
        end

        # Store value.
        contents[field] = value
    end

    # Check that all mandatory fields have been given.
    nerrs = 0
    for field in def.fields
        spec = def.spec[field]
        if ! haskey(contents, field) && ! spec.optional
            @warn("missing value for field \"$field\" in $dbname")
            nerrs += 1
        end
    end
    if nerrs > 0
        error("some fields are missing in $dbname")
    end
    return _DATABLOCKS[dbname](contents)
end

###################################################
# MANAGING THE DATABLOCK CONTAINER (THE "MASTER") #
###################################################

# Letter case and trailing spaces are insignificant according to FITS
# conventions.
fixname(name::AbstractString) = uppercase(rstrip(name))

# OIDataBlock and OIMaster can be used as iterators.
@static if isdefined(Base, :iterate)
    import Base: iterate

    iterate(iter::OIDataBlock) = iterate(iter.contents)
    iterate(iter::OIDataBlock, state) = iterate(iter.contents, state)

    iterate(iter::OIMaster) = iterate(iter.all)
    iterate(iter::OIMaster, state) = iterate(iter.all, state)
else
    import Base: start, done, next

    start(iter::OIDataBlock) = start(iter.contents)
    done(iter::OIDataBlock, state) = done(iter.contents, state)
    next(iter::OIDataBlock, state) = next(iter.contents, state)

    start(iter::OIMaster) = start(update!(iter).all)
    done(iter::OIMaster, state) = done(iter.all, state)
    next(iter::OIMaster, state) = next(iter.all, state)
end

function select(master::OIMaster, args::AbstractString...)
    datablocks = Array{OIDataBlock}(undef, 0)
    for db in master
        if get_dbname(db) ∈ args
            push!(datablocks, db)
        end
    end
    return datablocks
end

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_target(mst)
```
yields the `OITarget` data-block of `mst` if any, `nothing` otherwise.

Assuming `tgt` is an instance of `OITarget`, then:
```
    get_target(tgt)
```
yields the "TARGET" column of `tgt` which is an array of target names.
"""
get_target(master::OIMaster) = update!(master).tgt

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_array(mst, arrname)
```
yields the `OIArray` data-block of `mst` whose name is `arrname` if any,
`nothing` otherwise.

Assuming `db` is an instance of a sub-type of `OIDataBlock`, then:
```
    get_array(db)
```
yields corresponding instance of `OIArray` if any, `nothing` otherwise.
"""
function get_array(master::OIMaster, arrname::AbstractString)
    get(update!(master).arr, fixname(arrname), nothing)
end
get_array(db::OIDataBlock) = nothing
get_array(db::Union{OIVis,OIVis2,OIT3}) = db.arr
get_array(db::OIArray) = db

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_instrument(mst, insname)
```
yields the `OIWavelength` data-block of `mst` whose name is `insname` if any,
`nothing` otherwise.

Assuming `db` is an instance of a sub-type of `OIDataBlock`, then:
```
    get_instrument(db)
```
yields corresponding instance of `OIWavelength` if any, `nothing` otherwise.
"""
function get_instrument(master::OIMaster, insname::AbstractString)
    get(update!(master).ins, fixname(insname), nothing)
end
get_instrument(db::OIDataBlock) = nothing
get_instrument(db::Union{OIVis,OIVis2,OIT3}) = db.ins
get_instrument(db::OIWavelength) = db

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_targets(mst)
```
yields the names of the targets defined in `mst`.
"""
function get_targets(master::OIMaster)
    tgt = get_target(master)
    return tgt === nothing ? Array{String}(undef, 0) : get_target(tgt)
end

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_arrays(mst)
```
yields the names of the interferometric arrays defined in `mst`.
"""
get_arrays(master::OIMaster) = collect(keys(update!(master).arr))

"""
Assuming `mst` is an instance of `OIMaster`, then:
```
    get_instruments(mst)
```
yields the names of the instruments defined in `mst`.
"""
get_instruments(master::OIMaster) = collect(keys(update!(master).ins))

function new_master(datablocks::OIDataBlock...)
    master = OIMaster()
    for db in datablocks
        attach!(master, db)
    end
    return update!(master)
end

function new_master(datablocks::Array{OIDataBlock})
    master = OIMaster()
    for db in datablocks
        attach!(master, db)
    end
    return update!(master)
end

function attach!(master::OIMaster, db::OIDataBlock)
    db.owner === nothing || error("data-block already attached")
    if isa(db, OITarget)
        if master.tgt != nothing
            error("only one OI_TARGET data-block can be attached")
        end
        master.tgt = db
    elseif isa(db, OIWavelength)
        insname = fixname(db[:insname])
        if haskey(master.ins, insname)
            error("master already have an OI_WAVELENGTH data-block with INSNAME=\"$insname\"")
        end
        master.ins[insname] = db
    elseif isa(db, OIArray)
        arrname = fixname(db[:arrname])
        if haskey(master.arr, arrname)
            error("master already have an OI_ARRAY data-block with ARRNAME=\"$arrname\"")
        end
        master.arr[arrname] = db
    elseif isa(db, OICorrelation)
        corrname = fixname(db[:corrname])
        if haskey(master.corr, corrname)
            error("master already have an OI_CORR data-block with CORRNAME=\"$corrname\"")
        end
        master.corr[corrname] = db
    end
    push!(master.all, db)
    master.update_pending = true
    db.owner = master
    nothing
end

function update!(master::OIMaster)
    if master.update_pending
        if master.tgt === nothing
            error("missing mandatory OI_TARGET data-block")
        end
        for db in master.all
            if haskey(db, :insname) && ! isa(db, Union{OIWavelength,OIPolarization})
                insname = fixname(db[:insname])
                db.ins = get(master.ins, insname, nothing)
                if db.ins === nothing
                    error("OI_WAVELENGTH data-block with INSNAME=\"$insname\" not found in master")
                end
            end
            if haskey(db, :arrname) && ! isa(db, OIArray)
                arrname = fixname(db[:arrname])
                db.arr = get(master.arr, arrname, nothing)
                if db.arr === nothing
                    @warn("OI_ARRAY data-block with ARRNAME=\"$arrname\" not found in master")
                end
            end
            if haskey(db, :corrname) && ! isa(db, OICorrelation)
                corrname = fixname(db[:corrname])
                db.corr = get(master.corr, corrname, nothing)
                if db.corr === nothing
                    @warn("OI_CORR data-block with CORRNAME=\"$corrname\" not found in master")
                end
            end
        end
        master.update_pending = false
    end
    return master
end

"""
### Clone a data-block

This method clones an existing data-block.  The array data are shared between
the clones but not the links (owner, target, instrument, etc.).  Only the
fields defined by OI-FITS standard are cloned.  This method is intended when
a data-block from a given `OIMaster` instance is to be attached to another
`OIMaster` instance.
"""
function clone(db::OIDataBlock)
    (name, revn, defn) = get_descr(db)
    data = Dict{Symbol,Any}()
    for (key, val) in db.contents
        if haskey(defn.spec, key)
            data[key] = val
        end
    end
    return build_datablock(name, revn, data)
end

"""
### Select data for a given target

The method `select_target` selects a subset of data corresponding to a given
target.   The general syntax is:
```
    out = select_target(inp, tgt)
```
where `inp` is the input data (can be an instance of `OIMaster` or of any
`OIDataBlock` sub-types), `tgt` is the target number or name.

The result `out` may share part of its contents with the input data `inp`.

The result may be `nothing` if the input contains no data for the given target.
"""
select_target(db::OIDataBlock, target::Integer) = clone(db)

function select_target(db::OITarget, target::Integer)
    (name, revn, defn) = get_descr(db)
    k = findfirst(id -> id == target, get_target_id(db))
    k == 0 && return
    data = Dict{Symbol,Any}()
    for (key, val) in db.contents
        spec = get(defn.spec, key, nothing)
        spec != nothing || continue
        if spec.keyword
            data[key] = db[key]
        else
            data[key] = db[key][k:k]
        end
    end
    return build_datablock(name, revn, data)
end

function select_target(db::OIData, target::Integer)
    (name, revn, defn) = get_descr(db)
    target_id = get_target_id(db)
    sel = find(id -> id == target, target_id)
    length(sel) > 0 || return
    data = Dict{Symbol,Any}()
    cpy = (length(sel) == length(target_id)) # just copy?
    for (key, val) in db.contents
        spec = get(defn.spec, key, nothing)
        spec != nothing || continue
        if cpy || spec.keyword
            data[key] = db[key]
        else
            # Copy a sub-array corresponding to the selection.  (The target is
            # specified for each last index.)
            src = db[key]
            @assert(size(src)[end] == length(target_id))
            rank = ndims(src)
            dims = ntuple(i -> i == rank ? length(sel) : size(src, i), rank)
            dst = Array{eltype(src)}(undef, dims)
            data[key] = dst
            if rank == 1
                for j in 1:length(sel)
                    dst[j] = src[sel[j]]
                end
            elseif rank == 2
                for j in 1:length(sel)
                    dst[:,j] = src[:,sel[j]]
                end
            else
                error("unexpected rank $rank")
            end
        end
    end
    return build_datablock(name, revn, data)
end

function select_target(master::OIMaster, target::Integer)
    result = OIMaster()
    for db in master
        sel = select_target(db, target)
        sel != nothing && attach!(result, sel)
    end
    return update!(result)
end

function select_target(master::OIMaster, target::AbstractString)
    db = get_target(master)
    if db != nothing
        k = findfirst(name -> name == target, get_target(db))
        k != 0 && return select_target(master, get_target_id(db)[k])
    end
    error("target name not found")
end

"""
### Select data at given wavelengths

The method `select_wavelength` selects a subset of data on the basis of their
wavelength.   The general syntax is:
```
    out = select_wavelength(inp, sel)
```
where `inp` is the input data (can be an instance of `OIMaster` or of any
`OIDataBlock` sub-types), `sel` is a predicate function which takes a
wavelength (in meters) as argument and returns true if this wavelength is to be
selected and false otherwise.

An alternative is:
```
    out = select_wavelength(inp, wmin, wmax)
```
to select wavelengths `w` such that `wmin ≤ w ≤ wmax`.  The wavelength
bounds are in meters.

The result `out` may share part of its contents with the input data `inp`.

The result may be `nothing` if the input contains no data at selected
wavelengths.
"""
function select_wavelength(inp::Union{OIMaster,OIDataBlock},
                           wavemin::Real, wavemax::Real)
    wmin = convert(Cdouble, wavemin)
    wmax = convert(Cdouble, wavemax)
    select_wavelength(inp, w -> wmin ≤ w ≤ wmax)
end

select_wavelength(db::OIDataBlock, selector::Function) = clone(db)

function select_wavelength(db::Union{OIWavelength,OIVis,OIVis2,OIT3,OISpectrum},
                           selector::Function)
    (name, revn, defn) = get_descr(db)
    wave = get_eff_wave(db)
    sel = find(selector, wave)
    length(sel) > 0 || return
    data = Dict{Symbol,Any}()
    cpy = (length(sel) == length(wave)) # just copy?
    iswave = (typeof(db) == OIWavelength)
    for (key, val) in db.contents
        spec = get(defn.spec, key, nothing)
        spec != nothing || continue
        if cpy || spec.keyword || (spec.multiplier >= 0 && ! iswave)
            data[key] = db[key]
        else
            # Copy a sub-array corresponding to the selection.  (The wavelength
            # corresponds to the first index.)
            src = db[key]
            @assert(size(src, 1) == length(wave))
            rank = ndims(src)
            dims = ntuple(i -> i == 1 ? length(sel) : size(src, i), rank)
            dst = Array{eltype(src)}(undef, dims)
            data[key] = dst
            if rank == 1
                for j in 1:length(sel)
                    dst[j] = src[sel[j]]
                end
            elseif rank == 2
                for j in 1:length(sel)
                    dst[j,:] = src[sel[j],:]
                end
            else
                error("unexpected rank $rank")
            end
        end
    end
    return build_datablock(name, revn, data)
end

function select_wavelength(master::OIMaster, selector::Function)
    result = OIMaster()
    for db in master
        sel = select_wavelength(db, selector)
        sel != nothing && attach!(result, sel)
    end
    return update!(result)
end


text:
@testset "Part: 1" begin
    r = part(red, 1)
    g = part(green, 1)
    y = part(yellow, 1)
    b = part(blue, 1)

    @test color(r) == red
    @test color(g) == green
    @test color(y) == yellow
    @test color(b) == blue

    for p in [r, g, y, b]
        @test base(p) == empty_points
        @test Point(0, 0, 1) ∈ entrances(p)
        @test Point(0, 0, 0) ∈ exits(p)
        @test junctions(p) == empty_points
        @test shape(p) == 1
        @test !is_base(p)
        @test !is_junction(p)
        @test !is_leaf(p)
        @test !is_root(p)
    end
end

@testset "Part: 2" begin
    y = part(yellow, 2)
    b = part(blue, 2)

    @test color(y) == yellow
    @test color(b) == blue

    for p in [y, b]
        @test Point(0, 0, 0) ∈ base(p)
        @test entrances(p) == empty_points
        @test exits(p) == empty_points
        @test Point(0, 0, 1) ∈ junctions(p)
        @test shape(p) == 2
        @test is_base(p)
        @test is_junction(p)
        @test !is_leaf(p)
        @test !is_root(p)
    end
end