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

using Markdown
using InteractiveUtils

# ╔═╡ 6c8602e6-6996-11eb-1df7-b9a9535458ae
begin
	using StatsPlots
	using Random
	using TimeSeries
	using Dates
	using Statistics
	using DataFrames
end

# ╔═╡ 0edc349e-6996-11eb-17fe-75948b35db2b
md"""
# Introduction to Time Series II

> Author: Edwin Bedolla
>
> Date: Original, 6th April 2020. This version, 7th February 2021.

In this document, the main statistics such as the **mean function**, **autocovariance function**
and **autocorrelation function** will be described, along with some examples.

We will import all of the necessary modules first.

"""

# ╔═╡ 75a9e6e4-6996-11eb-2bb2-892253de4b62
gr();

# ╔═╡ 9a0f887c-6996-11eb-204e-1da9ecf118b7
# Ensure reproducibility of the results
rng = MersenneTwister(8092);

# ╔═╡ 9df97102-6996-11eb-23da-43ac18f5284a
md"""
## Descriptive statistics and measures

A full description of a given time series is always given by the **joint distribution function**
of the time series which is a multi-dimensional function that is very difficult to track for
most of the time series that are dealt with.

Instead, we usually work with what's known as the **marginal distribution function** defined as

```math
F_t(x) = P \{ x_t \leq x \}
```

where ``P \{ x_t \leq x \}`` is the probability that the *realization*
of the time series ``x_t``
at time ``t`` is less or equal that the value of ``x``.
Even more common is to use a related function known as the **marginal density function**

```math
f_t(x) = \frac{\partial F_t(x)}{\partial x}
```

and when both functions exist they can provide all the information needed to do meaningful
analysis of the time series.

### Mean function

With these functions we can now define one of the most important descriptive measures,
the **mean function** which is defined as

```math
\mu_{xt} = E(x_t) = \int_{-\infty}^{\infty} x f_t(x) dx
```

where ``E`` is the *expected value operator* found in classical statistics.

### Autocovariance and autocorrelation

We are also interested in analyzing the dependence or lack of between
realization values in different
time periods, i.e. ``x_t`` and ``x_s``; in that case we can use classical
statistics to define two very important and fundamental quantities.

The first one is known as the **autocovariance function** and it's defined as

```math
\gamma_{x} (s, t) = \text{cov}(x_s,x_t) = E[(x_s - \mu_s)(x_t - \mu_t)]
```

where ``\text{cov}`` is the [covariance](https://en.wikipedia.org/wiki/Covariance)
as defined in classical statistics. A simple way of defining the *autocovariance*
is the following

> The **autocovariance** tells us about the *linear* dependence between two points on the same
> time series observed at different times.

Normally, we know from classical statistics that if for a given time series
``x_t`` we should have
``\gamma_{x} (s, t) = 0`` then it means that there is no linear
dependence between ``x_t`` and ``x_s`` at time periods ``t`` and ``s``;
but this does not mean that there is **no** relation between them
at all. For that, we need another measure that we describe below.

We now introduce the **autocorrelation function** (ACF) and it's defined as

```math
\rho(s,t) = \frac{\gamma_{x} (s, t)}{\sqrt{\gamma_{x} (s, s) \gamma_{x} (t, t)}}
```

which is a measure of **predictability** and we can define it in words as follows

> The **autocorrelation** measures the *linear predictability* of a given time series ``x_t`` at
> time ``t``, using values from the same time series but at time ``s``.

This measure is very much related to[Pearson's correlation coefficient](https://en.wikipedia.org/wiki/Pearson_correlation_coefficient)
from classical statistics, which is a way to measure the relationship between values.

The range of values for the ACF is ``-1 \leq \rho(s,t) \leq 1``;
when ``\rho(s,t) = 1`` it means that a *linear model* can perfectly describe the
realization of the time series at time ``t``
provided with the realization at time ``s``, e.g. a trend goind upwards,
on the other hand, if ``\rho(s,t) = -1`` would mean that the realization
of the time series ``x_t`` decrease while the realization ``x_s`` is increasing.

"""

# ╔═╡ a25d2f06-6997-11eb-0c96-23650463738f
md"""
## Example

Let's look at an example for the particular case of the *moving average*. We will
be working out the analytic form of the *autocovariance function* and *ACF* for the moving average
while also providing the same results numerically using `Julia`.

Recall the 3-valued moving average to be defined as

```math
v_t = \frac{1}{3} \left( w_{t-1} + w_{t} + w_{t+1} \right)
```

Let's plot the moving average again. We will create a very big time series for the sake
of numerical approximation below.

First, we create the white noise time series.

"""

# ╔═╡ b9fcb5f0-6997-11eb-1027-dda87c151112
# Create a range of time for a year, spaced evenly every 1 minute
dates = DateTime(2018, 1, 1, 1):Dates.Minute(1):DateTime(2018, 12, 31, 24);

# ╔═╡ e5e97644-6997-11eb-0104-07a8556f3c96
# Build a TimeSeries object with the specified time range and white noise
ts = TimeArray(dates, randn(rng, length(dates)));

# ╔═╡ f32d34ee-6997-11eb-1629-a7f7bede2305
# Create a DataFrame of the TimeSeries for easier handling
df_ts = DataFrame(ts);

# ╔═╡ d47822d4-6997-11eb-3d00-39a0fe53f1e2
md"""
Then, as before, we compute the 3-valued moving average.
"""

# ╔═╡ dc19ba7a-6997-11eb-1c18-2f592d5120ec
# Compute the 3-valued moving average
moving_average = moving(mean, ts, 3);

# ╔═╡ e1567d34-6997-11eb-180d-99d7c6d4ebd1
# Create a DataFrame of the TimeSeries for easier handling
df_average = DataFrame(moving_average);

# ╔═╡ fdb064f4-6997-11eb-18ff-ddfa35e5e89f
md"""
Recall what these look like in a plot. We just plot the first 100 elements in
the time series to avoid having a very cluttered plot.
"""

# ╔═╡ 1698af3a-6998-11eb-297b-8176c30e0da8
# Indices to plot
idxs = 1:100;

# ╔═╡ 03658032-6998-11eb-0225-bfbce55888f5
begin
	@df df_ts plot(:timestamp[idxs], :A[idxs], label = "White noise")
	@df df_average plot!(:timestamp[idxs], :A[idxs], label = "Moving average")
end

# ╔═╡ 0da7f61a-6998-11eb-1381-ff73d8bf8d9f
md"""
We are now ready to do some calculations. First, we invoke the definition of the
*autocovariance function* and apply it to the moving average

```math
\gamma_v(s,t)=\text{cov}(v_s,v_t)=
\text{cov}\{\frac{1}{3} \left( w_{t-1} + w_{t} + w_{t+1} \right),
\frac{1}{3} \left( w_{s-1} + w_{s} + w_{s+1} \right)\}
```

and now we need to look at some special cases.

- When ``s = t`` we now have the following

```math
\gamma_v(t,t)=\text{cov}(v_t,v_t)=
\text{cov}\{\frac{1}{3} \left( w_{t-1} + w_{t} + w_{t+1} \right),
\frac{1}{3} \left( w_{t-1} + w_{t} + w_{t+1} \right)\}
```

then, by the property of
[covariance of linear combinations](https://en.wikipedia.org/wiki/Covariance#Covariance_of_linear_combinations)
we have the following simplification

```math
\gamma_v(t,t)=\text{cov}(v_t,v_t)=
\frac{1}{9}\{\text{cov}(w_{t-1},w_{t-1}) + \text{cov}(w_{t},w_{t})
+ \text{cov}(w_{t+1},w_{t+1})\}
```

and because $\text{cov}(U,U) = \text{var}(U)$ for a random variable $U$, for a white
noise random variable we have $\text{var}(w_t)=\sigma^2_{wt}$, thus

```math
\gamma_v(t,t)=\text{cov}(v_t,v_t)= \frac{3}{9} \sigma^2_{wt}
```

In this case, recall that our white noise is normally distributed
``w_t \sim \mathcal{N}(0,\sigma^2_{wt})`` with ``\sigma^2_{wt}`` so the true expected
value is the following

"""

# ╔═╡ 67cf2bc2-6998-11eb-1c9f-f9b0bde636fe
true_γ = 3 / 9

# ╔═╡ 882be004-6998-11eb-058c-61714e97aa01
md"""
We will try to compute the *autocovariance function* using classical statistics
by means of the `cov` function in `Julia`. We need to pass it the time series like
so
"""

# ╔═╡ 91347c6a-6998-11eb-3de5-fff4ccaaca50
γ_jl = cov(df_average[:, :A], df_average[:, :A])

# ╔═╡ c78e9cb0-6997-11eb-23fa-9f1936729792
md"""
And we can see that the value is quite similar. The error must come from the fact
that we may need a bigger ensemble of values, but this should suffice.
"""

# ╔═╡ acafb36a-6998-11eb-3081-f9c9868f502d
md"""
- When ``s = t + 1`` we now have the following

```math
\gamma_v(t+1,t)=\text{cov}(v_{t+1},v_t)=
\text{cov}\{\frac{1}{3} \left( w_{t} + w_{t+1} + w_{t+2} \right),
\frac{1}{3} \left( w_{t-1} + w_{t} + w_{t+1} \right)\} \\
\gamma_v(t+1,t)=\frac{1}{9}\{\text{cov}(w_{t},w_{t}) + \text{cov}(w_{t+1},w_{t+1})\} \\
\gamma_v(t+1,t)=\frac{2}{9} \sigma^2_{wt}
```

So the true value is now
"""

# ╔═╡ bc0badd2-6998-11eb-3c7f-71441b60348b
true_γ1 = 2 / 9

# ╔═╡ c0d3f3a6-6998-11eb-1f6f-bf3d9f53c907
md"""
To check this, we perform the same operations as before, but this time, we need
to *move* the time series one time step with respect to itself.
"""

# ╔═╡ d077547e-6998-11eb-2d7f-5be66624cbb2
# Remove the last element from the first and start with the second element
γ_jl1 = cov(df_average[1:(end-1), :A], df_average[2:end, :A])

# ╔═╡ da11ef28-6998-11eb-38b9-c32091494df2
md"""
Great! Within a tolerance value, this is quite a nice estimate. It turns out that
for the cases ``s = t + h`` where ``h \geq 2``, the value for the *autocovariance*
is zero. We'll check it numerically here.
"""

# ╔═╡ e87fcae2-6998-11eb-33bd-19d5fc7c0ab3
# Remove the last element from the first and start with the second element
γ_jl_zero = cov(df_average[1:(end-3), :A], df_average[4:end, :A])

# ╔═╡ f1188e62-6998-11eb-0be9-d90fc94b3bf2
md"""
It's actually true, a value very close to zero but, _¿why?_ It's easy to see
if one applies the *autocovariance function* definition and checks the case
``s = t + 3``, and so on.

Let's now focus on the **ACF** for a 3-valued moving average.
We have several cases, like before.
"""

# ╔═╡ 05fc8600-6999-11eb-3edc-43b7a149d975
md"""
- When ``s = t`` we now have the following

```math
\rho_v(t,t)=\frac{\gamma_v(t,t)}{\sqrt{\gamma_v(t,t)\gamma_v(t,t)}}\\
\rho_v(t,t)=\frac{\gamma_v(t,t)}{\gamma_v(t,t)} = 1
```

so it turns out that the true value is ``\rho_v(t,t)=1``, and we can check this using
the `cor` function to compute the correlation coefficient in `Julia` as an estimate
for the *ACF*
"""

# ╔═╡ 17613e36-6999-11eb-0aa5-57b2d7d2ce0c
ρ_est = cor(df_average[:, :A], df_average[:, :A])

# ╔═╡ 222001f4-6999-11eb-33ca-83ec4484a65b
md"""
- When ``s = t + 1`` we now have the following

```math
\rho_v(t+1,t)=\frac{\gamma_v(t+1,t)}{\sqrt{\gamma_v(t+1,t+1)\gamma_v(t,t)}}\\
```

recall from before that ``\gamma(t,t)=3/9 \sigma_{vt}^2``
for a white noise time series, and we also have
``\gamma(t+1,t)=2/9 \sigma_{vt}^2``, so the *ACF* is now

```math
\rho_v(t+1,t)=\frac{2/9 \sigma_{vt}^2}{\sqrt{(3/9 \sigma_{vt}^2)(3/9 \sigma_{vt}^2)}}\\
\rho_v(t+1,t)=\frac{18 \sigma_{vt}^2}{27 \sigma_{vt}^2}\\
\rho_v(t+1,t)=\frac{2}{3}
```

which is the true value
"""

# ╔═╡ 45212f3e-6999-11eb-2b9c-dd5199e5f42a
true_ρ2 = 2 / 3

# ╔═╡ 4b74f596-6999-11eb-1bd2-1d6fce0f49b9
md"""
and again, we can check this value numerically
"""

# ╔═╡ 4fb34158-6999-11eb-214b-dd4ae90bb907
ρ_est2 = cor(df_average[1:(end-1), :A], df_average[2:end, :A])

# ╔═╡ 57325360-6999-11eb-37a0-fd1ba4b66d6f
md"""
Lastly, like with the *autocovariance*, the *ACF* for the cases
``s = t + h`` where ``h \geq 2``
is zero as seen below
"""

# ╔═╡ 6606538c-6999-11eb-1490-83a7727526d0
ρ_est_zero2 = cor(df_average[1:(end-3), :A], df_average[4:end, :A])

# ╔═╡ Cell order:
# ╟─0edc349e-6996-11eb-17fe-75948b35db2b
# ╠═6c8602e6-6996-11eb-1df7-b9a9535458ae
# ╠═75a9e6e4-6996-11eb-2bb2-892253de4b62
# ╠═9a0f887c-6996-11eb-204e-1da9ecf118b7
# ╟─9df97102-6996-11eb-23da-43ac18f5284a
# ╟─a25d2f06-6997-11eb-0c96-23650463738f
# ╠═b9fcb5f0-6997-11eb-1027-dda87c151112
# ╠═e5e97644-6997-11eb-0104-07a8556f3c96
# ╠═f32d34ee-6997-11eb-1629-a7f7bede2305
# ╟─d47822d4-6997-11eb-3d00-39a0fe53f1e2
# ╠═dc19ba7a-6997-11eb-1c18-2f592d5120ec
# ╠═e1567d34-6997-11eb-180d-99d7c6d4ebd1
# ╟─fdb064f4-6997-11eb-18ff-ddfa35e5e89f
# ╠═1698af3a-6998-11eb-297b-8176c30e0da8
# ╠═03658032-6998-11eb-0225-bfbce55888f5
# ╟─0da7f61a-6998-11eb-1381-ff73d8bf8d9f
# ╠═67cf2bc2-6998-11eb-1c9f-f9b0bde636fe
# ╟─882be004-6998-11eb-058c-61714e97aa01
# ╠═91347c6a-6998-11eb-3de5-fff4ccaaca50
# ╟─c78e9cb0-6997-11eb-23fa-9f1936729792
# ╟─acafb36a-6998-11eb-3081-f9c9868f502d
# ╠═bc0badd2-6998-11eb-3c7f-71441b60348b
# ╟─c0d3f3a6-6998-11eb-1f6f-bf3d9f53c907
# ╠═d077547e-6998-11eb-2d7f-5be66624cbb2
# ╟─da11ef28-6998-11eb-38b9-c32091494df2
# ╠═e87fcae2-6998-11eb-33bd-19d5fc7c0ab3
# ╟─f1188e62-6998-11eb-0be9-d90fc94b3bf2
# ╟─05fc8600-6999-11eb-3edc-43b7a149d975
# ╠═17613e36-6999-11eb-0aa5-57b2d7d2ce0c
# ╟─222001f4-6999-11eb-33ca-83ec4484a65b
# ╠═45212f3e-6999-11eb-2b9c-dd5199e5f42a
# ╟─4b74f596-6999-11eb-1bd2-1d6fce0f49b9
# ╠═4fb34158-6999-11eb-214b-dd4ae90bb907
# ╟─57325360-6999-11eb-37a0-fd1ba4b66d6f
# ╠═6606538c-6999-11eb-1490-83a7727526d0

max_issues_repo_name: BioLogs/TimeSeriesConcepts


text:
# Medidas de conjunto (ensemble) y combinación lineal

"""
    EnsembleInflationFunction <: InflationFunction <: InflationFunction

Tipo abstracto para conformar conjuntos de medidas de inflación
"""
abstract type EnsembleInflationFunction <: InflationFunction end

## EnsembleFunction

"""
    EnsembleFunction{N} <: EnsembleInflationFunction

    EnsembleFunction(inflfn1, inflfn2 [, ...])

Función de inflación para computar un conjunto de `N` de medidas de inflación simultáneamente utilizando las funciones `inflfn1, inflfn2, ...`.
"""
struct EnsembleFunction{N} <: EnsembleInflationFunction
    functions::NTuple{N, F where {F <: InflationFunction}}
end

"""
    InflationEnsemble <: EnsembleInflationFunction <: InflationFunction

Alias para [`EnsembleFunction`](@ref).
"""
const InflationEnsemble = EnsembleFunction

# Constructor a partir de funciones de inflación
EnsembleFunction(functions::Vararg{InflationFunction}) = EnsembleFunction(functions)

# Número de medidas del conjunto
num_measures(ensfn::EnsembleFunction) = sum(num_measures(inflfn) for inflfn in ensfn.functions)

# Nombres de las medidas del conjunto
function measure_name(ensfn::EnsembleFunction)
    reduce(vcat, [measure_name(inflfn) for inflfn in ensfn.functions])
end

# Método para obtener las trayectorias de inflación del conjunto
# Se computan para cada medida y se concatenan horizontalmente
function (ensfn::EnsembleFunction)(cs::CountryStructure, ::CPIVarInterm)
    mapreduce(inflfn -> inflfn(cs, CPIVarInterm()), hcat, ensfn.functions)::Matrix{eltype(cs)}
end

function (ensfn::EnsembleFunction)(cs::CountryStructure)
    mapreduce(inflfn -> inflfn(cs), hcat, ensfn.functions)::Matrix{eltype(cs)}
end

## Métodos de ayuda

"""
    components(inflfn::EnsembleFunction)

Devuelve un `DataFrame` con las componentes del `EnsembleFunction`.
"""
function components(inflfn::EnsembleFunction)
    DataFrame(measure = measure_name(inflfn))
end

ext: jl


text:
# name: titlepage

# class: middle, centre

# # Literate.jl
# ## Simple package for literate programming in Julia
# Fredrik Ekre, JuliaCon 2019
# .right[![](logo.png)]

# ---

# # What is literate programming?

# Definition from Wikipedia:

# > Literate programming is a programming paradigm [...] in which a program is given as an
# > explanation of the program logic in a natural language, interspersed
# > with snippets of macros and traditional source code, from which a compilable source
# > code can be generated.

# --

# Loosely translated by me:

# > A program with readable comments that explain the code.

# ---

# name: motivation-want

# # Motivation -- What do I want?
#
# - Documenting packages by examples, cf. "example notebooks"
#
#
# - Examples should be easily testable on CI (no outdated examples)
#
#
# - Multiple output formats:
#   - Markdown/HTML pages
#   - Jupyter notebooks
#   - Julia script files
#
#
# - Minimal maintenance -- one source, one place to make changes
#

# ---

# name: motivation-problem

# # Motivation -- What's the problem?

# - Maintaining multiple outputs is difficult and they will diverge
# - Notebooks are *really* bad in version control like `git`
#   - small changes result in unreadable and unreviewable diffs
#     ![](nb-diff.png)
#   - too "rich" format, includes large objects like images
#     (often in multiple formats)

# ---

# name: goals

# # Design goals

# A `Literate.jl` source file should
#
# - have simple syntax that Julia users are familiar with
#
#
# - be valid Julia syntax, e.g. a `.jl` file that can be `include`d
#
#
# - be able to produce multiple outputs -- markdown files, notebooks and regular scripts

# ---

# name: goals

# # Design goals

# A `Literate.jl` source file should
#
# - have simple syntax that Julia users are familiar with: Julia syntax!
#
#
# - be valid Julia syntax, e.g. a `.jl` file that can be `include`d: Everything that is not code behind Julia comments (`# `)
#
#
# - be able to produce multiple outputs -- markdown files, notebooks and regular scripts: Multiple source transformations

# ---

# name: offer

# # What does `Literate.jl` offer?
#
# `Literate.jl` provides three functions:
#
# - `Literate.markdown`:
#   Transforms the source file to a markdown file, to be used with e.g. Documenter
#
#
# - `Literate.notebook`:
#   Transforms the source file to a Jupyter notebook
#
#
# - `Literate.script`:
#   Transforms the source file to a "pure" Julia script file
#

# ---

# # Example of a source file
# The following code block is an example source file, and also the source block for the upcoming slide
# ```julia
# # # Example of a `Literate.jl` source file
# #
# # This is a Julia comment, so this is treated as markdown input.
# # Here is some **bold text**, and some *text in italic*.
# # Let's create a function in a Julia code block
#
# f(x) = 2 * exp(x)
# g(x) = 2 * log(x)
#
# # If we plot `f` and `g` with e.g. Plots.jl, the output will be captured
#
# using Plots
# x = range(0.0, stop = 2.0, length = 100)
# plot(x, [f.(x) g.(x)]; labels = ["f(x) = 2exp(x)" "g(x) = 2log(x)"], legend = :topleft)
# ```
#
# The upcoming two slides are the output of the the above example.

# ---

# # Example of a `Literate.jl` source file
#
# This is a Julia comment, so this is treated as markdown input.
# Here is some **bold text**, and some *text in italic*.
# Let's create a function in a Julia code block

f(x) = 2 * exp(x)
g(x) = 2 * log(x)

# ---

# If we plot `f` and `g` with e.g. Plots.jl, the output will be captured

using Plots
x = range(0.0, stop = 2.0, length = 100)
plot(x, [f.(x) g.(x)]; labels = ["f(x) = 2exp(x)" "g(x) = 2log(x)"], legend = :topleft)

# ---

# # Usage in the wild -- `Optim.jl`
#
# `Optim.jl` (package for minimization) uses Literate for some examples in the documentation.
#
# - [Source file](https://github.com/JuliaNLSolvers/Optim.jl/blob/master/docs/src/examples/ipnewton_basics.jl)
#
#
# - [HTML output](https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/)
#

# ---

# # Usage in the wild -- `@cormullion`'s blog
#
# Generated with Literate.jl + Thibaut Lienart's (`@tlienart`) [JuDoc.jl](https://github.com/tlienart/JuDoc.jl) package
#
# - [Source file](https://github.com/cormullion/cormullion.github.io/blob/master/src/source/pentachoron.jl)
#
#
# - [HTML output](https://cormullion.github.io/pub/2019-04-03-pentachoron.html)
#

# ---

# # Usage in the wild -- This presentation
#
# This presentation is generated with `Literate.jl`, and turned into slides
# using Pietro Vertechi's (`@piever`) [Remark.jl](https://github.com/piever/Remark.jl) package.

# ---

# # Thanks for listening!

# ---


text:
# This file contains functions for the Profitability Efficiency DEA model
"""
    ProfitabilityDEAModel
An data structure representing a profitability DEA model.
"""
struct ProfitabilityDEAModel <: AbstractProfitabilityDEAModel
    n::Int64
    m::Int64
    s::Int64
    alpha::Float64
    eff::Vector
    lambda::SparseMatrixCSC{Float64, Int64}
    crseff::Vector
    vrseff::Vector
    scaleff::Vector
    alloceff::Vector
end

"""
    deaprofitability(X, Y, W, P)
Compute profitability efficiency using data envelopment analysis for
inputs `X`, outputs `Y`, price of inputs `W`, and price of outputs `P`.

# Optional Arguments
- `alpha=0.5`: alpha to use for the generalized distance function.

# Examples
```jldoctest
julia> X = [5 3; 2 4; 4 2; 4 8; 7 9.0];

julia> Y = [7 4; 10 8; 8 10; 5 4; 3 6.0];

julia> W = [2 1; 2 1; 2 1; 2 1; 2 1.0];

julia> P = [3 2; 3 2; 3 2; 3 2; 3 2.0];

julia> deaprofitability(X, Y, W, P)
Profitability DEA Model
DMUs = 5; Inputs = 2; Outputs = 2
alpha = 0.5; Returns to Scale = VRS
─────────────────────────────────────────────────────────
   Profitability       CRS      VRS     Scale  Allocative
─────────────────────────────────────────────────────────
1       0.38796   0.636364  0.68185  0.93329     0.609651
2       1.0       1.0       1.0      1.0         1.0
3       0.765217  1.0       1.0      1.0         0.765217
4       0.25      0.25      0.25     1.0         1.0
5       0.15879   0.26087   0.36     0.724638    0.608696
─────────────────────────────────────────────────────────
```
"""
function deaprofitability(X::Matrix, Y::Matrix, W::Matrix, P::Matrix; alpha::Float64 = 0.5)::ProfitabilityDEAModel
    # Check parameters
    nx, m = size(X)
    ny, s = size(Y)

    nw, mw = size(W)
    np, sp = size(P)

    if nx != ny
        error("number of observations is different in inputs and outputs")
    end
    if nw != nx
        error("number of observations is different in input prices and inputs")
    end
    if np != ny
        error("number of observations is different in output prices and outputs")
    end
    if mw != m
        error("number of input prices and intputs is different")
    end
    if sp != s
        error("number of output prices and otuputs is different")
    end

    # Compute efficiency for each DMU
    n = nx

    pefficiency = zeros(n)
    plambdaeff = spzeros(n, n)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]
        w0 = W[i,:]
        p0 = P[i,:]

        # Create the optimization model
        deamodel = Model(with_optimizer(Ipopt.Optimizer, print_level= 0 ))
        @variable(deamodel, eff, start = 1.0)
        @variable(deamodel, lambda[1:n] >= 0)

        @NLobjective(deamodel, Min, eff)

        @NLconstraint(deamodel, sum(sum(W[t,mi] * X[t,mi] for mi in 1:m) / sum(P[t,si] * Y[t,si] for si in 1:s) * lambda[t] for t in 1:n) == eff * sum(w0[j] * x0[j] for j in 1:m ) / sum(p0[j] * y0[j] for j in 1:s))

        @constraint(deamodel, sum(lambda) == 1)

        # Optimize and return results
        JuMP.optimize!(deamodel)

        pefficiency[i]  = JuMP.objective_value(deamodel)
        plambdaeff[i,:] = JuMP.value.(lambda)

    end

    # Technical, scale and allocative efficiency
    crsefficiency = efficiency(deagdf(X, Y, alpha, rts = :CRS, slack = false))
    vrsefficiency = efficiency(deagdf(X, Y, alpha, rts = :VRS, slack = false))
    scalefficiency = crsefficiency ./ vrsefficiency
    allocefficiency = pefficiency ./ crsefficiency

    return ProfitabilityDEAModel(n, m, s, alpha, pefficiency, plambdaeff, crsefficiency, vrsefficiency, scalefficiency, allocefficiency)

end

function deaprofitability(X::Vector, Y::Matrix, W::Vector, P::Matrix; alpha::Float64 = 0.5)::ProfitabilityDEAModel
    X = X[:,:]
    W = W[:,:]
    return deaprofitability(X, Y, W, P, alpha = alpha)
end

function deaprofitability(X::Matrix, Y::Vector, W::Matrix, P::Vector; alpha::Float64 = 0.5)::ProfitabilityDEAModel
    Y = Y[:,:]
    P = P[:,:]
    return deaprofitability(X, Y, W, P, alpha = alpha)
end

function deaprofitability(X::Vector, Y::Vector, W::Vector, P::Vector; alpha::Float64 = 0.5)::ProfitabilityDEAModel
    X = X[:,:]
    W = W[:,:]
    Y = Y[:,:]
    P = P[:,:]
    return deaprofitability(X, Y, W, P, alpha = alpha)
end

function Base.show(io::IO, x::ProfitabilityDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    eff = efficiency(x)
    crseff = efficiency(x, :CRS)
    vrseff = efficiency(x, :VRS)
    scaleeff = efficiency(x, :Scale)
    alloceff = efficiency(x, :Allocative)

    if !compact
        print(io, "Profitability DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "alpha = ", x.alpha)
        print(io, "; Returns to Scale = VRS")
        print(io, "\n")
        show(io, CoefTable(hcat(eff, crseff, vrseff, scaleeff, alloceff), ["Profitability", "CRS", "VRS", "Scale", "Allocative"], ["$i" for i in 1:n]))
    else

    end
end


text:
struct Model_Record #record struct to associate a log_Li with a saved, calculated model
    path::String
    log_Li::AbstractFloat
end

struct ICA_PWM_Model #Independent component analysis position weight matrix model
    name::String #designator for saving model to posterior
    origin::String #functions instantiating IPMs should give an informative desc of the means by which the model was generated
    sources::Vector{Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer}} #vector of PWM signal sources (LOG PROBABILITY!!!) tupled with an index denoting the position of the first PWM base on the prior matrix- allows us to permute length and redraw from the appropriate prior position
    source_length_limits::UnitRange{<:Integer} #min/max source lengths for init and permutation
    mix_matrix::BitMatrix # obs x sources bool matrix
    log_Li::AbstractFloat
    permute_blacklist::Vector{Function} #blacklist of functions that ought not be used to permute this model (eg. because to do so would not generate a different model for IPMs produced from fitting the mix matrix)
    function ICA_PWM_Model(name, origin, sources, source_length_limits, mix_matrix, log_Li, permute_blacklist=Vector{Function}())
        verify_srcs(name, origin, sources)
        new(name, origin, sources, source_length_limits, mix_matrix, log_Li, permute_blacklist)
    end
end

#ICA_PWM_Model FUNCTIONS
ICA_PWM_Model(name::String, source_priors::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},<:Bool}}, mix_prior::Tuple{BitMatrix,<:AbstractFloat}, bg_scores::AbstractArray{<:AbstractFloat}, observations::AbstractArray{<:Integer}, source_length_limits::UnitRange{<:Integer}) = init_IPM(name, source_priors,mix_prior,bg_scores,observations,source_length_limits)

function verify_srcs(name, ori, sources::Vector{<:Tuple{<:AbstractMatrix{<:AbstractFloat},<:Integer}})
    for (n,(pwm,pi)) in enumerate(sources)
        !all(isprobvec(exp.(pwm[pos,:])) for pos in 1:size(pwm,1)) && throw(DomainError("Bad probvec in model $name, source $n, origin $(ori):$(exp.(pwm))"))
    end
end

#MODEL INIT
function init_IPM(name::String, source_priors::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},<:Bool}}, mix_prior::Tuple{BitMatrix,<:AbstractFloat}, bg_scores::AbstractArray{<:AbstractFloat}, observations::AbstractArray{<:Integer}, source_length_limits::UnitRange{<:Integer})
    assert_obs_bg_compatibility(observations,bg_scores)
    T,O = size(observations)
    S=length(source_priors)
    obs_lengths=[findfirst(iszero,observations[:,o])-1 for o in 1:size(observations)[2]]
    sources=init_logPWM_sources(source_priors, source_length_limits)
    mix=init_mix_matrix(mix_prior,O,S)
    log_lh = IPM_likelihood(sources, observations, obs_lengths, bg_scores, mix)

   return ICA_PWM_Model(name, "init_IPM", sources, source_length_limits, mix, log_lh)
end
                #init_IPM SUBFUNCS
                function assert_obs_bg_compatibility(obs, bg_scores)
                    T,O=size(obs)
                    t,o=size(bg_scores)
                    O!=o && throw(DomainError("Background scores and observations must have same number of observation columns!"))
                    T!=t+1 && throw(DomainError("Background score array must have the same observation lengths as observations!"))
                end

                function init_logPWM_sources(prior_vector::AbstractVector{<:Union{<:AbstractVector{<:Dirichlet{<:AbstractFloat}},<:Bool}}, source_length_limits::UnitRange{<:Integer})
                    srcvec = Vector{Tuple{Matrix{Float64},Int64}}()
                        for prior in prior_vector
                            if typeof(prior)<:AbstractVector{<:Dirichlet{<:AbstractFloat}}
                                if length(prior) < source_length_limits[1]
                                    length_dist=DiscreteNonParametric(
                                        [source_length_limits...],
                                        [PRIOR_LENGTH_MASS,ones(length(source_length_limits[2:end]))/length(source_length_limits[2:end]).*(1-PRIOR_LENGTH_MASS)])
                                    PWM_length=rand(length_dist)
                                elseif source_length_limits[1] < length(prior) < source_length_limits[end]
                                    length_dist=DiscreteNonParametric(
                                        [source_length_limits...],
                                        [
                                            (ones(length(prior)-source_length_limits[1]+1)/(length(prior)-source_length_limits[1]+1).*PRIOR_LENGTH_MASS)...,
                                            (ones(source_length_limits[end]-length(prior))/(source_length_limits[end]-length(prior)).*(1-PRIOR_LENGTH_MASS))...
                                        ]
                                    )
                                    PWM_length=rand(length_dist)
                                else
                                    PWM_length=rand(source_length_limits)
                                end

                                if PWM_length>length(prior)
                                    prior_coord=rand(-(PWM_length-length(prior)):1)
                                else
                                    prior_coord=rand(1:length(prior)-PWM_length+1)
                                end

                                PWM = zeros(PWM_length,4)

                                curr_pos=prior_coord
                                for pos in 1:PWM_length
                                    curr_pos < 1 || curr_pos > length(prior) ? dirichlet=Dirichlet(ones(4)/4) :
                                        dirichlet=prior[curr_pos]
                                    PWM[pos, :] = rand(dirichlet)
                                end
                                push!(srcvec, (log.(PWM), prior_coord)) #push the source PWM to the source vector with the prior coord idx to allow drawing from the appropriate prior dirichlets on permuting source length
                            elseif prior==false
                                PWM_length=rand(source_length_limits)
                                PWM=zeros(PWM_length,4)
                                dirichlet=Dirichlet(ones(4)/4)
                                for pos in 1:PWM_length
                                    PWM[pos,:] = rand(dirichlet)
                                end
                                push!(srcvec, (log.(PWM), 1))
                            else
                                throw(ArgumentError("Bad prior supplied for ICA_PWM_Model!"))
                            end
                        end
                    return srcvec
                end

                function init_mix_matrix(mix_prior::Tuple{BitMatrix,<:AbstractFloat}, no_observations::Integer, no_sources::Integer)
                    inform,uninform=mix_prior
                    if size(inform,2) > 0
                        @assert size(inform,1)==no_observations && size(inform,2)<=no_sources "Bad informative mix prior dimensions!"
                    end
                    @assert 0.0 <= uninform <=1.0 "Uninformative mix prior not between 0.0 and 1.0!"
                    mix_matrix = falses(no_observations, no_sources)
                    if size(inform,2)>0
                        mix_matrix[:,1:size(inform,2)]=inform
                    end
                    for index in CartesianIndices((1:no_observations,size(inform,2)+1:no_sources))
                        rand() <= uninform && (mix_matrix[index] = true)
                    end
                    return mix_matrix
                end

function Base.show(io::IO, m::ICA_PWM_Model; nsrc::Integer=length(m.sources), progress=false)
    nsrc == 0 && (nsrc=length(m.sources))
    nsrc>length(m.sources) && (nsrc=length(m.sources))
    nsrc==length(m.sources) ? (srcstr="All") : (srcstr="Top $nsrc")

    printidxs,printsrcs,printfreqs=sort_sources(m,nsrc)

    printstyled(io, "ICA PWM Model $(m.name) w/ logLi $(m.log_Li)\n", bold=true)
    println(io, srcstr*" sources:")
    for src in 1:nsrc
        print(io, "S$(printidxs[src]), $(printfreqs[src]*100)%: ")
        pwmstr_to_io(io, printsrcs[src])
        println(io)
    end
    println(m.origin)

    progress && return(nsrc+4)
end

function sort_sources(m, nsrc)
    printidxs=Vector{Integer}()
    printsrcs=Vector{Matrix{Float64}}()
    printfreqs=Vector{Float64}()

    freqs=vec(sum(m.mix_matrix,dims=1)); total=size(m.mix_matrix,1)
    sortfreqs=sort(freqs,rev=true); sortidxs=sortperm(freqs,rev=true)
    for srcidx in 1:nsrc
        push!(printidxs, sortidxs[srcidx])
        push!(printsrcs, m.sources[sortidxs[srcidx]][1])
        push!(printfreqs, sortfreqs[srcidx]/total)
    end
    return printidxs,printsrcs,printfreqs
end


text:
using DataStructures
using DimensionalData

import Base: getindex, keys, haskey, convert, show, IO, MIME

"""
DimensionalData provides the equivalent of an xarray `DataArray`, this wrapper
adds the functionality of an xarray `DataSet` by containing multiple
`DimensionalArray`s.

Object should be created by calling the relevant
`convert($(FUNCTIONNAME), dataset::DataSet)` method on a [`DataSet`](@ref DataSet).

# See also

[`Backend`](@ref Backend), [`ArrayWrapper`](@ref ArrayWrapper)
"""
struct DimensionalArrayWrapper <: ArrayWrapper
    "Dictionary of `DimensionName::String => DimensionLength::Int`"
    dimensions::OrderedDict
    "Named tuple of `DatasetName::Symbol => Dataset::DimensionalArray`"
    datasets::T where T <: NamedTuple
    "Dictionary of `AttributeName::String => AttributeValue::Any`"
    attributes::OrderedDict
    "Dictionary containing encoding information (usually `source`, `filter_by_keys`, and `ecode_cf`)"
    encoding::Dict

    #  Manually define inner constructor here so that it does not appear twice
    #  in the docs
    DimensionalArrayWrapper(dimensions, datasets, attributes, encoding) = new(dimensions, datasets, attributes, encoding)
end

Base.getindex(obj::DimensionalArrayWrapper, key) = getfield(obj, :datasets)[key]
Base.keys(obj::DimensionalArrayWrapper) = keys(getfield(obj, :datasets))
Base.haskey(obj::DimensionalArrayWrapper, key) = key in keys(obj)

function Base.getproperty(obj::DimensionalArrayWrapper, key::Symbol)
    if key in keys(obj)
        return getindex(obj, key)
    else
        return getfield(obj, key)
    end
end

function Base.convert(
    ::Type{DimensionalArrayWrapper},
    dataset::CfGRIB.DataSet
)::DimensionalArrayWrapper
    dimensions = dataset.dimensions
    attributes = dataset.attributes
    encoding = dataset.encoding

    multidimensional_idx =  (
        [size(v.data) for v in values(dataset.variables)]
        .|> length
        .|> x -> x > 1)
    multidimensional_keys = collect(keys(dataset.variables))[multidimensional_idx]
    multidimensional_values = [dataset.variables[k] for k in multidimensional_keys]

    shared_dimensions = [dataset.variables[k] for k in keys(dataset.dimensions)]
    shared_axis = [Dim{Symbol(k)}(dataset.variables[k].data) for k in keys(dataset.dimensions)]

    datasets = NamedTuple{Tuple(Symbol.(multidimensional_keys))}((
        DimensionalArray(CfGRIB.convert(Array, dataset.variables[k].data), Tuple(shared_axis))
        for k in multidimensional_keys
    ))

    return DimensionalArrayWrapper(dimensions, datasets, attributes, encoding)
end

convert(::Type{DimensionalArray}, dataset::DataSet) = convert(DimensionalArrayWrapper, dataset)::DimensionalArrayWrapper


text:
module Bandana

include("paisley.jl")

using .Paisley

function entryway(args...)

  local epoch = strftime("h%s", time())

  local attune = "guitar"

  local layout = eadgbe

  function fabric(yarn::String="z0")
    try
      local seal = symbol(yarn)
      if sentinel(yarn) && haskey(codex, seal)
        local diadem = "$attune-$yarn-$epoch"
        local specie = layout(seal)
        println("\n\t", diadem)
        for course in specie
          println("\t", course)
        end
      else
        println("\n\t", yarn, '?')
      end
    catch anomaly
      println("\nCause $anomaly")
    end
  end

  function atrium(cargo::Array)
    for yarn in cargo
      fabric(yarn)
    end
  end

  function gamut()
    local art = collect(keys(codex))
    sort!(art)
    for its in art
      fabric(string(its))
    end
  end

  local argots = String[]

  for word in collect(args)
    push!(argots, lowercase(word))
  end

  local dexter = length(argots)

  if dexter > 0
    local orchid = argots[1]

    if dexter == 1
      fabric(orchid)

    elseif dexter > 1

      if orchid == "aug4th" ||
         orchid == "a4" ||
         orchid == "b5" ||
         orchid == "dim5th" ||
         ismatch(r"^(bf)+b?$", orchid) ||
         ismatch(r"^triton[ae]?[ls]?$", orchid)
        attune, layout = orchid, bfbfb
        shift!(argots)

      elseif orchid == "cello" ||
             orchid == "p5" ||
             ismatch(r"^.*gda.*$", orchid) ||
             ismatch(r"^mando.*", orchid) ||
             ismatch(r"^viol.*", orchid)
        attune, layout = orchid, cgdae
        shift!(argots)

      elseif orchid == "gtr" ||
             ismatch(r"^g.*it.*r$", orchid) ||
             ismatch(r"^.*dgbe.*$", orchid) ||
             ismatch(r"^.*gcea.*$", orchid) ||
             ismatch(r"^uk[ue](le)*$", orchid)
        attune, layout = orchid, eadgbe
        shift!(argots)

      elseif orchid == "fkbjdn" ||
             orchid == "maj3rd" ||
             orchid == "m3"
        attune, layout = orchid, fkbjdn
        shift!(argots)

      elseif orchid == "bass" ||
             orchid == "p4" ||
             ismatch(r"^dra[cg]o.*$", orchid) ||
             ismatch(r"^.*eadg.*$", orchid)
        attune, layout = orchid, dragon
        shift!(argots)

      else
        nothing
      end

      local cargo = copy(argots)

      local head = first(cargo)

      isequal(head, "every") ||
      isequal(head, "gamut") ? gamut() : atrium(cargo)
    end

  else
    catalog()
  end

end

if length(ARGS) > 0
  entryway(ARGS...)
else
  entryway()
end

println()

end # module

max_stars_repo_name: reidiiius/masquerade

max_issues_repo_name: reidiiius/masquerade

max_forks_repo_name: reidiiius/masquerade


text:
function get_current_state(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "getCurrentState", FieldODEStateAndDerivative, ())
end

function get_global_current_state(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "getGlobalCurrentState", FieldODEStateAndDerivative, ())
end

function get_global_previous_state(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "getGlobalPreviousState", FieldODEStateAndDerivative, ())
end

function get_interpolated_state(obj::AbstractFieldODEStateInterpolator, arg0::RealFieldElement)
    return jcall(obj, "getInterpolatedState", FieldODEStateAndDerivative, (RealFieldElement,), arg0)
end

function get_previous_state(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "getPreviousState", FieldODEStateAndDerivative, ())
end

function is_current_state_interpolated(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "isCurrentStateInterpolated", jboolean, ())
end

function is_forward(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "isForward", jboolean, ())
end

function is_previous_state_interpolated(obj::AbstractFieldODEStateInterpolator)
    return jcall(obj, "isPreviousStateInterpolated", jboolean, ())
end

function restrict_step(obj::AbstractFieldODEStateInterpolator, arg0::FieldODEStateAndDerivative, arg1::FieldODEStateAndDerivative)
    return jcall(obj, "restrictStep", AbstractFieldODEStateInterpolator, (FieldODEStateAndDerivative, FieldODEStateAndDerivative), arg0, arg1)
end

function taylor(::Type{AdamsFieldStateInterpolator}, arg0::FieldEquationsMapper, arg1::FieldODEStateAndDerivative, arg2::RealFieldElement, arg3::RealFieldElement, arg4::Vector{RealFieldElement}, arg5::Array2DRowFieldMatrix)
    return jcall(AdamsFieldStateInterpolator, "taylor", FieldODEStateAndDerivative, (FieldEquationsMapper, FieldODEStateAndDerivative, RealFieldElement, RealFieldElement, Vector{RealFieldElement}, Array2DRowFieldMatrix), arg0, arg1, arg2, arg3, arg4, arg5)
end


text:
@testset "Ensure that there are no non-stdlib dependencies" begin
    package_directories = String[
        dirname(dirname(pathof(DockerSandbox))),
        dirname(@__DIR__),
        dirname(dirname(@__FILE__)),
    ]
    for package_directory in package_directories
        package_project_filename = joinpath(package_directory, "Project.toml")
        predictmd_project = TOML.parsefile(package_project_filename)
        predictmd_direct_deps = predictmd_project["deps"]
        for (name, uuid) in pairs(predictmd_direct_deps)
            is_stdlib = Pkg.Types.is_stdlib(Base.UUID(uuid))
            if !is_stdlib
                @error("There is a non-stdlib dependency", name, uuid)
            end
            @test is_stdlib
        end
    end
end


text:
"""
Simple representation of the results of an optimization method.
"""
struct Response
    minimizer
    minimum::AbstractFloat
    iterations::Int
    converged::Bool
end

struct TreeNode
    value
    depth::Int

    function TreeNode(value, depth=0)
        new(value, depth)
    end
end

max_stars_repo_name: camirmas/DesignAutomation

max_issues_repo_name: camirmas/DesignAutomation

max_forks_repo_name: camirmas/DesignAutomation


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Ghostscript"
version = v"9.53.3"

# Collection of sources required to build
sources = [
    ArchiveSource(
        "https://github.com/ArtifexSoftware/ghostpdl-downloads/releases/download/gs$(version.major)$(version.minor)$(version.patch)/ghostscript-$(version).tar.gz", # URL
        "6eaf422f26a81854a230b80fd18aaef7e8d94d661485bd2e97e695b9dce7bf7f" # SHA256 hash
    ),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/ghostscript*

if [[ "${target}" == *-mingw* ]]; then
    # Patches adapted from
    # https://github.com/msys2/MINGW-packages/tree/d87e68ae356d773901e5c477854312e5de0548cf/mingw-w64-ghostscript
    # (we don't need to patch configure.ac)
    atomic_patch -p1 ../patches/001-mingw-build.patch
    atomic_patch -p1 ../patches/003-libspectre.patch
    atomic_patch -p1 ../patches/004-FT_CALLBACK_DEF-deprecated.patch
    # Patch adapted from
    # https://github.com/mxe/mxe/commit/3fee463ee24536144b053cc519cf963afb08dfb4
    atomic_patch -p1 ../patches/005-fix-conflicting-types-gp_local_arg_encoding_get_codepoint.patch
fi

# Specify the native compiler for the programs that need to be run on the host
export CCAUX=${CC_BUILD}

# configure the Makefiles
./configure --prefix=${prefix} \
    --build=${MACHTYPE} \
    --host=${target} \
    --without-x \
    --disable-contrib \
    --disable-cups

# create the binaries
make -j${nproc}

# install to prefixes
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

products = [
    ExecutableProduct("gs", :gs),
    # These are shell wrappers around gs, not binary executables
    FileProduct("bin/dvipdf", :dvipdf),
    FileProduct("bin/eps2eps", :eps2eps),
    FileProduct("bin/gsbj", :gsbj),
    FileProduct("bin/gsdj", :gsdj),
    FileProduct("bin/gsdj500", :gsdj500),
    FileProduct("bin/gslj", :gslj),
    FileProduct("bin/gslp", :gslp),
    FileProduct("bin/gsnd", :gsnd),
    FileProduct("bin/pdf2dsc", :pdf2dsc),
    FileProduct("bin/pdf2ps", :pdf2ps),
    FileProduct("bin/pf2afm", :pf2afm),
    FileProduct("bin/pfbtopfa", :pfbtopfa),
    FileProduct("bin/pphs", :pphs),
    FileProduct("bin/printafm", :printafm),
    FileProduct("bin/ps2ascii", :ps2ascii),
    FileProduct("bin/ps2epsi", :ps2epsi),
    FileProduct("bin/ps2pdf", :ps2pdf),
    FileProduct("bin/ps2pdf12", :ps2pdf12),
    FileProduct("bin/ps2pdf13", :ps2pdf13),
    FileProduct("bin/ps2pdf14", :ps2pdf14),
    FileProduct("bin/ps2pdfwr", :ps2pdfwr),
    FileProduct("bin/ps2ps", :ps2ps),
    FileProduct("bin/ps2ps2", :ps2ps2),
]

dependencies = Dependency[
]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Distributed
# addprocs(2)
# nprocs()
@everywhere using FRGRealTime
@everywhere import FRGRealTime.ThresholdFunctions as TF
@everywhere using DifferentialEquations
@everywhere using HCubature
@everywhere using Dierckx
using LaTeXStrings
using Plots
using BenchmarkTools
@everywhere path2="/home/tyy/Documents/FRGRealTime.jl/Applications/O(N) t-chanel Scheme2 Eu/"
@everywhere include(joinpath(path2,"struct.jl"))
@everywhere include(joinpath(path2,"flow.jl"))

include(joinpath(path2,"solver.jl"))

@everywhere config_spec.dtmax=50
@time fourpointdataeu=tchanelSolveFourPointParallel(
    0.1,
    150.0, lpakrang,
    x -> sol1(x)[1],
    x -> sol1(x)[2],
    RelambdaflowEutest,
    6400*2,
    u0 = [6*8.0],
    config=config_spec,
)

writedlm(joinpath(path2,"Data/k=100.dat"),hcat(fourpointdataeu...))

# yaxis = :log,
plot(fourpointdataeu[1],-fourpointdataeu[2],label =L"\mathrm{Im}\lambda",title="k=200MeV",xlabel=L"p_0")

max_stars_repo_path: Applications/O(N) t-chanel Scheme2 Eu/main.jl


text:
"""
Variables from a posterior distribution (typically generated by MCMC).

`keys` and `vars` are sorted by `keys`. Access outside this library
should always use keys.
"""
@auto_hash_equals immutable PosteriorDraws
    len::Int
    keys::Vector{Symbol}
    vars::Vector{Any}
    function PosteriorDraws(len::Int, keys::Vector{Symbol}, vars::Vector)
        length(keys) == length(vars) ||
            error(ArgumentError("Non-conforming keys and vars."))
        all(length(var) == len for var in vars) ||
            error("Non-conformable lengths")
        vars = Any[isa(v, Vector) ? v : convert(Vector, v) for v in vars]
        p = sortperm(keys)
        sorted_keys = keys[p]
        all(a < b for (a,b) in partition(sorted_keys, 2, 1)) ||
            error(ArgumentError("Duplicate or non-sorted keys."))
        new(len, sorted_keys, vars[p])
    end
end

# function PosteriorDraws(len::Int, keys, vars)
#     PosteriorDraws(len, convert(Vector{Symbol}, keys), convert(Vector{Any}, vars))
# end

PosteriorDraws(len::Int) = PosteriorDraws(len, Symbol[], Any[])

function PosteriorDraws(keys::Vector{Symbol}, vars::Vector)
    length(vars) ≥ 1 ||
        error(ArgumentError("Need at least one variable to determine length."))
    PosteriorDraws(length(vars[1]), keys, vars)
end

function PosteriorDraws{T}(key_var_pairs::Pair{Symbol, T}...)
    PosteriorDraws(first.([key_var_pairs...]), last.([key_var_pairs...]))
end

function show(io::IO, pd::PosteriorDraws)
    println(io, "PosteriorDraws with $(pd.len) observations")
    for (key,var) in zip(pd.keys, pd.vars)
        println(io, "    $(key) => $(eltype(var))")
    end
end

keys(pd::PosteriorDraws) = pd.keys

size(pd::PosteriorDraws) = (length(pd.keys), pd.len)
size(pd::PosteriorDraws, index) = size(pd)[index]

"Lookup `key` in `pd`, returning the corresponding indexes."
function _key2index(pd::PosteriorDraws, key::Symbol)
    index = searchsortedfirst(pd.keys, key)
    index > length(pd.keys) && error(ArgumentError("Variable $(key) not found"))
    index
end

function _key2index(pd::PosteriorDraws, keys::Vector{Symbol})
    [_key2index(pd, key) for key in keys]
end

_key2index(pd::PosteriorDraws, ::Colon) = indices(pd.keys,1)

function getindex(pd::PosteriorDraws, keys)
    keyinds = _key2index(pd, keys)
    if isa(keyinds, Int)
        pd.vars[keyinds]
    else
        PosteriorDraws(pd.len, pd.keys[keyinds], pd.vars[keyinds])
    end
end

function getindex(pd::PosteriorDraws, keys, drawinds)
    keyinds = _key2index(pd, keys)
    if isa(keyinds, Int)
        pd.vars[keyinds][drawinds]
    else
        PosteriorDraws(pd.keys[keyinds],
                       [var[drawinds] for var in pd.vars[keyinds]])
    end
end

ndims(pd::PosteriorDraws) = 2


text:
fast_scalar_indexing(::Type{<:CuArrays.CuArray}) = false
@inline allowed_getindex(x::CuArrays.CuArray,i...) = CuArrays.@allowscalar(x[i...])
@inline allowed_setindex!(x::CuArrays.CuArray,v,i...) = (CuArrays.@allowscalar(x[i...] = v))

function Base.setindex(x::CuArrays.CuArray,v,i::Int)
  _x = copy(x)
  allowed_setindex!(_x,v,i)
  _x
end

function restructure(x::CuArrays.CuArray,y)
  reshape(Adapt.adapt(parameterless_type(x),y),size(x)...)
end


text:
using H5PLEXOS

zipfiles = ["Model Base_8200 Solution.zip",
            "Model Base_8200 NoInterval Solution.zip",
            "Model DAY_AHEAD_NO_TX Solution.zip",
            "Model DAY_AHEAD_NO_TX Stochastic Solution.zip",
            "Model DAY_AHEAD_ALL_TX Solution.zip",
            "Model DA_h2hybrid_SCUC_select_lines_Test_1day Solution.zip",
            "Model DAY_AHEAD_PRAS Solution.zip"]

zipfiles = vcat(zipfiles, ["Model Test$i Solution.zip" for i in 1:9])

testfolder = dirname(@__FILE__) * "/"

for zipfile in zipfiles
    println(zipfile)
    zippath = testfolder * zipfile
    # TODO: Actually test things
    process(zippath, replace(zippath, ".zip" => ".h5"))
end


text:
import Base.*, Base./

"""
a closed-loop transfer function that relates an output `Y` and an input `U` in a feedback loop.

the resulting closed-loop transfer function is:

```
 Y      top
--- = --------
 U    1 + g_ol
```

# example

```julia
g_ol = 4 / (s + 1) * 2 / (s + 2)
top = 5 / (s + 4)
g = ClosedLoopTransferFunction(top, g_ol)
```

# Arguments
* `top::TransferFunction`: numerator
* `g_ol::TransferFunction`: open-loop transfer function
"""
struct ClosedLoopTransferFunction
	top::TransferFunction
	g_ol::TransferFunction
end

# algebra
*(tf_cl::ClosedLoopTransferFunction, tf::TransferFunction) = ClosedLoopTransferFunction(tf_cl.top * tf, tf_cl.g_ol)
*(tf::TransferFunction, tf_cl::ClosedLoopTransferFunction) = *(tf_cl::ClosedLoopTransferFunction, tf::TransferFunction)
*(tf_cl::ClosedLoopTransferFunction, x::Number) = ClosedLoopTransferFunction(tf_cl.top * x, tf_cl.g_ol)
*(x::Number, tf_cl::ClosedLoopTransferFunction) = *(tf_cl::ClosedLoopTransferFunction, x::Number)
/(tf_cl::ClosedLoopTransferFunction, tf::TransferFunction) = ClosedLoopTransferFunction(tf_cl.top / tf, tf_cl.g_ol)
/(tf_cl::ClosedLoopTransferFunction, x::Number) = ClosedLoopTransferFunction(tf_cl.top / x, tf_cl.g_ol)

# representation of a closed loop system in standard form.
# useful for converting into state space representation.
#        p_a(s) e^ {-θs}
#  ---------------------------
#   p_b(s) + p_c(s) e^ {-ϕs}
struct CLTFStandard
    # numerator
    p_a::Poly
    θ::Union{Float64, Int}
    # denominator
    p_b::Poly
    p_c::Poly
    ϕ::Union{Float64, Int}
end

function CLTFStandard(cl::ClosedLoopTransferFunction)
    # numerator
    p_a = cl.top.numerator * cl.g_ol.denominator
    θ = cl.top.time_delay
    # denominator
    p_b = cl.g_ol.denominator * cl.top.denominator
    p_c = cl.g_ol.numerator   * cl.top.denominator
    ϕ = cl.g_ol.time_delay
    return CLTFStandard(p_a, θ, p_b, p_c, ϕ)
end

order(cl::CLTFStandard) = degree(cl.p_b)
order(cl::ClosedLoopTransferFunction) = degree(CLTFStandard(cl).p_b)

function strictly_proper(cl::CLTFStandard)
    return degree(cl.p_b) > max(degree(cl.p_a), degree(cl.p_c))
end
strictly_proper(cl::ClosedLoopTransferFunction) = strictly_proper(CLTFStandard(cl))

"""
 dx
---- = A * x(t) + B * u(t) + C * x(t - ϕ)
 dt
y(t) = D * x(t)
we compute the matrices A, B, C, D here.
only works if strictly proper.
"""
function tf_to_ss(cl::CLTFStandard)
    @assert strictly_proper(cl) "closed-loop system not strictly proper!"

    a_i(i::Int) = cl.p_a[i]
    b_i(i::Int) = cl.p_b[i]
	c_i(i::Int) = cl.p_c[i]

    n = order(cl)

    A = zeros(n, n)
    for i = 1:n-1
        A[i, i + 1] = 1.0
    end
    for i = 1:n
        A[n, i] = - b_i(i-1) / b_i(n)
    end

    B = zeros(n, 1)
    B[n, 1] = 1 / b_i(n)

    C = zeros(n, n)
    for i = 1:n
        C[n, i] = c_i(i-1) / b_i(n)
    end

	D = zeros(1, n)
	for i = 1:n
		D[i] = a_i(i-1)
	end

    return A, B, C, D
end

# see sim.jl for doc string
function simulate(cl::ClosedLoopTransferFunction, final_time::Float64; nb_time_points::Int=100)
	cls = CLTFStandard(cl)

    if ! strictly_proper(cls)
        error("closed loop system is not strictly proper...")
    end

    # convert tf to state space form
    A, B, C, D = tf_to_ss(cls)

    x0 = deepcopy(B) # initial condition
	# h is history
	h(p, t) = zeros(order(cls))
    f(x, h, p, t) = A * x - C * h(p, t - cls.ϕ) # RHS of ODE (ignore p for params)
    prob = DDEProblem(f, x0, h, (0.0, final_time))
    sol = solve(prob, d_discontinuities=[0.0, cls.ϕ, cls.θ, cls.θ + cls.ϕ])

	t = vcat([-0.05 * final_time, -1e-5],
	range(1e-5, final_time, length=nb_time_points - 2))
	y = [NaN for i = 1:nb_time_points]
	for (i, t_i) in enumerate(t)
		if t_i < cls.θ
			y[i] = 0.0
		else
			y[i] = (D * sol(t_i - cls.θ))[1] # put in shifted time
		end
	end
    return DataFrame(t=t, output=y)
end


text:
len = length

using PyCall,DataFrames,RCall
unshift!(PyVector(pyimport("sys")["path"]), "")
@pyimport ncdata
filename = "dsmacc_out.nc"
filename = "methane.nc"
data = ncdata.get(filename)

specs = names!(DataFrame(data["spec"]),[Symbol(i)for i in data["sc"]])
rates = names!(DataFrame(data["rate"]),[Symbol(i)for i in data["rc"]])

rates = rates[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["rc"])]]
specs = specs[:,[Symbol(i) for i in filter(x->!ismatch(r"EMISS|DUMMY",x),data["sc"])]]

rates = rates[7:len(rates)]
rates = rates[:,[Symbol(i) for i in filter(x->sum(rates[x])>0,names(rates))]]

# reading in data and cleaning
joint = "\n"*(join(names(rates),"\n"))*"\n"

products= [split(i[4:len(i)],"+") for i in matchall(r"-->([A-z0-9+]*)",joint)]
reactants = [split(i[2:len(i)-1],"+") for i in matchall(r"\n([A-z0-9+]{1,60})([-->]{0,1})",joint)]
#=
products[789]
1-element Array{SubString{UTF8String},1}:
"DUMMY"

=#

node_set = Set([ i[2:len(i)-1] for i in matchall(r"[\n+>-]{1}[A-z0-9]+[\n+>-]{1}",joint) ])
nodes = [i for i in node_set]
dict = Dict( [[i,x] for (x,i) in enumerate(nodes)])

#df = DataFrame(source = Char[], target = Char[], weight = Float16[])

flux = []
edges = []
for i in 1:length(reactants)
  rcol=[]
  #println(i)
  for j in reactants[i]
    coeff = match( r"([\d\s\.]*)(\D[\d\D]*)", j)
    dummy = specs[Symbol(coeff[2])]
    try
      push!(rcol,parse(Float32,coeff[1]))
    catch
      push!(rcol,1.0f0)
    end

    prod = 1
    for k in rcol
      prod*=k
    end

    push!(flux, prod*rates[i])
    # make edges array
    for l in products[i]
        push!(edges,[match( r"([\d\s\.]*)(\D[\d\D]*)", j)[2],match( r"([\d\s\.]*)(\D[\d\D]*)", l)[2],i])
    end

  end
end

normalise(x) = x/norm(x)
t = 180

links = filter(i -> flux[i][t]>0 , 1:len(flux))
tflux = [log10(flux[i][t]) for i in links]
weight = 1+normalise(tflux)

newflux = [flux[i][t] for i in 1:len(flux)]
counter = 0
for i in links
  newflux[i] = weight[counter+=1]
end

edge = filter(i -> newflux[i[3]]>0 , edges)
source = [i[1] for i in edge]
target=[i[2] for i in edge]
weighted = [newflux[i[3]]+0.0001 for i in edge]

@rput source
@rput target
@rput weighted

R"library(igraph)"

R"el <- structure(list(V1 = source, V2 = target, weight = weighted), .Names = c('V1',
'V2', 'weight'), class = 'data.frame', row.names = c(1:$(length(edges))
))"

R"g <- graph.data.frame(el)"

R"net = g"

density = R"ecount(net)/(vcount(net)*(vcount(net)-1))"#edge_density(net, loops=F)"

ec = R"eigen_centrality(g, directed=TRUE, weights=E(g)$weight)$vector"

#liklyhood of returned pairs. reciprocity
reciprocity = R"reciprocity(net)"
# # Mutual, asymmetric, and null node pairs
#2*dyad_census(net)$mut/ecount(net) # Calculating reciprocity

transivity_global = R"transitivity(net, type='global')"# net is treated as an undirected networktransitivity(as.undirected(net, mode="collapse")) # same as above

trasivity_pernode = R"transitivity(net, type='local')"

triad_census =R"triad_census(net)"

diameter = R"diameter(net, directed=T)" #longest shortest path

diameter_nodes = R"get_diameter(net, directed=T)"

degree = R"degree(net, mode='all')"

deg_dist= R"degree_distribution(net, cumulative=T, mode='all')"#plot 0:max(degree) against 1-this

indegree = R"degree(net, mode='in')"
indegree_centralisation = R"centr_degree(net, mode='in', normalized=T)"

outdegree = R"degree(net, mode='out')"
outegree_centralisation = R"centr_degree(net, mode='out', normalized=T)"

closeness = R"closeness(net, mode='all')# ,weights=NA)"
closeness_centalisation =R"centr_clo(net, mode="all", normalized=T)"

ec = R"eigen_centrality(net, directed=TRUE, weights=E(net)$weight)$vector"
ec_centralisation = R"centr_eigen(net, directed=TRUE, normalized=T)"

betweenness = R"betweenness(net, directed=T)"
betweenness_edge = R"edge_betweenness(net, directed=T)"
betweenness_centralisation = R"centr_betweenness(net, directed=T)"

hubs = R"hub_score(net)$vector"
authorities = R"authority_score(net, weights=NA)$vector"

mean_distance = R"mean_distance(net, directed=T)"
all_shortest = R"distances(net)"

cocitation = R"cocitation(net)"

#######subgrahs and cliwues

R"net.sym <- as.undirected(net, mode= "collapse",edge.attr.comb=list(weight="sum", "ignore"))"#sum weight ignore all other attr
print("maybe dont sum edges but subtract?")

cliques = R"cliques(net.sym)"
cliques_sizes= R"sapply(cliques(net.sym), length)"
cliques_largest=R"largest_cliques(net.sym)"

###commmunity detection
#make membership a function

edge_betweenness =R"cluster_edge_betweenness(net)"#dendPlot(ceb, mode="hclust")
#edge_betweenness_classification = R"membership($(edge_betweenness))"
#modularity(ceb)
#crossing(ceb, net) # boolean vector: TRUE for edges across communities

kcore = R"coreness(net, mode="all")"

assortativity_degree = R"assortativity_degree(net, directed=T)"

pagerank= R"page_rank(g, algo = c('prpack', 'arpack', 'power'), vids = V(g),
     directed = TRUE, damping = 0.85, personalized = NULL, weights = NULL,
     options = NULL)$vector"

####clustering below:

label = R"a = cluster_label_prop(net)"
label = R"a = cluster_walktrap(g,step= 20000 )"
label = R"a = cluster_spinglass(g, spins=200)"
label = R"a = cluster_optimal(g)"

#label = R"a = cluster_louvain(as.undirected(net, mode= "collapse",edge.attr.comb=list(weight="mean", "ignore")))"
#label = R"a = cluster_leading_eigen(g, steps = -1, weights = NULL, start = NULL,  options = arpack_defaults, callback = NULL, extra = NULL,  env = parent.frame())"

label = R"a = cluster_edge_betweenness(g, weights = E(g)$weight, directed = TRUE,
  edge.betweenness = TRUE, merges = TRUE, bridges = TRUE,
  modularity = TRUE, membership = TRUE)"

f = open("julia.out","w")
for i in 1:length(label)
R"str = paste(unlist(a[$i]),collapse='\t')"
@rget str
if str != ""
write(f,str)
write(f,"\n")
end
end
close(f)

function cent(x)

  f = open("centrality.out","w")
  write(f,[string(i)*" " for i in names(x)])
  write(f,"\n")
  write(f,[string(x[i])*" " for i in 1:length(x)])
  close(f)

end

#=

We can extract the distances to a node or set of nodes we are interested in. Here we will get the
distance of every media from the New York Times.
dist.from.NYT <- distances(net, v=V(net)[media=="NY Times"], to=V(net), weights=NA)
 Set colors to plot the distances:
oranges <- colorRampPalette(c("dark red", "gold"))
col <- oranges(max(dist.from.NYT)+1)
col <- col[dist.from.NYT+1]
plot(net, vertex.color=col, vertex

We can also find the shortest path between specific nodes. Say here between MSNBC and the New
York Post:
news.path <- shortest_paths(net,
from = V(net)[media=="MSNBC"],
to = V(net)[media=="New York Post"],
output = "both") # both path nodes and edges
# Generate edge color variable to plot the path:
ecol <- rep("gray80", ecount(net))
ecol[unlist(news.path$epath)] <- "orange"
# Generate edge width variable to plot the path:
ew <- rep(2, ecount(net))
ew[unlist(news.path$epath)] <- 4
# Generate node color variable to plot the path:
vcol <- rep("gray40", vcount(net))
vcol[unlist(news.path$vpath)] <- "gold"
plot(net, vertex.color=vcol, edge.color=ecol,
edge.width=ew, edge.arrow.mode=0)

=#

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "CFITSIOBuilder"
version = v"3.45.0"

# Collection of sources required to build CFITSIO
sources = [
    "http://heasarc.gsfc.nasa.gov/FTP/software/fitsio/c/cfitsio3450.tar.gz" =>
    "bf6012dbe668ecb22c399c4b7b2814557ee282c74a7d5dc704eb17c30d9fb92e",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cfitsio
if [[ "${target}" == *freebsd* ]]; then
    ./configure --prefix=$prefix --host=$target --enable-reentrant CC=gcc
else
    ./configure --prefix=$prefix --host=$target --enable-reentrant
fi
make -j shared
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libcfitsio", :libcfitsio)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
#=
The code in here has been extracted from the JuliaParser.jl package
with license:

The JuliaParser.jl package is licensed under the MIT "Expat" License:

> Copyright (c) 2014: Jake Bolewski.
>
> Permission is hereby granted, free of charge, to any person obtaining
> a copy of this software and associated documentation files (the
> "Software"), to deal in the Software without restriction, including
> without limitation the rights to use, copy, modify, merge, publish,
> distribute, sublicense, and/or sell copies of the Software, and to
> permit persons to whom the Software is furnished to do so, subject to
> the following conditions:
>
> The above copyright notice and this permission notice shall be
> included in all copies or substantial portions of the Software.
>
> THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
> EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
> MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
> IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
> CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
> TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
> SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

const EOF_CHAR = typemax(Char)

function is_identifier_char(c::Char)
  c == EOF_CHAR && return false
    if (((c >= 'A') & (c <= 'Z')) | ((c >= 'a') & (c <= 'z')) | (c == '_') | ((c >= '0') & (c <= '9')) )
        return true
    end
    return false
end

function is_identifier_start_char(c::Char)
    c == EOF_CHAR && return false
    if (((c >= 'A') & (c <= 'Z')) | ((c >= 'a') & (c <= 'z')) | (c == '_'))
        return true
    end
    return false
end

eof(io::IO) = Base.eof(io)
eof(c::Char) = c === EOF_CHAR
readchar(io::IO) = eof(io) ? EOF_CHAR : read(io, Char)


text:
module FatoracaoLU

include("fat1.jl")

end # module


text:
"""
    Diskmargin

The notation follows "An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet

# Fields:
`α`: The disk margin
`ω0`: The worst-case frequency
`f0`: The destabilizing perturbation `f0` is a complex number with simultaneous gain and phase variation. This critical perturbation causes an instability with closed-loop pole on the imaginary axis at the critical frequency ω0
`δ0`: The uncertain element generating f0.
`γmin`: The lower real-axis intercept of the disk (analogous to classical lower gain margin).
`γmax`: The upper real-axis intercept of the disk (analogous to classical upper gain margin).
`ϕm`: is analogous to the classical phase margin.
`σ`: The skew parameter that was used to calculate the margin

Note, `γmax` and `ϕm` are in smaller than the classical gain and phase margins sicne the classical margins do not consider simultaneous perturbations in gain and phase.

The "disk" margin becomes a half plane for `α = 2` and an inverted circle for `α > 2`. In this case, the upper gain margin is infinite. See the paper for more details, in particular figure 6.
"""
struct Diskmargin
    α
    ω0
    f0
    δ0
    γmin
    γmax
    σ
    ϕm
    L
end

function Diskmargin(α, σ=0; ω0=mising, f0=missing, δ0=missing, L=nothing)
    d = Disk(; α, σ)
    γmin = d.γmin
    γmax = d.γmax
    ϕm = d.ϕm
    Diskmargin(α, ω0, f0, δ0, γmin, γmax, σ, ϕm, L)
end

function Base.getproperty(dm::Diskmargin, s::Symbol)
    s ∈ fieldnames(typeof(dm)) && return getfield(dm, s)
    if s ∈ (:margin, :alpha, :diskmargin)
        return dm.α
    elseif s ∈ (:gm, :gainmargin)
        return [dm.γmin, dm.γmax]
    elseif s === :phasemargin
        return dm.ϕm
    else
        throw(ArgumentError("$(typeof(dm)) has no property named $s"))
    end
end

Base.propertynames(dm::Diskmargin) = (fieldnames(typeof(dm))..., :margin, :gainmargin, :phasemargin)

function Base.show(io::IO, dm::Diskmargin)
    println(io, "Disk margin with:")
    println(io, "Margin: ", dm.α)
    println(io, "Frequency: ", dm.ω0)
    if dm.γmax < dm.γmin # In this case, we have an "inverted circle"
        println(io, "Gain margins: [$(dm.γmin), Inf]")
    else
        println(io, "Gain margins: [$(dm.γmin), $(dm.γmax)]")
    end
    println(io, "Phase margin: ", dm.ϕm)
    println(io, "Skew: ", dm.σ)
    println(io, "Worst-case perturbation: ", dm.f0)
end

"""
    Disk

Represents a perturbation disc in the complex plane. `Disk(0.5, 2)` represents all perturbations in the circle centered at 1.25 with radius 0.75, or in other words, a gain margin of 2 and a pahse margin of 36.9 degrees.

A disk can be converted to a Nyquist exclusion disk by `nyquist(disk)` and plotted using `plot(disk)`.

# Arguments:
- `γmin`: Lower intercept
- `γmax`: Upper intercept
- `c`: Center
- `r`: Radius
- `ϕm`: Angle of tangent line through origin.

If γmax < γmin the disk is inverted.
See [`diskmargin`](@ref) for disk margin computations.
"""
struct Disk{T}
    γmin::T
    γmax::T
    c::T
    r::T
    ϕm::T
end

Disk(args...) = Disk(promote(args...)...)

center_radius(γmin, γmax) = 1/2 * (γmax + γmin), 1/2 * (γmax - γmin)

function Disk(γmin, γmax)
    c, r = center_radius(γmin, γmax)
    Disk(γmin, γmax, c, r)
end

function Disk(γmin, γmax, c, r)
    if !isfinite(γmax)
        ϕm = 90.0
    else
        ϕm = (1 + γmin*γmax) / (γmin + γmax)
        ϕm = ϕm >= 1 ? 0 : ϕm <= -1 ? Inf : rad2deg(acos(ϕm))
    end
    Disk(γmin, γmax, c, r, ϕm)
end

function Disk(; α, σ)
    γmin = (2 - α*(1-σ)) / (2 + α*(1+σ))
    γmax = (2 + α*(1-σ)) / (2 - α*(1+σ))
    Disk(γmin, γmax)
end

Disk(dm::Diskmargin) = Disk(dm.γmin, dm.γmax)

"""
    nyquist(d::Disk)

Transform a `Disk` representing a diskmargin to a exclusion disk in the Nyquist plane. This can be useful for visualizing a diskmargin in the Nyquist plane.
"""
ControlSystems.nyquist(d::Disk) = Disk(-inv(d.γmin), -inv(d.γmax)) # translate the disk to a nyquist exclusion disk

"""
    diskmargin(L, σ = 0)
    diskmargin(L, σ::Real, ω)

Calculate the disk margin of LTI system `L`. `L` is supposed to be a loop-transfer function, i.e., it should be square. If `L = PC` the disk margin for output perturbations is computed, whereas if `L = CP`, input perturbations are considered. If the method `diskmargin(P, C, args...)` is used, both are computed.

The implementation and notation follows
"An Introduction to Disk Margins", Peter Seiler, Andrew Packard, and Pascal Gahinet
https://arxiv.org/abs/2003.04771

The margins are aviable as fields of the returned objects, see [`Diskmargin`](@ref).

# Arguments:
- `L`: A loop-transfer function.
- `σ`: If little is known about the distribution of gain variations then σ = 0 is a reasonable choice as it allows for a gain increase or decrease by the same relative amount. *The choice σ < 0* is justified if the gain can decrease by a larger factor than it can increase. Similarly, *the choice σ > 0* is justified when the gain can increase by a larger factor than it can decrease. *If σ = −1* then the disk margin condition is αmax = inv(MT). This margin is related to the robust stability condition for models with multiplicative uncertainty of the form P (1 + δ). If σ = +1 then the disk margin condition is αmax = inv(MS)
- `kwargs`: Are sent to the [`hinfnorm`](@ref) calculation
- `ω`: If a vector of frequencies is supplied, the frequency-dependent disk margin will be computed, see example below.

# Example:
```
L = tf(25, [1,10,10,10])
dm = diskmargin(L, 0)
plot(dm) # Plot the disk margin to illustrate maximum allowed simultaneous gain and phase variations.

nyquistplot(L)
plot!(dm, nyquist=true) # plot a nyquist exclusion disk. The Nyquist curve will be tangent to this disk at `dm.ω0`
nyquistplot!(dm.f0*L) # If we perturb the system with the worst-case perturbation `f0`, the curve will pass through the critical point -1.

## Frequency-dependent margin
w = exp10.(LinRange(-2, 2, 500))
dms = diskmargin(L, 0, w)
plot(w, dms)
```

See also [`ncfmargin`](@ref).
"""
function diskmargin(L::LTISystem, σ::Real=0; kwargs...)
    issiso(L) || return sim_diskmargin(L, σ)
    M = 1/(1 + L) + (σ-1)/2
    n,ω0 = hinfnorm(M; kwargs...)
    diskmargin(L, σ, ω0)
end

"""
    diskmargin(L::LTISystem, σ::Real, ω)

Calculate the diskmargin at a particular frequency or vector of frequencies. If `ω` is a vector, you get a frequency-dependent diskmargin plot if you plot the returned value.
See also [`ncfmargin`](@ref).
"""
diskmargin(L::LTISystem, σ::Real, ω::AbstractArray) = map(w->diskmargin(L, σ, w), ω)

function diskmargin(L::LTISystem, σ::Real, ω0::Real)
    issiso(L) || return sim_diskmargin(L, σ, [ω0])[]
    M = 1/(1 + L) + (σ-1)/2
    freq = isdiscrete(L) ? cis(ω0*L.Ts) : complex(0, ω0)
    Sω = evalfr(M, freq)[]
    αmax = 1/abs(Sω)
    δ0 = inv(Sω)
    dp = Disk(; α = αmax, σ)
    if δ0 == 2/(σ+1)  # S = 1, L = 0
        Diskmargin(αmax, ω0, Inf, δ0, 0, Inf, σ, dp.ϕm, L)
    else
        f0 = (2 + δ0*(1-σ)) / (2 - δ0*(1+σ))
        Diskmargin(αmax, ω0, f0, δ0, dp.γmin, dp.γmax, σ, dp.ϕm, L)
    end
end

## Plotting ====================================================================

γϕcurve(dm::Diskmargin; kwargs...) = γϕcurve(dm.α, dm.σ; kwargs...)

"""
    γϕcurve(α, σ; N = 200)

Internal function. Get the curve of extremal stable gain and phase perturbations. This function is called with a single diskmargin object is plotted.

- `N`: Number of points on the curve.
"""
function γϕcurve(α, σ; N = 200)
    θ = LinRange(0, π, N)
    f = @. (2 - α*cis(θ)*(1-σ)) / (2 + α*cis(θ)*(1+σ))
    @. (abs(f), abs(rad2deg(angle(f))))
end

@recipe function plot(d::Disk; nyquist = false)
    θ = LinRange(0, 2pi, 200)
    re, im = @. cos(θ), sin(θ)
    if nyquist
        d = ControlSystems.nyquist(d)
    end
    c,r = d.c, d.r
    @series begin
        fill --> true
        fillalpha --> 0.5
        @. r*re+c, r*im
    end
end

@recipe function plot(dm::Diskmargin; nyquist=false)
    if nyquist
        @series begin
            label --> "σ = $(dm.σ)"
            ControlSystems.nyquist(Disk(dm))
        end
    else
        γ, ϕ = γϕcurve(dm)
        @series begin
            title --> "Stable region for combined gain and phase variation"
            xguide --> "Gain variation"
            yguide --> "Phase variation"
            label --> "σ = $(dm.σ)"
            fill --> true
            fillalpha --> 0.5
            γ, ϕ
        end
    end
end

@recipe function plot(dm::AbstractVector{<:Diskmargin}; lower=true, phase=true)
    w = [dm.ω0 for dm in dm]
    layout --> (phase ? 2 : 1, 1)
    link --> :x
    @series begin
        subplot --> 1
        title --> "Gain margin"
        label --> (lower ? ["Lower" "Upper"] : "Upper")
        # xguide --> "Frequency"
        xscale --> :log10
        yscale --> :log10
        gma = getfield.(dm, :γmax)
        gmi = getfield.(dm, :γmin)
        # ylims --> (0, min(10, maximum(gma)))
        data = (lower ? [gmi gma] : gma)
        data = max.(0, data)
        replace!(data, 0 => -Inf)
        w, data
    end
    phase || return
    @series begin
        subplot --> 2
        title --> "Phase margin"
        xguide --> "Frequency"
        xscale --> :log10
        label --> ""
        ϕm = getfield.(dm, :ϕm)
        w, ϕm
    end
end

@recipe function plot(dm::AbstractVector{<:AbstractVector{Diskmargin}})
    reduce(hcat, dm)
end

@recipe function plot(dm::AbstractMatrix{Diskmargin}; lower=true, phase=true)
    w = getfield.(dm[:,1], :ω0)
    ny = size(dm, 2)
    length(w) == size(dm, 1) || throw(ArgumentError("Frequency vector and diskmargin vector must have the same lengths."))
    layout --> (phase ? 2 : 1, 1)
    link --> :x
    neg2inf(x) = x <= 0 ? Inf : x
    for i = 1:ny
        is = ny == 1 ? "" : "$i"
        @series begin
            subplot --> 1
            title --> "Gain margin"
            # label --> permutedims([["Lower $i" for i = 1:ny]; ["Upper $i" for i = 1:ny]])
            label --> (lower ? ["Upper"*is "Lower"*is] : "Upper"*is)
            # xguide --> "Frequency"
            xscale --> :log10
            yscale --> :log10
            gma = getfield.(dm[:, i], :γmax) .|> neg2inf
            gmi = getfield.(dm[:, i], :γmin) .|> neg2inf
            replace!(gmi, 0 => -Inf)
            # ylims --> (0, min(10, maximum(gma)))
            w, (lower ? [gmi gma] : gma)
        end
        phase || continue
        @series begin
            subplot --> 2
            title --> "Phase margin"
            xguide --> "Frequency"
            xscale --> :log10
            label --> is#permutedims(["$i" for i = 1:ny])
            ϕm = getfield.(dm[:, i], :ϕm)
            w, ϕm
        end
    end
end


text:
using SegyIO
using Test

include("test_read.jl")
include("test_SeisBlock.jl")
include("test_SeisCon.jl")
include("test_read_ensemble.jl")


text:
######################################################################
# Unit Tests for KirUtil's XCopy pattern
# -----
# Licensed under MIT License
module TestXCopy
using Test
using KirUtil

struct XCopyTest1
    counter::Int
    flag::Bool
end

struct XCopyTest2
    counter::Int
    flag::Bool
end

@xcopy XCopyTest1

@xcopy XCopyTest2
@xcopy_override XCopyTest2 :counter tpl.counter + 1
@xcopy_override XCopyTest2 :flag true

@testset "KirUtil XCopy" begin
    @testset "simple" begin
        let test = XCopyTest1(1, false)
            @test (test = xcopy(test)) == XCopyTest1(1, false)
            @test (test = xcopy(test)) == XCopyTest1(1, false)
        end
    end

    @testset "overridden" begin
        let test = XCopyTest2(1, false)
            @test (test = xcopy(test)) == XCopyTest2(2, true)
            @test (test = xcopy(test)) == XCopyTest2(3, true)
        end
    end
end
end # module TestXCopy