text:
# ==============================================================================
# Algorithme de Kung (extrait S_N d'un ensemble statique de points S de IR^2)

function getNonDominatedPoints(XFeas, YFeas)
    S = []
    for i=1:length(XFeas)
        push!(S, (XFeas[i] , YFeas[i]) )
    end
    sort!(S, by = x -> x[1])
    SN=[] ; push!(SN, S[1]) ; minYFeas = S[1][2]
    for i=2:length(XFeas)
        if S[i][2] < minYFeas
            push!(SN, S[i]) ; minYFeas = S[i][2]
        end
    end
    return SN
end

# ==============================================================================
# Extraction de l'ensemble bornant primal + sa frontiere de l'ensemble des points realisables

function ExtractEBP(XFeas, YFeas)
    X_EBP_frontiere = (Int64)[];   Y_EBP_frontiere = (Int64)[]
    X_EBP = (Int64)[] ;            Y_EBP = (Int64)[]

    if length(XFeas) > 0
        SN = getNonDominatedPoints(XFeas, YFeas)

        push!(X_EBP_frontiere, SN[1][1]) ;
        push!(Y_EBP_frontiere, ceil(Int64, 1.1 * maximum(YFeas)))

        for i in 1:length(SN)-1
            push!(X_EBP_frontiere, SN[i][1]);     push!(Y_EBP_frontiere, SN[i][2])
            push!(X_EBP_frontiere, SN[i+1][1]);   push!(Y_EBP_frontiere, SN[i][2])
            push!(X_EBP, SN[i][1]);               push!(Y_EBP, SN[i][2])
        end
        push!(X_EBP_frontiere, SN[end][1]);       push!(Y_EBP_frontiere, SN[end][2])

        push!(X_EBP_frontiere, ceil(Int64, 1.1 * maximum(XFeas)))
        push!(Y_EBP_frontiere, SN[end][2])

        push!(X_EBP, SN[end][1]);   push!(Y_EBP, SN[end][2])
    end
    return X_EBP_frontiere, Y_EBP_frontiere,   X_EBP, Y_EBP
end

max_stars_repo_name: ducanh-le/GravityMachine

max_issues_repo_name: ducanh-le/GravityMachine

max_forks_repo_name: ducanh-le/GravityMachine


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	function transmission(EP::Model, inputs::Dict, UCommit::Int, NetworkExpansion::Int, CapacityReserveMargin::Int)
This function establishes decisions, expressions, and constraints related to transmission power flows between model zones and associated transmission losses (if modeled).
The function adds transmission reinforcement or construction costs to the objective function. Transmission reinforcement costs are equal to the sum across all lines of the product between the transmission reinforcement/construction cost, $pi^{TCAP}_{l}$, times the additional transmission capacity variable, $\bigtriangleup\varphi^{max}_{l}$.
```math
\begin{aligned}
& \sum_{l \in \mathcal{L}}\left(\pi^{TCAP}_{l} \times \bigtriangleup\varphi^{max}_{l}\right)
\end{aligned}
```
Note that fixed O\&M and replacement capital costs (depreciation) for existing transmission capacity is treated as a sunk cost and not included explicitly in the GenX objective function.
Power flow and transmission loss terms are also added to the power balance constraint for each zone:
```math
\begin{aligned}
&	- \sum_{l\in \mathcal{L}}{(\varphi^{map}_{l,z} \times \Phi_{l,t})} - \frac{1}{2} \sum_{l\in \mathcal{L}}{(\varphi^{map}_{l,z} \times \beta_{l,t}(\cdot))}
\end{aligned}
```
Power flows, $\Phi_{l,t}$, on each line $l$ into or out of a zone (defined by the network map $\varphi^{map}_{l,z}$), are considered in the demand balance equation for each zone. By definition, power flows leaving their reference zone are positive, thus the minus sign is used for this term. Losses due to power flows increase demand, and one-half of losses across a line linking two zones are attributed to each connected zone. The losses function $\beta_{l,t}(\cdot)$ will depend on the configuration used to model losses (see below).
**Accounting for Transmission and Network Expansion Between Zones**
Transmission flow constraints are modeled using a 'transport method','' where power flow, $\Phi_{l,t}$, on each line (or more likely a `path' aggregating flows across multiple parallel lines) is constrained to be less than or equal to the line's maximum power transfer capacity, $\varphi^{max}_{l}$, plus any transmission capacity added on that line (for lines eligible for expansion in the set $\mathcal{E}$). The additional transmission capacity, $\bigtriangleup\varphi^{max}_{l} $, is constrained by a maximum allowed reinforcement, $\overline{\bigtriangleup\varphi^{max}_{l}}$, for each line $l \in \mathcal{E}$.
```math
\begin{aligned}
	% trasmission constraints
	&-\varphi^{max}_{l} \leq  \Phi_{l,t} \leq \varphi^{max}_{l} , &\quad \forall l \in (\mathcal{L} \setminus \mathcal{E} ),\forall t  \in \mathcal{T}\\
	% trasmission expansion
	&-(\varphi^{max}_{l} + \bigtriangleup\varphi^{max}_{l} ) \leq  \Phi_{l,t} \leq (\varphi^{max}_{l} + \bigtriangleup\varphi^{max}_{l} ) , &\quad \forall l \in \mathcal{E},\forall t  \in \mathcal{T}\\
	& \bigtriangleup\varphi^{max}_{l}  \leq \overline{\bigtriangleup\varphi^{max}_{l}}, &\quad \forall l \in \mathcal{E}
\end{aligned}
```
**Accounting for Transmission Losses**
Transmission losses due to power flows can be accounted for in three different ways. The first option is to neglect losses entirely, setting the value of the losses function to zero for all lines at all hours. The second option is to assume that losses are a fixed percentage, $\varphi^{loss}_{l}$, of the magnitude of power flow on each line, $\mid \Phi_{l,t} \mid$ (e.g., losses are a linear function of power flows). Finally, the third option is to calculate losses, $\ell_{l,t}$, by approximating a quadratic-loss function of power flow across the line using a piecewise-linear function with total number of segments equal to the size of the set $\mathcal{M}$.
```math
\begin{aligned}
%configurable losses formulation
	& \beta_{l,t}(\cdot) = \begin{cases} 0 & \text{if~} \text{losses.~0} \\ \\ \varphi^{loss}_{l}\times \mid \Phi_{l,t} \mid & \text{if~} \text{losses.~1} \\ \\ \ell_{l,t} &\text{if~} \text{losses.~2} \end{cases}, &\quad \forall l \in \mathcal{L},\forall t  \in \mathcal{T}
\end{aligned}
```
For the second option, an absolute value approximation is utilized to calculate the magnitude of the power flow on each line (reflecting the fact that negative power flows for a line linking nodes $i$ and $j$ represents flows from node $j$ to $i$ and causes the same magnitude of losses as an equal power flow from $i$ to $j$). This absolute value function is linearized such that the flow in the line must be equal to the subtraction of the auxiliary variable for flow in the positive direction, $\Phi^{+}_{l,t}$, and the auxiliary variable for flow in the negative direction, $\Phi^{+}_{l,t}$, of the line. Then, the magnitude of the flow is calculated as the sum of the two auxiliary variables. The sum of positive and negative directional flows are also constrained by the maximum line flow capacity.
```math
\begin{aligned}
% trasmission losses simple
	&\Phi_{l,t} =  \Phi^{+}_{l,t}  - \Phi^{-}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\mid \Phi_{l,t} \mid =  \Phi^{+}_{l,t}  + \Phi^{-}_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\Phi^{+}_{l,t}  + \Phi^{-}_{l,t} \leq \varphi^{max}_{l}, &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}
\end{aligned}
```
If discrete unit commitment decisions are modeled, ``phantom losses'' can be observed wherein the auxiliary variables for flows in both directions ($\Phi^{+}_{l,t}$ and $\Phi^{-}_{l,t}$) are both increased to produce increased losses so as to avoid cycling a thermal generator and incurring start-up costs or opportunity costs related to minimum down times. This unrealistic behavior can be eliminated via inclusion of additional constraints and a set of auxiliary binary variables, $ON^{+}_{l,t} \in {0,1} \forall l \in \mathcal{L}$. Then the following additional constraints are created:
```math
\begin{aligned}
	\Phi^{+}_{l,t} \leq TransON^{+}_{l,t},  &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	\Phi^{-}_{l,t} \leq (\varphi^{max}_{l} + \bigtriangleup\varphi^{max}_{l}) -TransON^{+}_{l,t}, &\quad  \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}
\end{aligned}
```
where $TransON^{+}_{l,t}$ is a continuous variable, representing the product of the binary variable $ON^{+}_{l,t}$ and the expression, $(\varphi^{max}_{l} + \bigtriangleup\varphi^{max}_{l})$. This product cannot be defined explicitly, since it will lead to a bilinear expression involving two variables. Instead, we enforce this definition via the Glover's Linearization as shown below (also referred McCormick Envelopes constraints for bilinear expressions, which is exact when one of the variables is binary).
```math
\begin{aligned}
	TransON^{+}_{l,t} \leq  (\varphi^{max}_{l} + \overline{\bigtriangleup\varphi^{max}_{l}}) \times TransON^{+}_{l,t},  &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T} \\
	TransON^{+}_{l,t} \leq  (\varphi^{max}_{l} + \bigtriangleup\varphi^{max}_{l}),  &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T} \\
	TransON^{+}_{l,t} \leq (\varphi^{max}_{l} + \bigtriangleup\varphi^{max}_{l}) - (\varphi^{max}_{l} + \overline{\bigtriangleup\varphi^{max}_{l}}) \times(1- TransON^{+}_{l,t}),  &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T} \\
\end{aligned}
```
These constraints permit only the positive or negative auxiliary flow variables to be non-zero at a given time period, not both.
For the third option, losses are calculated as a piecewise-linear approximation of a quadratic function of power flows. In order to do this, we represent the absolute value of the line flow variable by the sum of positive stepwise flow variables $(\mathcal{S}^{+}_{m,l,t}, \mathcal{S}^{-}_{m,l,t})$, associated with each partition of line losses computed using the corresponding linear expressions. This can be understood as a segmentwise linear fitting (or first order approximation) of the quadratic losses function. The first constraint below computes the losses a the accumulated sum of losses for each linear stepwise segment of the approximated quadratic function, including both positive domain and negative domain segments. A second constraint ensures that the stepwise variables do not exceed the maximum size per segment. The slope and maximum size for each segment are calculated as per the method in \cite{Zhang2013}.
```math
\begin{aligned}
	& \ell_{l,t} = \frac{\varphi^{ohm}_{l}}{(\varphi^{volt}_{l})^2}\bigg( \sum_{m \in \mathcal{M}}( S^{+}_{m,l}\times \mathcal{S}^{+}_{m,l,t} + S^{-}_{m,l}\times \mathcal{S}^{-}_{m,l,t}) \bigg), &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T} \\
	& \text{\quad Where:} \\
	& \quad S^{+}_{m,l} = \frac{2+4 \times \sqrt{2}\times (m-1)}{1+\sqrt{2} \times (2 \times M-1)} (\varphi^{max}_{l} + \overline{\bigtriangleup\varphi^{max}_{l}}) &\quad \forall m \in [1 \colon M], l \in \mathcal{L}  \\
	& \quad S^{-}_{m,l} = \frac{2+4 \times \sqrt{2}\times (m-1)}{1+\sqrt{2} \times (2 \times M-1)} (\varphi^{max}_{l} + \overline{\bigtriangleup\varphi^{max}_{l}}) &\quad \forall m \in [1 \colon M], l \in \mathcal{L}\\
	& \\
	& \mathcal{S}^{+}_{m,l,t}, \mathcal{S}^{-}_{m,l,t} <= \overline{\mathcal{S}_{m,l}} &\quad \forall m \in [1:M], l \in \mathcal{L}, t \in \mathcal{T} \\
	& \text{\quad Where:} \\
	& \quad \overline{S_{l,z}} =  \begin{cases} \frac{(1+\sqrt{2})}{1+\sqrt{2} \times (2 \times M-1)}  (\varphi^{max}_{l} + \overline{\bigtriangleup\varphi^{max}_{l}}) & \text{if~} m = 1 \\
	\frac{2 \times \sqrt{2} }{1+\sqrt{2} \times (2 \times M-1)} (\varphi^{max}_{l} + \overline{\bigtriangleup\varphi^{max}_{l}}) & \text{if~} m > 1 \end{cases}
\end{aligned}
```
Next, a constraint ensures that the sum of auxiliary segment variables ($m \geq 1$) minus the "zero" segment (which allows values to go into the negative domain) from both positive and negative domains must total the actual power flow across the line, and a constraint ensures that the sum of negative and positive flows do not exceed the maximum flow for the line.
```math
\begin{aligned}
	&\sum_{m \in [1:M]} (\mathcal{S}^{+}_{m,l,t}) - \mathcal{S}^{+}_{0,l,t} =  \Phi_{l,t}, &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\sum_{m \in [1:M]} (\mathcal{S}^{-}_{m,l,t}) - \mathcal{S}^{-}_{0,l,t}  =  - \Phi_{l,t}
\end{aligned}
```
As with losses option 2, this segment-wise approximation of a quadratic loss function also permits 'phantom losses' to avoid cycling thermal units when discrete unit commitment decisions are modeled. In this case, the additional constraints below are also added to ensure that auxiliary segments variables do not exceed maximum value per segment and that they are filled in order; i.e., one segment cannot be non-zero unless prior segment is at its maximum value. Binary constraints deal with absolute value of power flow on each line. If the flow is positive, $\mathcal{S}^{+}_{0,l,t}$ must be zero; if flow is negative, $\mathcal{S}^{+}_{0,l,t}$ must be positive and takes on value of the full negative flow, forcing all $\mathcal{S}^{+}_{m,l,t}$ other segments ($m \geq 1$) to be zero. Conversely, if the flow is negative, $\mathcal{S}^{-}_{0,l,t}$ must be zero; if flow is positive, $\mathcal{S}^{-}_{0,l,t}$ must be positive and takes on value of the full positive flow, forcing all $\mathcal{S}^{-}_{m,l,t}$ other segments ($m \geq 1$) to be zero. Requiring segments to fill in sequential order and binary variables to ensure variables reflect the actual direction of power flows are both necessary to eliminate ``phantom losses'' from the solution. These constraints and binary decisions are ommited if the model is fully linear.
```math
\begin{aligned}
	&\mathcal{S}^{+}_{m,l,t} <=    \overline{\mathcal{S}_{m,l}} \times ON^{+}_{m,l,t}, &\quad \forall m \in [1:M], \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\mathcal{S}^{-}_{m,l,t} <=    \overline{\mathcal{S}_{m,l}} \times ON^{-}_{m,l,t},  &\quad \forall m \in[1:M], \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\mathcal{S}^{+}_{m,l,t} \geq ON^{+}_{m+1,l,t} \times \overline{\mathcal{S}_{m,l}}, &\quad \forall m \in [1:M], \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\mathcal{S}^{-}_{m,l,t} \geq ON^{-}_{m+1,l,t} \times \overline{\mathcal{S}_{m,l}} , &\quad \forall m \in [1:M], \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\mathcal{S}^{+}_{0,l,t} \leq \varphi^{max}_{l} \times (1- ON^{+}_{1,l,t}), &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}\\
	&\mathcal{S}^{-}_{0,l,t} \leq \varphi^{max}_{l} \times (1- ON^{-}_{1,l,t}), &\quad \forall l \in \mathcal{L}, \forall t  \in \mathcal{T}
\end{aligned}
```
"""
function transmission(EP::Model, inputs::Dict, UCommit::Int, NetworkExpansion::Int, CapacityReserveMargin::Int)

	println("Transmission Module")

	dfGen = inputs["dfGen"]

	T = inputs["T"]     # Number of time steps (hours)
	Z = inputs["Z"]     # Number of zones
	L = inputs["L"]     # Number of transmission lines
	SEG = inputs["SEG"] # Number of load curtailment segments

	## sets and indices for transmission losses and expansion
	TRANS_LOSS_SEGS = inputs["TRANS_LOSS_SEGS"] # Number of segments used in piecewise linear approximations quadratic loss functions - can only take values of TRANS_LOSS_SEGS =1, 2
	LOSS_LINES = inputs["LOSS_LINES"] # Lines for which loss coefficients apply (are non-zero);
	if NetworkExpansion == 1
		# Network lines and zones that are expandable have non-negative maximum reinforcement inputs
		EXPANSION_LINES = inputs["EXPANSION_LINES"]
		NO_EXPANSION_LINES = inputs["NO_EXPANSION_LINES"]
	end

	### Variables ###

	# Power flow on each transmission line "l" at hour "t"
	@variable(EP, vFLOW[l=1:L,t=1:T]);

	if NetworkExpansion == 1
		# Transmission network capacity reinforcements per line
		@variable(EP, vNEW_TRANS_CAP[l in EXPANSION_LINES] >= 0)
	end

  	if (TRANS_LOSS_SEGS==1)  #loss is a constant times absolute value of power flow
		# Positive and negative flow variables
		@variable(EP, vTAUX_NEG[l in LOSS_LINES,t=1:T] >= 0)
		@variable(EP, vTAUX_POS[l in LOSS_LINES,t=1:T] >= 0)

		if UCommit == 1
			# Single binary variable to ensure positive or negative flows only
			@variable(EP, vTAUX_POS_ON[l in LOSS_LINES,t=1:T],Bin)
			# Continuous variable representing product of binary variable (vTAUX_POS_ON) and avail transmission capacity
			@variable(EP, vPROD_TRANSCAP_ON[l in LOSS_LINES,t=1:T]>=0)
		end
	else # TRANS_LOSS_SEGS>1
		# Auxiliary variables for linear piecewise interpolation of quadratic losses
		@variable(EP, vTAUX_NEG[l in LOSS_LINES, s=0:TRANS_LOSS_SEGS, t=1:T] >= 0)
		@variable(EP, vTAUX_POS[l in LOSS_LINES, s=0:TRANS_LOSS_SEGS, t=1:T] >= 0)
		if UCommit == 1
			# Binary auxilary variables for each segment >1 to ensure segments fill in order
			@variable(EP, vTAUX_POS_ON[l in LOSS_LINES, s=1:TRANS_LOSS_SEGS, t=1:T], Bin)
			@variable(EP, vTAUX_NEG_ON[l in LOSS_LINES, s=1:TRANS_LOSS_SEGS, t=1:T], Bin)
		end
    	end

	# Transmission losses on each transmission line "l" at hour "t"
	@variable(EP, vTLOSS[l in LOSS_LINES,t=1:T] >= 0)

	### Expressions ###

	## Transmission power flow and loss related expressions:
	# Total availabile maximum transmission capacity is the sum of existing maximum transmission capacity plus new transmission capacity
	if NetworkExpansion == 1
		@expression(EP, eAvail_Trans_Cap[l=1:L],
			if l in EXPANSION_LINES
				inputs["pTrans_Max"][l] + vNEW_TRANS_CAP[l]
			else
				inputs["pTrans_Max"][l] + EP[:vZERO]
			end
		)
	else
		@expression(EP, eAvail_Trans_Cap[l=1:L], inputs["pTrans_Max"][l] + EP[:vZERO])
	end

	# Net power flow outgoing from zone "z" at hour "t" in MW
    	@expression(EP, eNet_Export_Flows[z=1:Z,t=1:T], sum(inputs["pNet_Map"][l,z] * vFLOW[l,t] for l=1:L))

	# Losses from power flows into or out of zone "z" in MW
    	@expression(EP, eLosses_By_Zone[z=1:Z,t=1:T], sum(abs(inputs["pNet_Map"][l,z]) * vTLOSS[l,t] for l in LOSS_LINES))

	## Objective Function Expressions ##

	if NetworkExpansion == 1
		@expression(EP, eTotalCNetworkExp, sum(vNEW_TRANS_CAP[l]*inputs["pC_Line_Reinforcement"][l] for l in EXPANSION_LINES))
		EP[:eObj] += eTotalCNetworkExp
    	end

	## End Objective Function Expressions ##

	## Power Balance Expressions ##

	@expression(EP, ePowerBalanceNetExportFlows[t=1:T, z=1:Z],
		-eNet_Export_Flows[z,t])
	@expression(EP, ePowerBalanceLossesByZone[t=1:T, z=1:Z],
		-(1/2)*eLosses_By_Zone[z,t])

	EP[:ePowerBalance] += ePowerBalanceLossesByZone
	EP[:ePowerBalance] += ePowerBalanceNetExportFlows

	# Capacity Reserves Margin policy
	if CapacityReserveMargin > 0
		if Z > 1
			@expression(EP, eCapResMarBalanceTrans[res=1:inputs["NCapacityReserveMargin"], t=1:T], sum(inputs["dfTransCapRes_excl"][l,res] * inputs["dfDerateTransCapRes"][l,res]* EP[:vFLOW][l,t] for l in 1:L))
			EP[:eCapResMarBalance] -= eCapResMarBalanceTrans
		end
	end

	### Constraints ###

  	## Power flow and transmission (between zone) loss related constraints

	# Maximum power flows, power flow on each transmission line cannot exceed maximum capacity of the line at any hour "t"
	# Allow expansion of transmission capacity for lines eligible for reinforcement
	@constraints(EP, begin
		cMaxFlow_out[l=1:L, t=1:T], vFLOW[l,t] <= eAvail_Trans_Cap[l]
		cMaxFlow_in[l=1:L, t=1:T], vFLOW[l,t] >= -eAvail_Trans_Cap[l]
	end)

	# If network expansion is used:
	if NetworkExpansion == 1
		# Transmission network related power flow and capacity constraints
		# Constrain maximum line capacity reinforcement for lines eligible for expansion
		@constraint(EP, cMaxLineReinforcement[l in EXPANSION_LINES], vNEW_TRANS_CAP[l] <= inputs["pMax_Line_Reinforcement"][l])
	end
	#END network expansion contraints

	# Transmission loss related constraints - linear losses as a function of absolute value
	if TRANS_LOSS_SEGS == 1

		@constraints(EP, begin
			# Losses are alpha times absolute values
			cTLoss[l in LOSS_LINES, t=1:T], vTLOSS[l,t] == inputs["pPercent_Loss"][l]*(vTAUX_POS[l,t]+vTAUX_NEG[l,t])

			# Power flow is sum of positive and negative components
			cTAuxSum[l in LOSS_LINES, t=1:T], vTAUX_POS[l,t]-vTAUX_NEG[l,t] == vFLOW[l,t]

			# Sum of auxiliary flow variables in either direction cannot exceed maximum line flow capacity
			cTAuxLimit[l in LOSS_LINES, t=1:T], vTAUX_POS[l,t]+vTAUX_NEG[l,t] <= eAvail_Trans_Cap[l]
		end)

		if UCommit == 1
			# Constraints to limit phantom losses that can occur to avoid discrete cycling costs/opportunity costs due to min down
			@constraints(EP, begin
				cTAuxPosUB[l in LOSS_LINES, t=1:T], vTAUX_POS[l,t] <= vPROD_TRANSCAP_ON[l,t]

				# Either negative or positive flows are activated, not both
				cTAuxNegUB[l in LOSS_LINES, t=1:T], vTAUX_NEG[l,t] <= eAvail_Trans_Cap[l]-vPROD_TRANSCAP_ON[l,t]

				# McCormick representation of product of continuous and binary variable
				# (in this case, of: vPROD_TRANSCAP_ON[l,t] = eAvail_Trans_Cap[l] * vTAUX_POS_ON[l,t])
				# McCormick constraint 1
				[l in LOSS_LINES,t=1:T], vPROD_TRANSCAP_ON[l,t] <= inputs["pTrans_Max_Possible"][l]*vTAUX_POS_ON[l,t]

				# McCormick constraint 2
				[l in LOSS_LINES,t=1:T], vPROD_TRANSCAP_ON[l,t] <= eAvail_Trans_Cap[l]

				# McCormick constraint 3
				[l in LOSS_LINES,t=1:T], vPROD_TRANSCAP_ON[l,t] >= eAvail_Trans_Cap[l]-(1-vTAUX_POS_ON[l,t])*inputs["pTrans_Max_Possible"][l]
			end)
		end

	end # End if(TRANS_LOSS_SEGS == 1) block

	# When number of segments is greater than 1
	if (TRANS_LOSS_SEGS > 1)
		## between zone transmission loss constraints
		# Losses are expressed as a piecewise approximation of a quadratic function of power flows across each line
		# Eq 1: Total losses are function of loss coefficient times the sum of auxilary segment variables across all segments of piecewise approximation
		# (Includes both positive domain and negative domain segments)
		@constraint(EP, cTLoss[l in LOSS_LINES, t=1:T], vTLOSS[l,t] ==
							(inputs["pTrans_Loss_Coef"][l]*sum((2*s-1)*(inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)*vTAUX_POS[l,s,t] for s=1:TRANS_LOSS_SEGS)) +
							(inputs["pTrans_Loss_Coef"][l]*sum((2*s-1)*(inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)*vTAUX_NEG[l,s,t] for s=1:TRANS_LOSS_SEGS)) )
		# Eq 2: Sum of auxilary segment variables (s >= 1) minus the "zero" segment (which allows values to go negative)
		# from both positive and negative domains must total the actual power flow across the line
		@constraints(EP, begin
			cTAuxSumPos[l in LOSS_LINES, t=1:T], sum(vTAUX_POS[l,s,t] for s=1:TRANS_LOSS_SEGS)-vTAUX_POS[l,0,t]  == vFLOW[l,t]
			cTAuxSumNeg[l in LOSS_LINES, t=1:T], sum(vTAUX_NEG[l,s,t] for s=1:TRANS_LOSS_SEGS) - vTAUX_NEG[l,0,t]  == -vFLOW[l,t]
		end)
		if UCommit == 0 || UCommit == 2
			# Eq 3: Each auxilary segment variables (s >= 1) must be less than the maximum power flow in the zone / number of segments
			@constraints(EP, begin
				cTAuxMaxPos[l in LOSS_LINES, s=1:TRANS_LOSS_SEGS, t=1:T], vTAUX_POS[l,s,t] <= (inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)
				cTAuxMaxNeg[l in LOSS_LINES, s=1:TRANS_LOSS_SEGS, t=1:T], vTAUX_NEG[l,s,t] <= (inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)
			end)
		else # Constraints that can be ommitted if problem is convex (i.e. if not using MILP unit commitment constraints)
			# Eqs 3-4: Ensure that auxilary segment variables do not exceed maximum value per segment and that they
			# "fill" in order: i.e. one segment cannot be non-zero unless prior segment is at it's maximum value
			# (These constraints are necessary to prevents phantom losses in MILP problems)
			@constraints(EP, begin
				cTAuxOrderPos1[l in LOSS_LINES, s=1:TRANS_LOSS_SEGS, t=1:T], vTAUX_POS[l,s,t] <=  (inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)*vTAUX_POS_ON[l,s,t]
				cTAuxOrderNeg1[l in LOSS_LINES, s=1:TRANS_LOSS_SEGS, t=1:T], vTAUX_NEG[l,s,t] <=  (inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)*vTAUX_NEG_ON[l,s,t]
				cTAuxOrderPos2[l in LOSS_LINES, s=1:(TRANS_LOSS_SEGS-1), t=1:T], vTAUX_POS[l,s,t] >=  (inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)*vTAUX_POS_ON[l,s+1,t]
				cTAuxOrderNeg2[l in LOSS_LINES, s=1:(TRANS_LOSS_SEGS-1), t=1:T], vTAUX_NEG[l,s,t] >=  (inputs["pTrans_Max_Possible"][l]/TRANS_LOSS_SEGS)*vTAUX_NEG_ON[l,s+1,t]
			end)

			# Eq 5: Binary constraints to deal with absolute value of vFLOW.
			@constraints(EP, begin
				# If flow is positive, vTAUX_POS segment 0 must be zero; If flow is negative, vTAUX_POS segment 0 must be positive
				# (and takes on value of the full negative flow), forcing all vTAUX_POS other segments (s>=1) to be zero
				cTAuxSegmentZeroPos[l in LOSS_LINES, t=1:T], vTAUX_POS[l,0,t] <= inputs["pTrans_Max_Possible"][l]*(1-vTAUX_POS_ON[l,1,t])

				# If flow is negative, vTAUX_NEG segment 0 must be zero; If flow is positive, vTAUX_NEG segment 0 must be positive
				# (and takes on value of the full positive flow), forcing all other vTAUX_NEG segments (s>=1) to be zero
				cTAuxSegmentZeroNeg[l in LOSS_LINES, t=1:T], vTAUX_NEG[l,0,t] <= inputs["pTrans_Max_Possible"][l]*(1-vTAUX_NEG_ON[l,1,t])
			end)
		end
	end # End if(TRANS_LOSS_SEGS > 0) block

	return EP
end

ext: jl

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
export PopObj, PopData, PopDataInfo, show, Genotype, GenoArray, SNP, MSat
export getindex, getproperty

"""
    AbstractType PopObj
Generic AbstractType for use in PopGen.jl
"""
abstract type PopObj end

"""
```
mutable struct PopDataInfo
    samples::Int64
    sampeinfo::DataFrame
    loci::Int64
    locusinfo::DataFrame
    populations::Int64
    ploidy::Union{Int8, Vector{Int8}}
    biallelic::Bool
end
```
The data struct used internally as `PopData.metadata` fields to store basic information
about the `PopData` for easy access.
"""
mutable struct PopDataInfo
    samples::Int64
    sampleinfo::DataFrame
    loci::Int64
    locusinfo::DataFrame
    populations::Int64
    ploidy::Union{Int8, Vector{Int8}}
    biallelic::Bool
end

# constructor FORMAT just the genodata dataframe
function PopDataInfo(genodf::DataFrame)
    sampleinfo = unique(dropmissing(genodf), :name)
    sampleinfo.ploidy = [ismissing(geno) ? Int8(0) : Int8(length(geno)) for geno in sampleinfo.genotype]
    select!(sampleinfo, :name => collect => :name, :population, :ploidy)
    ploidy = unique(sampleinfo.ploidy)
    ploidy = length(ploidy) == 1 ? ploidy[1] : ploidy
    PopDataInfo(
        length(genodf.name.pool),
        sampleinfo,
        length(genodf.locus.pool),
        DataFrame(:chromosome => PooledArray(fill("0",length(genodf.locus.pool)) , compress = true), :locus => genodf.locus.pool, :cm => Int8(0), :bp => Int8(0)),
        length(genodf.population.pool),
        ploidy,
        isbiallelic(genodf)
    )
end

function Base.show(io::IO, data::PopDataInfo)
    println(io, " ploidy:       ", join(data.ploidy, ","))
    println(io, " loci:         ", data.loci)
    println(io, " samples:      ", data.samples)
    println(io, " populations:  ", data.populations)
    println(io, " biallelic:    ", data.biallelic)
    println(io, " other fields: locusinfo, sampleinfo")
end

"""
```
PopData
    metadata::PopDataInfo
    genodata::DataFrame
```
The data struct used for the PopGen population genetics ecosystem. You are
**strongly** discouraged from manually creating tables to pass into a PopData,
and instead should use the provided file importers and utilities.

    - `metadata` PopDataInfo of  data information
        - `samples` - the number of samples in the data
        - `sampleinfo` - DataFrame of sample names,populations, ploidy, etc.
        - `loci` - the number of loci in the data
        - `locusinfo` - DataFrame of locus names, chromosome, physical position, etc.
        - `populations` - the number of populations in the data
        - `ploidy` - the ploidy (or ploidies) present in the data
        - `biallelic` - if all the markers are biallelic
    - `genodata` DataFrame of sample genotype records
        - `name` - the individual/sample names [`PooledArray`]
        - `population` - population names [`PooledArray`]
        - `locus` - locus names [`PooledArray`]
        - `genotype` - genotype values [`NTuple{N,Signed}`]

"""
struct PopData <: PopObj
    metadata::PopDataInfo
    genodata::DataFrame
    function PopData(meta::PopDataInfo, loci::DataFrame)
        if "ploidy" ∉ names(meta.sampleinfo)
            insertcols!(meta.sampleinfo, 3, :ploidy => Int8(0))
        end
        metacols = ["name", "population", "ploidy"]
        metacheck = intersect(metacols, names(meta.sampleinfo))
        if metacheck != metacols
            throw(error("metadata sampleinfo missing columns $(symdiff(metacheck, metacols))"))
        end
        genocols = ["name", "population", "locus", "genotype"]
        genocheck = intersect(genocols, names(loci))
        if genocheck != genocols
            throw(error("genodata missing columns $(symdiff(genocheck, genocols))"))
        end
        sort(meta.sampleinfo.name) != sort(loci.name.pool) && throw(ArgumentError("metadata and genodata do not contain the same sample names"))
        sort(unique(meta.sampleinfo.population)) != sort(loci.population.pool) && throw(ArgumentError("metadata and genodata do not contain the same population names"))
        new(meta, loci)
    end
end

PopData(data::DataFrame) = PopData(PopDataInfo(data), data)

# method to update PopDataInfo from PopData, all in one swoop
function PopDataInfo!(data::PopData)
    data.metadata.samples = length(data.genodata.name.pool)
    data.metadata.loci = length(data.genodata.locus.pool)
    data.metadata.populations = length(data.genodata.population.pool)
    filter!(:name => x -> x ∈ data.genodata.name.pool,  data.sampleinfo)
    filter!(:locus => x -> x ∈ data.genodata.locus.pool,  data.locusinfo)
    if "ploidy" ∈ names(data.sampleinfo)
        ploidy = unique(data.sampleinfo.ploidy)
        ploidy = length(ploidy) == 1 ? Int8(ploidy[1]) : Int8.(ploidy)
    else
        ploidy = Int8(0)
    end
    data.metadata.ploidy = ploidy
    data.metadata.biallelic = data.metadata.biallelic ? true : isbiallelic(data)
    return
end

# method to update preexisting PopDataInfo with new genodata
# useful for getindex and creating new PopData from that
function PopDataInfo!(popdatainfo::PopDataInfo, genodata::DataFrame)
    popdatainfo.samples = length(genodata.name.pool)
    popdatainfo.loci = length(genodata.locus.pool)
    popdatainfo.populations = length(genodata.population.pool)
    filter!(:name => x -> x ∈ genodata.name.pool,  popdatainfo.sampleinfo)
    filter!(:locus => x -> x ∈ genodata.locus.pool,  popdatainfo.locusinfo)
    if "ploidy" ∈ names(popdatainfo.sampleinfo)
        ploidy = unique(popdatainfo.sampleinfo.ploidy)
        ploidy = length(ploidy) == 1 ? Int8(ploidy[1]) : Int8.(ploidy)
    else
        ploidy = Int8(0)
    end
    popdatainfo.ploidy = ploidy
    popdatainfo.biallelic = popdatainfo.biallelic ? true : isbiallelic(genodata)
    return popdatainfo
end

"""
    Genotype::DataType
For convenience purposes, an alias for `NTuple{N, <:Signed} where N`, which is
the type describing individual genotypes in PopData. Specifically, there exist
`SNP` as an alias for `NTuple{N, Int8}` and `MSat` for `NTuple{N, Int16}`
"""
const Genotype = NTuple{N, <:Signed} where N

"""
    SNP::DataType
An alias for `NTuple{N, Int8}`
"""
const SNP = NTuple{N, Int8} where N
_SNP(geno) = all(ismissing.(geno)) ? missing : SNP(geno)

"""
MSat::DataType
An alias for `NTuple{N, Int16}`
    """
const MSat = NTuple{N, Int16} where N
_MSat(geno) = all(ismissing.(geno)) ? missing : MSat(geno)

"""
    GenoArray::DataType
An alias for an `AbstractVector` of elements `Missing`
and `Genotype`, which itself is of type `NTuple{N, <:Integer} where N`.
The definition as an `AbstractVector` adds flexibility for `SubArray`
cases.
"""
const GenoArray = AbstractVector{S} where S<:Union{Missing,Genotype}
#const SnpArray = PooledVector{Union{Missing, Tuple{Int8, Int8}}, UInt8, Vector{UInt8}}

function _ploidy2text(ploidy::Int8)
    ploidy == 0 ? "" :
    ploidy == 1 ? "Haploid" :
    ploidy == 2 ? "Diploid" :
    ploidy == 3 ? "Triploid" :
    ploidy == 4 ? "Tetraploid" :
    ploidy == 5 ? "Pentaploid" :
    ploidy == 6 ? "Hexaploid" : "Octaploid"
end

_ploidy2text(ploidy::Vector{Int8}) = "Mixed-ploidy"

function Base.show(io::IO, data::PopData)
    if occursin("Int16", string(eltype(data.genodata.genotype)))
        marker = "Microsatellite"
    else
        marker = "SNP"
    end
    ploidy = _ploidy2text(data.metadata.ploidy)
    ploidy = ploidy == "" ? "" : ploidy * ", "
    n_loc = data.metadata.loci
    println(io, "PopData", "{" * ploidy, n_loc, " " , marker * " loci}")
    println(io, "  Samples: $(data.metadata.samples)")
    print(io, "  Populations: $(data.metadata.populations)")
    allcols = vcat(names(data.sampleinfo), names(data.genodata)) |> unique
    extracols = symdiff(allcols, ["name", "population", "ploidy", "locus", "genotype"])
    if !isempty(extracols)
        print(io, "\n  Other Info: ", extracols)
    end
end

function Base.getindex(data::PopData, idx::Symbol)
    if idx == :sampleinfo
        return data.sampleinfo
    elseif idx == :genodata
        return data.genodata
    elseif idx == :name
        return data.metadata..sampleinfo.name
    elseif idx == :population
        return data.sampleinfo.population
    elseif idx == :ploidy
        return data.sampleinfo.ploidy
    elseif idx == :locus
        return data.genodata.locus.pool
    elseif idx == :coordinates
        return data.sampleinfo[:, [:longitude, :latitude]]
    else
        throw(ArgumentError("Cannot directly index PopData with the \':$idx\' field"))
    end
end

function Base.getindex(data::PopData, args)
    geno = getindex(data.genodata, args, :)
    transform!(
        geno,
        1 => (i -> PooledArray(i, compress = true)) => :name,
        2 => (i -> PooledArray(i, compress = true)) => :population,
        3 => (i -> PooledArray(i, compress = true)) => :locus,
        4
    )
    pdinfo = deepcopy(data.info)
    out = PopDataInfo!(pdinfo, geno)
    PopData(out, geno)
end

function Base.getindex(data::PopData, expression, cols)
    getindex(data.genodata, expression, cols)
end

# shortcut methods for convenience and less verbose typing
"""
    getproperty(data::PopData, field::Symbol)
    data.field

A convenience method to access certain elements in a `PopData` with fewer keystrokes.
Essentially a standard `getproperty` call, except `sampleinfo` accesses `metadata.sampleinfo`,
`locusinfo` accesses `metadata.locusinfo`, and `info` is an alias for `metadata`.

## Example
```julia
cats = @nancycats ;
cats.metadata == cats.info
cats.metadata.sampleinfo == cats.sampleinfo
cats.metadata.locusinfo == cats.locusinfo
```
"""
function Base.getproperty(data::PopData, field::Symbol)
    if field == :sampleinfo
        return data.metadata.sampleinfo
    elseif field == :locusinfo
        return data.metadata.locusinfo
    elseif field == :info
        return data.metadata
    else
        return getfield(data, field)
    end
end


text:
module ALFFA

import ..LABELS
import ..LOCATIONS
import ..NON_SPEECH_UNITS

using Glob
using NaturalSort

const URL_REPO = "https://github.com/getalp/ALFFA_PUBLIC"
const LANGUAGES = Set([:amharic, :swahili, :wolof])
const WAVPATTERN = Dict(
    :amharic => "wav/*wav",
    :swahili => "wav/*/*wav",
)
const SUBDIR = Dict(
    :amharic => joinpath("ASR", "AMHARIC"),
    :swahili=> joinpath("ASR", "SWAHILI"),
)

function prepare(datadir, lang)
    lang ∈ LANGUAGES || error("$lang is unknown, see ALFFA.LANGUAGES for supported language")

    localdir = mkpath(joinpath(datadir, LOCATIONS[:local]))

    repopath = joinpath(localdir, "alffa-github")
    if ! ispath(repopath)
        @info "downloading the data, this may take a while..."
        run(`git clone $URL_REPO $repopath`)
    else
        @info "data aleady downloaded"
    end

    for dataset in ["train", "test"]
        datasetdir = mkpath(joinpath(datadir, "$lang", dataset))
        @info "preparing $datasetdir..."

        uttids = Set()

        dir = joinpath(repopath, SUBDIR[lang], "data", dataset)
        open(joinpath(datasetdir, LOCATIONS[:wavs]), "w") do f

            for path in glob(WAVPATTERN[lang], dir)
                uttid = splitext(basename(path))[1]
                push!(uttids, uttid)
                println(f, uttid, "\t", path)
            end
        end

        open(joinpath(datasetdir, LOCATIONS[:uttids]), "w") do f
            for uttid in sort(collect(uttids), lt=natural)
                println(f, uttid)
            end
        end

        cp(joinpath(dir, "text"), joinpath(datasetdir, LOCATIONS[:trans]*".wrd"), force = true)
        cp(joinpath(dir, "utt2spk"), joinpath(datasetdir, LOCATIONS[:utt2spk]), force = true)
    end

    langdir = mkpath(joinpath(datadir, "$lang", LOCATIONS[:lang]))
    src = joinpath(repopath, SUBDIR[lang], "lang", "lexicon.txt")
    dest = joinpath(langdir, LOCATIONS[:lexicon])
    run(pipeline(`cat $src`, `grep -v SIL`, `grep -v unk`, `grep -v UNK`, dest))

    open(dest, "a") do f
        println(f, LABELS[:sil], "\t", LABELS[:nonspeechunit])
    end

    words = joinpath(langdir, LOCATIONS[:words])
    run(pipeline(`cut -f1 $dest`, "$words"))

    open(joinpath(langdir, LOCATIONS[:units]), "w") do f
        println(f, LABELS[:sil], "\t", LABELS[:nonspeechunit])
        open(joinpath(repopath, SUBDIR[lang], "lang", "nonsilence_phones.txt"), "r") do f2
            for line in eachline(f2)
                println(f, line, "\t", LABELS[:speechunit])
            end
        end
    end

    return nothing
end

end # module

max_issues_repo_name: lucasondel/SpeechDatasets

max_forks_repo_name: lucasondel/SpeechDatasets


text:
function leap_it()

    # Prompt the user for input and convert it to an int from a string
    year = parse(Int64, (Base.prompt("Which year would you like to check if it's a leap year?")))

    # If the remainder of dividing the year by 4 is zero (which is the basic math to check a leap year)
    if (year % 4) == 0
        println("$year is a leap year...")
    else

        # Calculate how long till the next leap year, if the entered year is not one
        years_away = 4 - (year % 4)
        println("$year is not a leap year, you have to wait $years_away year...")
    end
end


text:
using JuMPTutorials
using Test
using Statistics

cd("introduction")
@testset "Introduction" begin
include("introduction/an_introduction_to_julia.jl")
include("introduction/getting_started_with_JuMP.jl")
@test objective_value(model) ≈ 205
include("introduction/variables_constraints_objective.jl")
@test objective_value(vector_model) ≈ 4.9230769230769225
include("introduction/solvers_and_solutions.jl")
end
cd("..")

cd("using_JuMP")
@testset "Using JuMP" begin
include("using_JuMP/working_with_data_files.jl")
@test objective_value(model) == 23
@test countryindex == [1, 5, 9, 10, 38, 39, 55, 63, 64, 75, 78, 81, 89, 104, 107, 130, 138, 158, 162, 167, 182, 188, 190]
include("using_JuMP/problem_modification.jl")
include("using_JuMP/nonlinear_modelling.jl")
end
cd("..")

cd("optimization_concepts")
@testset "Optimization Concepts" begin
include("optimization_concepts/benders_decomposition.jl")
@test objective_value(master_problem_model) == -4.0
@test value.(x) == [0.0, 1.0]
include("optimization_concepts/integer_programming.jl")
include("optimization_concepts/conic_programming.jl")
end
cd("..")

cd("modelling")
@testset "Modelling Examples" begin
include("modelling/sudoku.jl")
@test sol == [
    5  3  4  6  7  8  9  1  2;
    6  7  2  1  9  5  3  4  8;
    1  9  8  3  4  2  5  6  7;
    8  5  9  7  6  1  4  2  3;
    4  2  6  8  5  3  7  9  1;
    7  1  3  9  2  4  8  5  6;
    9  6  1  5  3  7  2  8  4;
    2  8  7  4  1  9  6  3  5;
    3  4  5  2  8  6  1  7  9]

include("modelling/problems_on_graphs.jl")
@test value.(y) == [0.0, 1.0, 1.0, 1.0, 0.0, 0.0]
@test value.(x) == [0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 1.0]
@test value.(m) == [
    0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0;
    0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0;
    0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0;
    0.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0;
    1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
    0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0;
    0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0;
    0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0
]
@test value(sum(z)) ≈ 3
@test sum(c) ≈ nv(g)

include("modelling/network_flows.jl")
@test objective_value(shortest_path) == 55.0
@test value.(x) == [
    0.0 0.0 1.0 0.0 0.0;
    0.0 0.0 0.0 0.0 0.0;
    0.0 0.0 0.0 1.0 0.0;
    0.0 1.0 0.0 0.0 0.0;
    0.0 0.0 0.0 0.0 0.0]
@test objective_value(assignment) ≈ 20.0
@test value.(y) == [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 1.0 0.0 0.0 0.0; 0.0 0.0 0.0 1.0]
@test objective_value(max_flow) ≈ 6.0

include("modelling/finance.jl")
@test objective_value(financing) ≈ 92.49694915254233
@test objective_value(auction) ≈ 21.0
@test value.(y) == [1.0, 1.0, 1.0, 0.0, 0.0, 0.0]
@test objective_value(portfolio) ≈ 22634.41784988414
@test value.(x) ≈ [497.0455298498642, 0.0, 502.95448015948074]
include("modelling/power_systems.jl")
@test g_opt == [1000.0, 300.0]
@test w_opt == 200.0
@test w_f - w_opt == 0
@test obj == 90000.0
include("modelling/geometric_problems.jl")
@test value.(p) ≈ [
     0.44790964261631827    0.0468981793661497;
    -0.03193526635198919   -0.6706136210384356;
     0.404335805799056     -0.45130815913688466;
    -0.39429534726904925   -0.13282535401213758;
     0.02532703978422118    0.4124207687120701;
    -0.0016520566419420052 -0.43954821308159137;
     1.0         1.0;
     1.0        -1.0;
    -1.0        -1.0;
    -1.0         1.0;
     1.0        -0.5;
    -1.0        -0.2;
    -0.2        -1.0;
     0.1         1.0]
include("modelling/rocket_control.jl")
@test objective_value(rocket) ≈ 1.0128340648308058
end
cd("..")


text:
export replace_all

"""
    replace_all(path::AbstractString)
    replace_all(file::File)
    replace_all(dir::Dir)

Replace headers, bullets, blurbs, math inline and blocks and footnotes.
"""
replace_all(path::AbstractString) = replace_all(Path(path))

function replace_all(file::File)
    replace_lines(file)
    replace_blurbs(file)
    replace_math(file)
    replace_codeblocks(file)
    replace_footnotes(file)
end

replace_all(dir::Dir) = visitdir(replace_all, dir.path)

max_stars_repo_name: ordovician/AsciiDocTools

max_issues_repo_name: ordovician/AsciiDocTools

max_forks_repo_name: ordovician/AsciiDocTools


text:
export elu, gelu, hardsigmoid, hardswish, relu, selu, sigm, swish, tanh_
using AutoGrad: AutoGrad, @primitive1, @zerograd
using SpecialFunctions: erf

# AutoGrad does not support broadcasted functions with keywords (https://github.com/denizyuret/AutoGrad.jl/issues/124)
# In Ops21 we will discourage use of broadcasted notation for activations and use f(x::Array) instead of f.(x::Array)

"""
    elu(x)

Return `(x > 0 ? x : exp(x)-1)`.

Reference: Fast and Accurate Deep Network Learning by Exponential Linear Units (ELUs) (https://arxiv.org/abs/1511.07289).

"""
elu(x::T) where {T<:Real} = (x >= 0 ? x : exp(x)-T(1))
eluback(x::T,y::T,dy::T) where {T<:Real} = (y >= 0 ? dy : dy * (T(1)+y))

elu(x::T) where {T<:Array} = elu.(x)
eluback(x::T, y::T, dy::T) where {T<:Array} = eluback.(x, y, dy)

@primitive1 elu(x),dy,y  eluback(x,y,dy)
#@primitive1 eluback(x::Real,y::Real,dy::Real),ddx  (ddx.*((y.>=0).+(y.<0).*(y.+1)))  (ddx.*dy.*(y.<0))
@zerograd eluback(x,y,dy)       # TODO: second derivatives

"""
    gelu(x)

Return `x * Φ(x)` where `Φ(x)` is the cumulative distribution function for the Gaussian
distribution.

References:
* [Hendrycks and Gimpel 2016](https://arxiv.org/abs/1606.08415) Gaussian Error Linear Units (GELUs)

"""
gelu(x::T) where {T<:Real} = x * (1 + erf(x/T(sqrt(2)))) / 2
geluback(x::T,y::T,dy::T) where {T<:Real} = dy * ((1 + erf(x/T(sqrt(2)))) / 2 + x*exp(-x^2/2)/T(sqrt(2pi)))

gelu(x::T) where {T<:Array} = gelu.(x)
geluback(x::T, y::T, dy::T) where {T<:Array} = geluback.(x, y, dy)

@primitive1 gelu(x),dy,y  geluback(x,y,dy)
# @primitive1 geluback(x::Real,y::Real,dy::Real),ddx,dx  (x2=x.^2; ddx.*dy.*(exp.(-x2 ./ 2).*(2 .- x2)./eltype(x)(sqrt(2pi))))  (ddx .* dx ./ dy)
@zerograd geluback(x,y,dy)

# Other possible approximations given in the paper:
# gelu(x::T) where T = 0.5*x*(1 + tanh(T(sqrt(2/pi))*(x + T(0.044715)*x^3)))
# gelu(x::T) where T = x*sigm(T(1.702)*x)

"""
    hardsigmoid(x)

Return `(x <= -3 ? 0 : x >= 3 ? 1 : (x+3)/6)`.

References:
* [Howard+ 2019](https://arxiv.org/abs/1905.02244) Searching for MobileNetV3

"""
hardsigmoid(x::T) where {T<:Real} = (x <= T(-3) ? T(0) : x >= T(3) ? T(1) : (x+3)/6)
hardsigmoidback(x::T,y::T,dy::T) where {T<:Real} = dy * (x <= T(-3) ? T(0) : x >= T(3) ? T(0) : T(1/6))

hardsigmoid(x::T) where {T<:Array} = hardsigmoid.(x)
hardsigmoidback(x::T, y::T, dy::T) where {T<:Array} = hardsigmoidback.(x, y, dy)

@primitive1 hardsigmoid(x),dy,y  hardsigmoidback(x,y,dy)
# @primitive1 hardsigmoidback(x::Real,y::Real,dy::Real),ddx,dx  ??
@zerograd hardsigmoidback(x,y,dy)

"""
    hardswish(x)

Return `(x <= -3 ? 0 : x >= 3 ? x : x*(x+3)/6)`.

References:
* [Howard+ 2019](https://arxiv.org/abs/1905.02244) Searching for MobileNetV3

"""
hardswish(x::T) where {T<:Real} = (x <= T(-3) ? T(0) : x >= T(3) ? x : x*(x+3)/6)
hardswishback(x::T,y::T,dy::T) where {T<:Real} = dy * (x <= T(-3) ? T(0) : x >= T(3) ? T(1) : (2x+3)/6)

hardswish(x::T) where {T<:Array} = hardswish.(x)
hardswishback(x::T, y::T, dy::T) where {T<:Array} = hardswishback.(x, y, dy)

@primitive1 hardswish(x),dy,y  hardswishback(x,y,dy)
# @primitive1 hardswishback(x::Real,y::Real,dy::Real),ddx,dx  ??
@zerograd hardswishback(x,y,dy)

"""
    relu(x; max_value=Inf, negative_slope=0, threshold=0)

Return the following:

    (x >= max_value ? max_value :
     x >= threshold ? x :
     negative_slope * (x - threshold))
"""
function relu(x::T; max_value=Inf, negative_slope=0, threshold=0) where {T<:Real}
    (x >= max_value ? oftype(x, max_value) :
     x >= threshold ? x :
     negative_slope == 0 ? zero(T) :
     negative_slope * (x - oftype(x, threshold)))
end

function reluback(x::T,y::T,dy::T; max_value=Inf, negative_slope=0, threshold=0) where {T<:Real}
    dy * (x >= max_value ? zero(T) :
          x >= threshold ? one(T) :
          oftype(x, negative_slope))
end

relu(x::T; o...) where {T<:Array} = relu.(x; o...)
reluback(x::T, y::T, dy::T; o...) where {T<:Array} = reluback.(x, y, dy; o...)

@primitive1 relu(x; o...),dy,y   reluback(x,y,dy; o...)
#@primitive1 reluback(x::Real,y::Real,dy::Real; o...),ddx,dx  zero(dy)
@zerograd reluback(x,y,dy; o...)

"""
    selu(x)

Return `λ01 * (x > 0 ? x : α01 * (exp(x)-1))` where `λ01=1.0507009873554805` and `α01=1.6732632423543778`.

Reference: Self-Normalizing Neural Networks (https://arxiv.org/abs/1706.02515).
"""
selu(x::T) where {T<:Real} = (x >= 0 ? T(λ01)*x : T(λα01)*(exp(x)-T(1)))
seluback(x::T,y::T,dy::T) where {T<:Real} = (y >= 0 ? dy * T(λ01) : dy * (y + T(λα01)))

selu(x::T) where {T<:Array} = selu.(x)
seluback(x::T, y::T, dy::T) where {T<:Array} = seluback.(x, y, dy)

@primitive1 selu(x),dy,y seluback(x,y,dy)
# @primitive1 seluback(x::Real,y::Real,dy::Real),ddx  (T=eltype(y); ddx.*((y.>=0).*T(λ01).+(y.<0).*(y.+T(λα01))))  (ddx.*dy.*(y.<0))
@zerograd seluback(x,y,dy)

const λ01 = 1.0507009873554805  # (1-erfc(1/sqrt(2))*sqrt(exp(1)))*sqrt(2pi)*(2*erfc(sqrt(2))*exp(2)+pi*erfc(1/sqrt(2))^2*exp(1)-2*(2+pi)*erfc(1/sqrt(2))*sqrt(exp(1))+pi+2)^(-0.5)
const α01 = 1.6732632423543778  # -sqrt(2/pi)/(erfc(1/sqrt(2))*exp(1/2)-1)
const λα01 = 1.7580993408473773 # λ01 * α01

"""
    sigm(x)

Return `1/(1+exp(-x))`.

Reference: Numerically stable sigm implementation from http://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick.
"""
sigm(x::T) where {T<:Real} = (x >= 0 ? T(1)/(T(1)+exp(-x)) : (z=exp(x); z/(T(1)+z)))
sigmback(x::T,y::T,dy::T) where {T<:Real} = (dy*y*(T(1)-y))

sigm(x::T) where {T<:Array} = sigm.(x)
sigmback(x::T, y::T, dy::T) where {T<:Array} = sigmback.(x, y, dy)

@primitive1 sigm(x),dy,y  sigmback(x,y,dy)
#@primitive1 sigmback(x::Real,y::Real,dy::Real),ddx  ddx.*y.*(1 .- y)  ddx.*dy.*(1 .- 2 .* y)
@zerograd sigmback(x,y,dy)

"""
    swish(x)

Return `x * sigm(x)`.

References:
* Ramachandran, P., Zoph, B., & Le, Q. V. (2017). Searching for activation functions. arXiv preprint arXiv:1710.05941.
"""
swish(x::T) where {T<:Real} = x * sigm(x)
swishback(x::T,y::T,dy::T) where {T<:Real} = dy * (x >= 0 ? (z=exp(-x);(z*x+z+1)/((z+1)^2)) : (z=exp(x);(z*(x+z+1)/((z+1)^2))))

swish(x::T) where {T<:Array} = swish.(x)
swishback(x::T, y::T, dy::T) where {T<:Array} = swishback.(x, y, dy)

@primitive1 swish(x),dy,y  swishback(x,y,dy)
# @primitive1 swishback(x::Real,y::Real,dy::Real),ddx,dx  (x2=x.^2; ddx.*dy.*(exp.(-x2 ./ 2).*(2 .- x2)./eltype(x)(sqrt(2pi))))  (ddx .* dx ./ dy)
@zerograd swishback(x,y,dy)

"""
    tanh_(x)

Return `tanh(x)`. The difference is `tanh_(x::Array)` broadcasts instead of computing a Matrix function.
"""
tanh_(x::T) where {T<:Real} = tanh(x)
tanh_back(x::T,y::T,dy::T) where {T<:Real} = dy*(1 - abs2(y))

tanh_(x::T) where {T<:Array} = tanh_.(x)
tanh_back(x::T, y::T, dy::T) where {T<:Array} = tanh_back.(x, y, dy)

@primitive1 tanh_(x),dy,y  tanh_back(x,y,dy)
# @primitive1 tanh_back(x::Real,y::Real,dy::Real),ddx,dx  (x2=x.^2; ddx.*dy.*(exp.(-x2 ./ 2).*(2 .- x2)./eltype(x)(sqrt(2pi))))  (ddx .* dx ./ dy)
@zerograd tanh_back(x,y,dy)


text:
using Dates
using Statistics
using Test

using MarketData

using TimeSeries

@testset "meta" begin

@testset "construction with and without meta field" begin
    nometa = TimeArray(timestamp(cl), values(cl), colnames(cl))

    @testset "default meta field to nothing" begin
        @test meta(nometa) == nothing
    end

    @testset "allow objects in meta field" begin
        @test meta(mdata) == "Apple"
    end
end

@testset "get index operations preserve meta" begin
    @testset "index by integer row" begin
        @test meta(mdata[1]) == "Apple"
    end

    @testset "index by integer range" begin
        @test meta(mdata[1:2]) == "Apple"
    end

    @testset "index by column name" begin
        @test meta(mdata[:Close]) == "Apple"
    end

    @testset "index by date range" begin
        @test meta(mdata[[Date(2000,1,3), Date(2000,1,14)]]) == "Apple"
    end
end

@testset "split operations preserve meta" begin
    @testset "when" begin
        @test meta(when(mdata, dayofweek, 1)) == "Apple"
    end

    @testset "from" begin
        @test meta(from(mdata, Date(2000,1,1))) == "Apple"
    end

    @testset "to" begin
        @test meta(to(mdata, Date(2000,1,1))) == "Apple"
    end
end

@testset "apply operations preserve meta" begin
    @testset "lag" begin
        @test meta(lag(mdata)) == "Apple"
    end

    @testset "lead" begin
        @test meta(lead(mdata)) == "Apple"
    end

    @testset "percentchange" begin
        @test meta(percentchange(mdata)) == "Apple"
    end

    @testset "moving" begin
        @test meta(moving(mean,mdata,10)) == "Apple"
    end

    @testset "upto" begin
        @test meta(upto(sum, mdata)) == "Apple"
    end
end

@testset "combine operations preserve meta" begin
    @testset "merge when both have identical meta" begin
        @test meta(merge(cl, op))         == "AAPL"
        @test meta(merge(cl, op, :left))  == "AAPL"
        @test meta(merge(cl, op, :right)) == "AAPL"
        @test meta(merge(cl, op, :outer)) == "AAPL"
    end

    @testset "merged meta field value concatenates when both objects' meta field values are strings" begin
        @test meta(merge(mdata, cl))         == "Apple_AAPL"
        @test meta(merge(mdata, cl, :left))  == "Apple_AAPL"
        @test meta(merge(mdata, cl, :right)) == "Apple_AAPL"
        @test meta(merge(mdata, cl, :outer)) == "Apple_AAPL"
    end

    @testset "merge when supplied with meta" begin
        @test meta(merge(mdata, mdata, meta=47))         == 47
        @test meta(merge(mdata, mdata, :left, meta=47))  == 47
        @test meta(merge(mdata, mdata, :right, meta=47)) == 47
        @test meta(merge(mdata, mdata, :outer, meta=47)) == 47
        @test meta(merge(mdata, cl, meta=47))            == 47
        @test meta(merge(mdata, cl, :left, meta=47))     == 47
        @test meta(merge(mdata, cl, :right, meta=47))    == 47
        @test meta(merge(mdata, cl, :outer, meta=47))    == 47
    end

    @testset "merged meta field value for disparate types in meta field defaults to Void" begin
        @test meta(merge(mdata, merge(cl, op, meta=47)))         == nothing
        @test meta(merge(mdata, merge(cl, op, meta=47), :left))  == nothing
        @test meta(merge(mdata, merge(cl, op, meta=47), :right)) == nothing
        @test meta(merge(mdata, merge(cl, op, meta=47), :outer)) == nothing
    end

    @testset "collapse" begin
        @test meta(collapse(mdata, week, first)) == "Apple"
    end
end

@testset "basecall operations preserve meta" begin
    @testset "basecall" begin
        @test meta(basecall(mdata, cumsum)) == "Apple"
    end
end

@testset "mathematical and comparison operations preserve meta" begin
    @testset ".+" begin
        @test meta(mdata .+ mdata) == "Apple"
        @test meta(mdata .+ cl)    == nothing
    end

    @testset ".<" begin
        @test meta(mdata .< mdata) == "Apple"
        @test meta(mdata .< cl)    == nothing
    end
end

@testset "readwrite accepts meta argument" begin
    @testset "Apple is present" begin
        @test meta(mdata) == "Apple"
    end
end

end  # @testset "meta"


text:
import .Plots: Plots

function Base.convert(t::Type{PngImage}, plot::Plots.Plot)
    pb = PipeBuffer()
    show(pb, MIME("image/png"), plot)
    return PngImage(pb)
end

preprocess(name, plot::Plots.Plot, data) = preprocess(name, convert(PngImage, plot), data)
preprocess(name, plots::AbstractArray{<:Plots.Plot}, data) = begin
    for (i, plot)=enumerate(plots)
        preprocess(name*"/$i", plot, data)
    end
    return data
end


text:
# Example for basic 2D modeling:
# The receiver positions and the source wavelets are the same for each of the four experiments.
# Author: Philipp Witte, pwitte@eos.ubc.ca
# Date: January 2017
#
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Update Date: July 2020

parsed_args = parse_commandline()

nlayer = parsed_args["nlayer"]
tti = parsed_args["tti"]
fs =  parsed_args["fs"]

### Model
model, model0, dm = setup_model(parsed_args["tti"], parsed_args["nlayer"])
q, srcGeometry, recGeometry, info = setup_geom(model)
dt = srcGeometry.dt[1]

m0 = model0.m
######################## WITH DENSITY ############################################

@testset "Jacobian test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin
    # Write shots as segy files to disk
    opt = Options(sum_padding=true, dt_comp=dt, free_surface=parsed_args["fs"])

    # Setup operators
    Pr = judiProjection(info, recGeometry)
    F = judiModeling(info, model; options=opt)
    F0 = judiModeling(info, model0; options=opt)
    Ps = judiProjection(info, srcGeometry)
    J = judiJacobian(Pr*F0*adjoint(Ps), q)

    # Linear modeling
    dobs = Pr*F0*Ps'*q
    dD = J*dm

    @test norm(J*(0f0.*dm)) == 0
    @test isapprox(dD, J*vec(dm.data); rtol=1f-6)

    # Jacobian test
    maxiter = 6
    h = 5f-2
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)

    for j=1:maxiter

        F.model.m = m0 + h*dm
        dpred = Pr*F*Ps'*q

        err1[j] = norm(dpred - dobs)
        err2[j] = norm(dpred - dobs - h*dD)
        j == 1 ? prev = 1 : prev = j - 1
        @printf("h = %2.2e, e1 = %2.2e, rate = %2.2e", h, err1[j], err1[prev]/err1[j])
        @printf(", e2 = %2.2e, rate = %2.2e \n", err2[j], err2[prev]/err2[j])

        h = h * .8f0
    end

    rate_1 = sum(err1[1:end-1]./err1[2:end])/(maxiter - 1)
    rate_2 = sum(err2[1:end-1]./err2[2:end])/(maxiter - 1)

    @test isapprox(rate_1, 1.25f0; rtol=1f-2)
    @test isapprox(rate_2, 1.5625f0; rtol=1f-2)
end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
@testset "Finite-differences" begin
    include("basis_functions.jl")
    include("scalar_operators.jl")
    include("inner_products.jl")

    @testset "Function interpolation" begin
        R = FiniteDifferences(20,1.0)
        R⁻¹ = pinv(R)
        r̃ = locs(R)
        χ = R[r̃,:]

        r = axes(R,1)

        y = Inclusion(0.0..5.0)
        @test_throws DimensionMismatch R \ y.^2

        fu = r -> r^2*exp(-r)
        u = R*(R\fu.(r))
        @test norm(χ * (R⁻¹*u) - fu.(r̃)) == 0

        fv = r -> r^6*exp(-r)
        v = R*(R\fv.(r))
        @test norm(χ * (R⁻¹*v) - fv.(r̃)) == 0
    end

    include("derivatives.jl")
end


text:
@testset "Vector" begin
    i = [1, 0, 0]
    j = [0, 1, 0]
    k = [0, 0, 1]

    @test i.x == 1
    @test i.y == 0
    @test i.z == 0
    @test_throws ErrorException i.w == 0

    @test j.x == 0
    @test j.y == 1
    @test j.z == 0
    @test_throws ErrorException j.w == 0

    @test k.x == 0
    @test k.y == 0
    @test k.z == 1
    @test_throws ErrorException k.w == 0

    v = 2i + 3j - 4k
    @test v.x == 2
    @test v.y == 3
    @test v.z == -4
    @test_throws ErrorException v.w == 0

    @test_throws AssertionError v.u == 0

    extend!(v)
    @test v.w == 0
    @test mag(v) == √(2^2 + 3^2 + (-4)^2)

    u = withangle(π, 2)
    @test u == 2 * [cos(π), sin(π)]
    @test angle(u) ≈ π

    u = withangle(π/2, abs(u))
    @test u == 2 * [cos(π/2), sin(π/2)]
    @test angle(u) == π/2

    l = [1, 1]
    @test limit(l, 1) == l / √2
end

@testset "Math" begin
    @test point([1, 2]) == Point(1, 2)

    for N in 0:10
        @test 0 <= random(N) <= N
    end

    for n in 0:10
        @test mapr(n, [0 10; 0 1]) == n / 10
        @test (n ⟶ [0 10; 0 5]) == mapr(n, [0 10; 0 5])
    end
end

@testset "SVG Path" begin
    path = "M 10 10 H 20 V 30 L 50 50 Z"

    points = getpoints(path)

    @test points == [
        Point(10, 10),
        Point(10, 20),
        Point(30, 20),
        Point(50, 50),
        Point(10, 10)
    ]
end


text:
#
# function to center the reference coordinates at the origin
#
# It modifies the vector x of coordinates
#

function tocm!( x :: Matrix{Float64}, xcm :: Vector{Float64} )
  for i in 1:length(x[:,1])
    x[i,1] = x[i,1] - xcm[1]
    x[i,2] = x[i,2] - xcm[2]
    x[i,3] = x[i,3] - xcm[3]
  end
end

max_stars_repo_name: mcubeg/JPleXL

max_issues_repo_name: mcubeg/JPleXL

max_forks_repo_name: mcubeg/JPleXL


text:
include("utils.jl")

function get_station_list(;LineCode::String = "All", IncludeAdditionalInfo::Bool = false)
    LineCode = verify_line_input(LineCode)

    if LineCode == "All"
        url = wmata.station_list_url
    else
        url = wmata.station_list_url * "?LineCode=" * LineCode
    end

    r = wmata_request(url)

    response_elements = [
        "Name",
        "Code",
        "LineCode1",
        "Lat",
        "Lon",
        "City",
        "State",
        "Street",
        "Zip"
    ]

    if IncludeAdditionalInfo == true
        append!(
            response_elements,
            ["LineCode2",
            "LineCode3",
            "LineCode4",
            "StationTogether1"]
        )
    else
        response_elements
    end

    station_list = DataFrame(
        map(response_elements) do id_col
            if id_col in ["City", "State", "Street", "Zip"]
                (id_col => [station["Address"][id_col] for station in r["Stations"]])
            else
                (id_col => [station[id_col] for station in r["Stations"]])
            end
        end
    )

    # rename some columns to be clearer
    new_names = Dict(
        "Name" => "StationName",
        "Code" => "StationCode",
        "LineCode1" => "LineCode",
        "Lat" => "Latitude",
        "Lon" => "Longitude"
    )

    for og_name in keys(new_names)
        rename!(station_list, Symbol(og_name) => Symbol(new_names[og_name]))
    end

    return station_list
end

function get_station_timings(;StationCode::String = "", StationName::String = "")
    if StationName != ""
        StationCode = get_station_code(StationName)
    else
        verify_station_input(StationCode)
    end

    url = wmata.station_timings_url * "?StationCode=" * StationCode

    r = wmata_request(url)

    days_of_week = ["Sunday" ,"Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday"]

    station_name = r["StationTimes"][1]["StationName"]
    station_code = r["StationTimes"][1]["Code"]

    opening_times = [r["StationTimes"][1][day]["OpeningTime"] for day in days_of_week]
    #=
    from messing about in the API, there are cases where first/last train info is coming back empty and
    resulting in an out of bounds error in Julia. I'm setting a drop condition so that if these
    successfully run, the function will return the full dataframe. Otherwise, first/last train info
    will be blank.
    =#
    drop_condition = "OK"
    try
        first_trains_destinations = [r["StationTimes"][:1][day]["FirstTrains"][:1]["DestinationStation"] for day in days_of_week]
    catch
        @warn "No First/Last train information available."
        drop_condition = "DROP"
    end

    if drop_condition != "DROP"
        first_trains_destinations = [r["StationTimes"][:1][day]["FirstTrains"][:1]["DestinationStation"] for day in days_of_week]
        first_trains_times = [r["StationTimes"][:1][day]["FirstTrains"][:1]["Time"] for day in days_of_week]
        last_trains_destinations = [r["StationTimes"][:1][day]["LastTrains"][:1]["DestinationStation"] for day in days_of_week]
        last_trains_times = [r["StationTimes"][:1][day]["LastTrains"][:1]["Time"] for day in days_of_week]
    else
        first_trains_destinations = ["--" for day in days_of_week]
        first_trains_times = ["--" for day in days_of_week]
        last_trains_times = ["--" for day in days_of_week]
        last_trains_destinations = ["--" for day in days_of_week]
    end

    return DataFrame(
        "StationName" => station_name,
        "StationCode" => station_code,
        "DayOfWeek" => days_of_week,
        "OpeningTime" => opening_times,
        "FirstTrainDestination" => first_trains_destinations,
        "FirstTrainTime" => first_trains_times,
        "LastTrainDestination" => last_trains_destinations,
        "LastTrainTime" => last_trains_times
    )
end

function get_rail_predictions(;StationCode::String = "All", StationName::String = "")
    if StationName != ""
        StationCode = get_station_code(StationName)
    else
        verify_station_input(StationCode)
    end

    url = wmata.rail_predictions_url * StationCode * "/"

    r = wmata_request(url)

    response_elements = [
        "LocationName",
        "LocationCode",
        "Line",
        "Car",
        "Destination",
        "Group",
        "Min"
    ]

    rail_predictions = DataFrame(
        map(
        id_col -> (id_col => [station[id_col] for station in r["Trains"]]),
        response_elements
        )
    )

    rename!(rail_predictions, :LocationName => :ArrivalStation)
    rail_predictions[!, :Min] = convert_arrival_times(rail_predictions[!, :Min])

    return rail_predictions
end

function get_path_between(;FromStationCode::String, ToStationCode::String)
    FromStationCode = verify_station_input(FromStationCode)
    ToStationCode = verify_station_input(ToStationCode)

    url = wmata.paths_url * "FromStationCode=" * FromStationCode * "&" * "ToStationCode=" * ToStationCode

    r = wmata_request(url)

    response_elements = [
        "SeqNum",
        "StationName",
        "StationCode",
        "LineCode",
        "DistanceToPrev"
    ]

    paths_between = DataFrame(
        map(
        id_col -> (id_col => [r["Path"][path_point][id_col] for path_point in 1:length(r["Path"])]),
        response_elements
        )
    )

    if nrow(paths_between) == 0
        @error "No path between stations. Did you choose stations on the same line?"
    else
        return paths_between
    end
end

function get_station_to_station(;FromStationCode::String = "", ToStationCode::String = "")
    if (FromStationCode == "" && ToStationCode == "")
        url = wmata.station_to_station_url
    else
        url = wmata.station_to_station_url * "?FromStationCode=" * FromStationCode * "&" * "ToStationCode=" * ToStationCode
    end

    r = wmata_request(url)

    response_elements = [
        "SourceStation",
        "DestinationStation",
        "CompositeMiles",
        "RailTime",
        "SeniorDisabled",
        "PeakTime",
        "OffPeakTime"
    ]

    return DataFrame(
       map(response_elements) do id_col
        if id_col in ["SeniorDisabled", "PeakTime", "OffPeakTime"]
            (id_col => [r["StationToStationInfos"][num]["RailFare"][id_col] for num in 1:length(r["StationToStationInfos"])])
        else
            (id_col => [r["StationToStationInfos"][num][id_col] for num in 1:length(r["StationToStationInfos"])])
        end
        end
    )
end

function get_train_positions()
    r = wmata_request(wmata.train_positions_url)
    train_positions = r["TrainPositions"]

    response_elements = [
    "CarCount",
    "DestinationStationCode",
    "DirectionNum",
    "LineCode",
    "SecondsAtLocation",
    "ServiceType",
    "TrainId",
    "TrainNumber",
    "CircuitId"
    ]

    return DataFrame(
        map(
        id_col -> (id_col => [train[id_col] for train in train_positions]),
        response_elements
        )
    )
end

function get_rail_incidents()
    r = wmata_request(wmata.rail_incidents_url)

    response_elements = [
    "IncidentID",
    "Description",
    "EndLocationFullName",
    "PassengerDelay",
    "LinesAffected",
    "IncidentType",
    "DateUpdated"
    ]

    DataFrame(
        map(response_elements) do id_col
        if id_col == "LinesAffected"
            lines_affected = [r["Incidents"][incident][id_col] for incident in 1:length(r["Incidents"])]
            ("LinesAffected" => map(x -> split(replace(x, " " => ""), ';', keepempty = false), lines_affected))
        else
            (id_col => [r["Incidents"][incident][id_col] for incident in 1:length(r["Incidents"])])
        end
        end
    )
end


text:
include("../functions/Functions.jl")
include("../evaluate/Evaluate.jl")
include("../opt/SAGA.jl")
include("../data/data.jl")

X, y = getOldData()

x0 = zeros(size(X,2),1)

H = logreg_hessian(X,y,x0)

(vec,val) = eig(H)
println(vec)
println(val)

max_stars_repo_name: zagazao/num-opt

max_issues_repo_name: zagazao/num-opt

max_forks_repo_name: zagazao/num-opt


text:
#!/usr/bin/env julia

function main()
    values = zeros(UInt64, 256)
    values[UInt8('A')+1] = values[UInt8('T')+1] = 1
    values[UInt8('G')+1] = values[UInt8('C')+1] = 1 << 32

    input = open("chry_multiplied.fa")
    buflen = 4096
    buf = zeros(UInt8, buflen)

    n = readbytes!(input, buf, buflen)

    count = UInt64(0)
    mask = UInt64(0)
    while n > 0
        if n < buflen
            resize!(buf, n)
        end

        for c in buf
            if c == UInt8('\n')
                mask = ~UInt64(0)
            elseif c == UInt8('>')
                mask = UInt64(0)
            else
                count += mask & @inbounds values[c+1]
            end
        end

        n = readbytes!(input, buf, buflen)
    end

    at_count = count & 0xffffffff
    gc_count = count >> 32

    gc = gc_count / (gc_count + at_count)
    println(gc)
end

main()

max_stars_repo_name: pditommaso/gccontent-benchmark

max_issues_repo_name: pditommaso/gccontent-benchmark

max_forks_repo_name: pditommaso/gccontent-benchmark


text:
function derivative(f::Function, ftype::Symbol, dtype::Symbol)
  if ftype == :scalar
    g(x::Number) = finite_difference(f, x, dtype)
  elseif ftype == :vector
    g(x::Vector) = finite_difference(f, x, dtype)
  else
    error("ftype must :scalar or :vector")
  end
  return g
end
derivative{T <: Number}(f::Function, x::Union(T, Vector{T}), dtype::Symbol) = finite_difference(f, x, dtype)
derivative{T <: Number}(f::Function, x::Union(T, Vector{T})) = finite_difference(f, x, :central)
derivative(f::Function, dtype::Symbol) = derivative(f, :scalar, dtype)
derivative(f::Function) = derivative(f, :scalar, :central)

gradient(f::Function, dtype::Symbol) = derivative(f, :vector, dtype)
gradient{T <: Number}(f::Function, x::Union(T, Vector{T}), dtype::Symbol) = finite_difference(f, x, dtype)
gradient(f::Function) = derivative(f, :vector, :central)
gradient{T <: Number}(f::Function, x::Union(T, Vector{T})) = finite_difference(f, x, :central)

ctranspose(f::Function) = derivative(f)

function jacobian{T <: Number}(f::Function, x::Vector{T}, dtype::Symbol)
    finite_difference_jacobian(f, x, dtype)
end
function jacobian(f::Function, dtype::Symbol)
    g(x::Vector) = finite_difference_jacobian(f, x, dtype)
    return g
end
jacobian(f::Function) = jacobian(f, :central)

function second_derivative(f::Function, g::Function, ftype::Symbol, dtype::Symbol)
  if ftype == :scalar
    h(x::Number) = finite_difference_hessian(f, g, x, dtype)
  elseif ftype == :vector
    h(x::Vector) = finite_difference_hessian(f, g, x, dtype)
  else
    error("ftype must :scalar or :vector")
  end
  return h
end
function second_derivative{T <: Number}(f::Function, g::Function, x::Union(T, Vector{T}), dtype::Symbol)
  finite_difference_hessian(f, g, x, dtype)
end
function hessian{T <: Number}(f::Function, g::Function, x::Union(T, Vector{T}), dtype::Symbol)
  finite_difference_hessian(f, g, x, dtype)
end
function second_derivative{T <: Number}(f::Function, g::Function, x::Union(T, Vector{T}))
  finite_difference_hessian(f, g, x, :central)
end
function hessian{T <: Number}(f::Function, g::Function, x::Union(T, Vector{T}))
  finite_difference_hessian(f, g, x, :central)
end
function second_derivative(f::Function, x::Number, dtype::Symbol)
  finite_difference_hessian(f, derivative(f), x, dtype)
end
function hessian(f::Function, x::Number, dtype::Symbol)
  finite_difference_hessian(f, derivative(f), x, dtype)
end
function second_derivative{T <: Number}(f::Function, x::Vector{T}, dtype::Symbol)
  finite_difference_hessian(f, gradient(f), x, dtype)
end
function hessian{T <: Number}(f::Function, x::Vector{T}, dtype::Symbol)
  finite_difference_hessian(f, gradient(f), x, dtype)
end
function second_derivative(f::Function, x::Number)
  finite_difference_hessian(f, derivative(f), x, :central)
end
function hessian(f::Function, x::Number)
  finite_difference_hessian(f, derivative(f), x, :central)
end
function second_derivative{T <: Number}(f::Function, x::Vector{T})
  finite_difference_hessian(f, gradient(f), x, :central)
end
function hessian{T <: Number}(f::Function, x::Vector{T})
  finite_difference_hessian(f, gradient(f), x, :central)
end
second_derivative(f::Function, g::Function, dtype::Symbol) = second_derivative(f, g, :scalar, dtype)
second_derivative(f::Function, g::Function) = second_derivative(f, g, :scalar, :central)
second_derivative(f::Function) = second_derivative(f, derivative(f), :scalar, :central)
hessian(f::Function, g::Function, dtype::Symbol) = second_derivative(f, g, :vector, dtype)
hessian(f::Function, g::Function) = second_derivative(f, g, :vector, :central)
hessian(f::Function) = second_derivative(f, gradient(f), :vector, :central)


text:
"""
Cost fileformat
See EUROCONTROL NEST Manual for Cost fileformat description

Example cost fileformat
139729486 ED 356.159
139729486 LB 230.25
139729486 LC 86.7351
139729486 LH 148.107
139729486 LO 234.767
139729486 LT 365.319
139729486 LY 190.189
"""
module Cost

export read

using CSV

const fileformat = Dict(1=>Int64, 2=>String, 3=>Float64)

const header_format = ["FLIGHTID", "COUNTRYCODE", "COST"]

function read(file)
    return CSV.read(file, delim=" ", types=fileformat, header=header_format,
    copycols=true)
end

end # module


text:
"""
    Indices(iter)

Construct an `Indices <: AbstractIndices` from an arbitrary Julia iterable. Lookup of
elements uses naive iteration.

If `iter` is a `Vector` then the result is `isinsertable` and `istokenizable`, making simple
and flexible indices using naive search that may be optimal for small collections.
"""
struct Indices{I, Inds} <: AbstractIndices{I}
	inds::Inds
end

Indices() = Indices([])
Indices{I}() where {I} = Indices{I}(I[])
Indices{I, Inds}() where {I, Inds} = Indices{I}(Inds())

Indices(iter) = Indices{eltype(iter), typeof(iter)}(iter) # EltypeUnknown...
Indices{I}(iter) where {I} = Indices{I, typeof(iter)}(iter) # There is a corner case where the elements of `iter` might not be of type `I`

# Basic interface
@propagate_inbounds function Base.iterate(i::Indices{I}, state...) where {I}
    # We haven't enforced uniqueness anywhere...
    iterate(i.inds, state...)::Union{Nothing, Tuple{I, Any}}
end

Base.in(i::I, inds::Indices{I}) where {I} = in(i, inds.inds)
Base.IteratorSize(i::Indices) = Base.IteratorSize(i.inds)
Base.length(i::Indices) = length(i.inds)

# Specialize for `Vector` elements. Satisfy the tokenization and insertion interface
const VectorIndices{I} = Indices{I, Vector{I}}

istokenizable(i::VectorIndices) = true
tokentype(::VectorIndices) = Int
@inline iteratetoken(inds::VectorIndices, s...) = iterate(keys(inds.inds), s...)
@inline function gettoken(inds::VectorIndices{I}, i::I) where {I}
    @inbounds for x in keys(inds.inds)
        if isequal(i, inds.inds[x])
            return (true, x)
        end
    end
    return (false, 0)
end
@propagate_inbounds gettokenvalue(inds::VectorIndices, x::Int) = inds.inds[x]

isinsertable(i::VectorIndices) = true
@inline function gettoken!(inds::VectorIndices{I}, i::I) where {I}
    @inbounds for x in keys(inds.inds)
        if isequal(i, inds.inds[x])
            return (true, x)
        end
    end
    push!(inds.inds, i)
    return (false, length(inds.inds))
end

@inline function deletetoken!(inds::VectorIndices, x::Int)
    deleteat!(inds.inds, x)
    return inds
end

function Base.empty!(inds::VectorIndices)
    empty!(inds.inds)
    return inds
end

Base.empty(inds::VectorIndices, ::Type{I}) where {I} = Indices{I, Vector{I}}(Vector{I}())

Base.copy(inds::VectorIndices) = Indices(copy(inds.inds))


text:
using Publish
using FastAI
using Pkg.Artifacts

# Download theme
download_artifact(
    Base.SHA1("6e4be8bec8da9323c18d777d7855ef79dddcf524"),
    "https://github.com/darsnack/flux-theme/releases/download/v0.2.1/flux-theme-0.2.1.tar.gz",
    "b00941248ecdb643a72960bbfda19a3128591eca3d7cbcb3bfb80a6ab1c7f99e",
)

Publish.Themes.default() = artifact"flux-theme"
ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

p = Publish.Project(FastAI)

rm("dev", recursive = true, force = true)
rm(p.env["version"], recursive = true, force = true)

deploy(FastAI; root = "/FastAI.jl", force = true, label = "dev")


text:
using Test

using MPI

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

# Not possible to PROD a Char (and neither Int8 with OpenMPI)
typs = setdiff([Base.uniontypes(MPI.MPIDatatype)...], [Char, Int8, UInt8])
for typ in typs
    val = convert(typ,rank+1)
    B = MPI.Exscan(val, MPI.PROD, comm)
    if rank > 0
        @test B[1] === convert(typ, prod(1:rank))
    end
end

MPI.Finalize()


text:
#concrete type of trees
abstract Decision
typealias DecisionNode Node{Decision}
typealias DecisionLeaf Leaf{Decision}

type DecisionTree <: Tree{Decision}
   root::DecisionNode
end

type DecisionBranch <: Branch{Decision}
    feature::Int
    value::Float64
    left::DecisionNode
    right::DecisionNode
end

type RegressionLeaf <: DecisionLeaf
  value::Float64
end

type ClassificationLeaf <: DecisionLeaf
  prob::Vector{Float64}
end

max_stars_repo_name: mathcg/BayesTree

max_issues_repo_name: mathcg/BayesTree

max_forks_repo_name: mathcg/BayesTree


text:
#Init the Net
function init_θ_pol_SAC(net_s::Array{Int,1})
    """
    This network will contain weights for actions and log_std
    """
    x      = net_s[1]
    hidden = net_s[2:end-1]
    ysize  = net_s[end]

    w = []
    for y in [hidden...]
        push!(w, xavier(y, x))
        push!(w, xavier(y, 1))
        x = y
    end
    push!(w, xavier(ysize, x)) #Weights outputting the mean
    push!(w, zeros(ysize, 1))

    push!(w, xavier(ysize, x)) #Weights outputting the log_std
    push!(w, zeros(ysize, 1))

    return w
end

#Init the Net
function init_θ_pol_SAC_norm(net_s::Array{Int,1})
    """
    This network will contain weights for actions and log_std
    """
    x      = net_s[1]
    hidden = net_s[2:end-1]
    ysize  = net_s[end]

    w = []
    for y in [hidden...]
        push!(w, xavier(y, x))
        push!(w, ones(y))
        push!(w, zeros(y))
        x = y
    end
    push!(w, xavier(ysize, x)) #Weights outputting the mean
    push!(w, zeros(ysize, 1))

    push!(w, xavier(ysize, x)) #Weights outputting the log_std
    push!(w, zeros(ysize, 1))

    return w
end

#Evaluate the Policy-Network
@everywhere function mlp_nlayer_policy_SAC(w::Any,x::AbstractVector)
    """
    Evaluation of the network
    """

    LOG_SIG_MAX = 2.0
    LOG_SIG_MIN = -20.0

    #Go over the hidden layers:
    for i = 1:2:length(w)-4
        x = relu.(w[i] * x .+ w[i+1])
    end

    mu_act   = vec(w[end-3] * x .+ w[end-2]) # without tanh here
    log_std  = vec(w[end-1] * x .+ w[end])

    std = exp.(log_std)

    return mu_act, std
end

#Evaluate the Policy-Network
@everywhere function mlp_nlayer_policy_SAC_norm(w::Any,x::AbstractVector)
    """
    Evaluation of the network
    """

    LOG_SIG_MAX = 2.0
    LOG_SIG_MIN = -20.0

    #Go over the hidden layers:
    for i = 1:3:length(w)-4
        #1. Get the new summed input
        x = w[i] * x

        #2. Calculate the mean and std of the layer
        x_mean, x_std = get_norm(x)

        #3. Normalize activation by applying the norm and mean
        x_norm = (x .- x_mean) / x_std

        #4. Reparametrize the output
        x = x_norm .* w[i+1] + w[i+2]

        #5. Apply the non-linearity
        x = relu.(x)
    end

    mu_act   = vec(w[end-3] * x .+ w[end-2]) # without tanh here
    log_std  = vec(w[end-1] * x .+ w[end])

    std = exp.(log_std)

    return mu_act, std
end

#Get the required action
@everywhere function sample_action_SAC(μ::AbstractVector, σ::AbstractVector, rand_vec::AbstractVector, params::Params)
    u_range = params.u_range
    μ       = convert(Array{Float64}, μ)
    sample  = μ .+ rand_vec .* σ
    actions = u_range.*tanh.(sample)
    #Return both values for probability calculation:
    return actions, sample
end

#Get the log(pi(a|s))
@everywhere function logpdf_SAC(a_t,out_t,std)
    fac = Float64(-log(sqrt(2pi)))
    r = (a_t-out_t) ./ std
    return -r.* r.* Float64(0.5) - (log.(std)) .+ fac
end

@everywhere function run_episode_SAC!(params::Params,Pl::Any, Vf::Any, traj::Trajectory, rms::RunMeanStd)
    dt          = params.dt
    t_horizon   = params.t_horizon
    s_init      = params.s_init + 0.1*sum(params.σ)*randn(length(params.s_init))
    l_traj      = params.l_traj
    σ           = params.σ
    u_range     = params.u_range
    l_out       = length(u_range)

    #Initialize :
    s_t1 = copy(s_init)
    s_t2 = similar(s_t1)
    totalr = 0

    for i_t = 1:l_traj
        #Normalize the Input
        #s_norm = normalize_rms(rms,s_t1) #select for normalization
        s_norm = s_t1 #select for unnormalized output

        #Evaluate the Network Output
        out_t, std = mlp_nlayer_policy_SAC_norm(Pl,s_norm)
        val_t      = mlp_nlayer_value_norm(Vf,s_norm,params)

        #Get Action with stochastic policy:
        rand_vec         = randn(l_out)
        loc_a_t, sample  = sample_action_SAC(out_t, std, rand_vec, params)
        p_a_t            = get_π_prob(out_t,std,sample)

        #Interact with the Environment:
        s_t2, loc_r_t, terminal = environment_call(s_t1,loc_a_t,params)

        #Save data in trajectory:
        traj.s_t[i_t,:]    .= s_norm
        traj.s_t2[i_t,:]   .= s_t2
        traj.a_t[i_t,:]    .= loc_a_t
        traj.p_a_t[i_t,:]   = p_a_t
        traj.r_t[i_t]       = 0.01 * loc_r_t  #* -1 #REWARD NEGATED HERE!!!!
        traj.out_t[i_t,:]   = out_t
        traj.val_t[i_t]     = val_t
        traj.len            = i_t
        totalr             += loc_r_t

        #Override the state:
        s_t1 .= s_t2;

        #Check if the terminal condition is fullfilled and break
        if terminal == 1
                break
        end

    end
    return totalr / traj.len
end

"""
Function enabling parallel execution of the episodes by setting up the channels
"""
function parallel_episodes_SAC(params::Params,Pl::Any,Vf::Any,all_traj::Traj_Data,rms_state::RunMeanStd)
        n_traj = params.n_traj
        active_workers = workers()

        # Create Channels storing results and starting valus
        s_init_channel = RemoteChannel(()->Channel{Vector{Float64}}(length(workers())+1)) # will contain s_init
        result_channel = RemoteChannel(()->Channel{Tuple}(length(workers())+1)) # will contain traj and gradients

        # Create problem on every worker
        begin
            for p in workers()
                remote_do(distributed_episodes_SAC, p, params, Pl, Vf, rms_state, s_init_channel, result_channel)
            end
        end

        # Start optimizations
        for i_traj = 1:n_traj
                s_init_local = params.s_init
                @async put!(s_init_channel, s_init_local)
        end

        # Collect the resulting data
        prec = 0
        res = 0
        prog = Progress(n_traj,1)
        for i_traj = 1:n_traj
                #Take from Channel:
                result = take!(result_channel)
                #Oder the Tuple:
                all_traj.traj[i_traj] = result[1]
                res  += result[2]
                prec += result[3]
        end

        # Tell workers that its done
        for j = 1:length(workers())
                NaN_mat = zeros(length(params.s_init))
                fill!(NaN_mat,NaN)
                @async put!(s_init_channel, NaN_mat)
        end

        #Close Channels
        sleep(0.1)
        close(s_init_channel)
        close(result_channel)

        #Return the Feedback
        return (res, prec)
end

# Function running on the worker
@everywhere function distributed_episodes_SAC(params::Params, Pl::Any, Vf::Any, rms_state::RunMeanStd, s_init_channel::RemoteChannel, result_channel::RemoteChannel)

    # Preallocate trajectory
    traj = init_traj(params,params.l_traj)
    s_init_local = zeros(length(params.s_init))
    # Start loop where work is done
    while true
        # break loop if NaN state received
        try
            s_init_local .= take!(s_init_channel)
            if (sum(isnan, s_init_local) > 0)
                break
            end

        catch
            println("Some problem occured, worker $(myid()) exiting trajectory optimization function.")
            break
        end

        #DO THE WORK:
        #-----------------------------------------------------------------------
        res = run_episode_SAC!(params,Pl, Vf, traj, rms_state)

        #Collect Precision
        prec = loss_value(Vf,params,traj)
        #-----------------------------------------------------------------------

        put!(result_channel, (traj, res, prec))
    end

    return nothing
end

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
using Test
using LinearAlgebra
using GeoParams

@testset "MeltingParam.jl" begin

    #Make sure structure is isbits
    x = MeltingParam_Caricchi()
    @test isbits(x)

    # This tests the various melting parameterizations
    CharUnits_GEO = GEO_units(; viscosity=1e19, length=10km)

    T = collect(250:100:1250) * K .+ 273.15K
    T_nd = Float64.(T / CharUnits_GEO.Temperature)

    # Caricchi parameterization [in ND numbers, which is anyways the typical use case]
    p = MeltingParam_Caricchi()
    phi_dim = zeros(size(T))
    args=(;T=ustrip.(T))
    compute_meltfraction!(phi_dim, p, args)

    phi_dim1 = zeros(size(phi_dim))
    compute_meltfraction!(phi_dim1, p, args) # in-place routine

    p_nd = p
    p_nd = nondimensionalize(p_nd, CharUnits_GEO)
    phi_nd = zeros(size(T))
    args=(;T=T_nd)
    compute_meltfraction!(phi_nd, p_nd, args)

    # Do this computation manually, using the actual expression of Caricchi
    T_C = collect(250:100:1250)    # in celcius
    Phi_solid = 1.0 .- 1.0 ./ (1.0 .+ exp.((800.0 .- T_C) ./ 23.0))
    Phi_anal = 1.0 .- Phi_solid

    @test sum(phi_dim - Phi_anal) < 1e-12
    @test sum(phi_dim1 - Phi_anal) < 1e-12
    @test sum(phi_nd - Phi_anal) < 1e-12

    # test derivative vs T
    dϕdT_dim = zeros(size(T))
    args=(;T=ustrip.(T))
    compute_dϕdT!(dϕdT_dim, p, args)
    @test sum(dϕdT_dim) ≈ 0.008102237679214096

    #------------------------------
    # 5th order polynomial
    p = MeltingParam_5thOrder()
    compute_meltfraction!(phi_dim, p, args)
    @test sum(phi_dim) ≈ 4.708427909521561

    # experimental data to create the fit
    data = [
        1115 100
        1050 90
        969 85
        932 80
        907 70
        880 60
        850 54
        825 51.7
        800 52.9
        775 46.3
        750 44.9
        725 29.9
        700 14.9
        690 0
    ]

    data[:, 2] = data[:, 2] / 100
    Tdata = data[:, 1] .+ 273.15
    phi = zeros(size(Tdata))
    args = (;T=ustrip.(Tdata))
    compute_meltfraction!(phi, p, args)

    @test norm(data[:, 2] - phi) ≈ 0.07151515017819135

    # test derivative vs T
    dϕdT_dim = zeros(size(T))
    args = (;T=ustrip.(T))
    compute_dϕdT!(dϕdT_dim, p, args)
    @test sum(dϕdT_dim) ≈ 0.006484458453421382

    #------------------------------

    #------------------------------
    # 4th order polynomial
    p = MeltingParam_4thOrder()
    compute_meltfraction!(phi_dim, p, args)
    @test sum(phi_dim) ≈ 4.853749635538406

    # experimental data to create the fit
    data = [
        1000 100
        990 100
        975 93
        950 89.2
        925 76.3
        900 69.6
        875 59
        850 54
        825 51.7
        800 52.9
        775 46.3
        750 44.9
        725 29.9
        700 14.9
        690 0
    ]

    data[:, 2] = data[:, 2] / 100
    Tdata = data[:, 1] .+ 273.15
    phi = zeros(size(Tdata))
    args = (;T=ustrip.(Tdata))
    compute_meltfraction!(phi, p, args)

    @test norm(data[:, 2] - phi) ≈ 0.0678052542705406

    # test derivative vs T
    dϕdT_dim = zeros(size(T))
    args = (;T=ustrip.(T))
    compute_dϕdT!(dϕdT_dim, p, args)
    @test sum(dϕdT_dim) ≈ 0.00830985782591842

    #------------------------------

    #------------------------------
    # Quadratic parameterisation
    p = MeltingParam_Quadratic()
    compute_meltfraction!(phi_dim, p, args)
    @test sum(phi_dim) ≈ 5.0894901144641

    dϕdT_dim = zeros(size(T))
    compute_dϕdT!(dϕdT_dim, p, args)
    @test sum(dϕdT_dim) ≈ 0.009365244536940681
    #------------------------------

    #------------------------------
    # Assimilation parameterisation
    p        =  MeltingParam_Assimilation();
    compute_meltfraction!(phi_dim, p, args)
    @test sum(phi_dim) ≈ 4.995
    dϕdT_dim =  zeros(size(T))
    compute_dϕdT!(dϕdT_dim, p, args)
    @test sum(abs.(dϕdT_dim)) ≈0.004605170185988078
    #------------------------------

    # Test computation of melt parameterization for the whole computational domain, using arrays
    MatParam = Vector{MaterialParams}(undef, 4)
    MatParam[1] = SetMaterialParams(;
        Name="Mantle", Phase=1, Melting=PerpleX_LaMEM_Diagram("test_data/Peridotite.in")
    )

    MatParam[2] = SetMaterialParams(;
        Name="Crust", Phase=2, Melting=MeltingParam_Caricchi()
    )

    # No melting parameterization for this phase
    MatParam[3] = SetMaterialParams(;
        Name="UpperCrust", Phase=3, Melting=MeltingParam_5thOrder(), Density=PT_Density()
    )

    # No melting parameterization for this phase
    MatParam[4] = SetMaterialParams(; Name="LowerCrust", Phase=4, Density=PT_Density())

    Mat_tup = Tuple(MatParam)

    # test computing material properties
    n = 100
    Phases = ones(Int64, n, n, n)
    Phases[:, :, 20:end] .= 2
    Phases[:, :, 80:end] .= 3
    Phases[:, :, 90:end] .= 4

    ϕ = zeros(size(Phases))
    dϕdT = zeros(size(Phases))
    T = ones(size(Phases)) * 1500
    P = ones(size(Phases)) * 10
    args = (P=P, T=T)
    compute_meltfraction!(ϕ, Mat_tup, Phases, args) #allocations coming from computing meltfraction using PhaseDiagram_LookupTable
    @test sum(ϕ) / n^3 ≈ 0.7463001302812086

    compute_dϕdT!(dϕdT, Mat_tup, Phases, args)
    @test sum(dϕdT) / n^3 ≈ 0.000176112129245805

    # test computing material properties when we have PhaseRatios, instead of Phase numbers
    PhaseRatio = zeros(n, n, n, 4)
    for i in CartesianIndices(Phases)
        iz = Phases[i]
        I = CartesianIndex(i, iz)
        PhaseRatio[I] = 1.0
    end

    compute_meltfraction!(ϕ, Mat_tup, PhaseRatio, args)
    @test sum(ϕ) / n^3 ≈ 0.7463001302812086

    compute_dϕdT!(dϕdT, Mat_tup, PhaseRatio, args)
    @test sum(dϕdT) / n^3 ≈ 0.000176112129245805

    # Test smoothening of the melting curves:
    p = SmoothMelting(p=MeltingParam_5thOrder())
    @test isbits(p)
    T = collect(250:100:1250) * K .+ 273.15K
    phi_dim = zeros(size(T))
    args=(;T=ustrip.(T))
    compute_meltfraction!(phi_dim, p, args)
    @test sum(phi_dim) ≈ 4.7084279086574226

    dϕdT = zeros(size(T))
    compute_dϕdT!(dϕdT, p, args)
    @test sum(dϕdT) ≈ 0.006484456307540648

    # try non-dimensionalisation
    p_nd = nondimensionalize(p, CharUnits_GEO)
    @test isbits(p)
    @test isdimensional(p_nd.p.a)==false
    @test p_nd.p.a ≈  6968.639721576996

    p1= dimensionalize(p_nd,CharUnits_GEO)
    @test isdimensional(p1.p.a)==true
    @test Value(p1.p.a) ≈ Value(p.p.a)

    #Mat_tup = ( SetMaterialParams( Name="Mantle", Phase=1, Melting=SmoothMelting(MeltingParam_4thOrder())),
    #            SetMaterialParams( Name="Crust", Phase=2, Melting=MeltingParam_5thOrder()),
    #            SetMaterialParams( Name="UpperCrust", Phase=3, Melting=SmoothMelting(MeltingParam_5thOrder()), Density=PT_Density()),
    #            SetMaterialParams( Name="LowerCrust", Phase=4, Density=PT_Density())
    #                )
    Mat_tup = ( SetMaterialParams( Name="Mantle",       Phase=1, Melting=SmoothMelting(MeltingParam_4thOrder())),
                SetMaterialParams( Name="Crust",        Phase=2, Melting=MeltingParam_5thOrder()),
                SetMaterialParams( Name="UpperCrust",   Phase=3, Melting=SmoothMelting(MeltingParam_5thOrder()), Density=PT_Density()),
                SetMaterialParams( Name="LowerCrust",   Phase=4, Density=PT_Density())
                )

    ϕ = zeros(size(Phases));
    dϕdT = zeros(size(Phases));
    T = ones(size(Phases)) * (800+273.15);
    P = ones(size(Phases)) * 10;
    args = (P=P, T=T);
    compute_meltfraction!(ϕ, Mat_tup, Phases, args); #allocation free
    @test sum(ϕ) / n^3 ≈ 0.4409766063389861

    compute_dϕdT!(dϕdT, Mat_tup, Phases, args) #allocation free
    @test sum(dϕdT) / n^3 ≈ 0.0006838372430250584

end


text:
#
# Computes the ordered internal distance matrix of the coordinates of
# a structure
#
# Input: x[n,3] the 3D coordinates of the atoms
#
# Returns: D, the internal distance matrix of the x
#          Dind, the indices of the elements of D, ordered from lower
#          the higher at each column
#

function intDmat(x :: Matrix{Float64})

  n = length(x[:,1])
  D = Matrix{Float64}(undef,n,n)
  Dind = Matrix{Int64}(undef,n,n)

  # Compute disprot, the internal square-distance matrix

  for i in 1:n
    D[i,i] = 0.
  end

  for j in 1:n-1
    for i in j+1:n
      dij = ( x[i,1] - x[j,1] )^2 + ( x[i,2] - x[j,2] )^2 + ( x[i,3] - x[j,3] )^2
      D[i,j] = dij
      D[j,i] = dij
    end
  end

  # Sort the distances to atom i of the protein

  for j in 1:n
    ind = [ i for i in 1:n ]
    Dind[:,j] = sort(ind, by = i -> D[i,j] )
  end

  return D, Dind

end


text:
mutable struct FTestResult <: TestResult
    levels::Union{Vector{String}, Nothing}
    n_obs::IntVector
    sds::FloatVector
    statistic::Float64
    dofs::IntVector
    pval::Float64
end

function show(io::IO, ftr::FTestResult)
    println("\n      F test to compare two variances")
    println("Data (± SD)")
    @printf(io, "  %s: ± %.2f\n",
            ftr.levels[1], ftr.sds[1])
    @printf(io, "  %s: ± %.2f\n",
            ftr.levels[2], ftr.sds[2])
    println(io, "\nResult:\n", coeftable(ftr))
end

function coeftable(ftr::FTestResult)
    CoefTable(
        [ftr.statistic, ftr.dofs[1], ftr.dofs[2], ftr.pval],
        ["F", "DF₁", "DF₂", "P-value"],
        [""],
        4, 1
    )
end

"""
    f_test(x1, x2)

F-test

The arguments `x1` and `x2` can be numeric `Vector`.
"""
function f_test(
    x1::NumVector,
    x2::NumVector;
    levels::Union{Vector{String}, Nothing}=nothing
)::FTestResult
    # use temporal sample names
    if isnothing(levels)
        @info "Using the first and second sample names as 'Sample1' and 'Sample2'"
        levels = ["Sample1", "Sample2"]
    end

    # check types
    if !(x1 isa FloatVector)
        @warn "$levels[1] is converted to Float64"
        x1 = convert(Vector{Float64}, x1)
    end
    if !(x2 isa FloatVector)
        @warn "$levels[2] is converted to Float64"
        x2 = convert(Vector{Float64}, x2)
    end

    # check length
    n1 = length(x1)
    n2 = length(x2)
    n1 > 1 || error("size of $levels[1] is not enough: actual size is $n1")
    n2 > 1 || error("size of $levels[2] is not enough: actual size is $n2")

    # perform f-test
    v1 = var(x1)
    v2 = var(x2)
    ν1 = v1 < v2 ? n2-1 : n1-1
    ν2 = v1 < v2 ? n1-1 : n2-1
    f_stat = v1 < v2 ? v2/v1 : v1/v2
    pval = 2.0 * D.ccdf(FDist(ν1, ν2), f_stat)
    res = FTestResult(levels, [n1, n2], [sqrt(v1), sqrt(v2)],
                      f_stat, [ν1, ν2], pval)
    return res
end

"""
    FTestModel

A `TestModel` type for f-test

# Members

- `X`: model matrix
- `y`: vector to compare its means
- `levels`: `String` vector of the unique levels (only 2 elements are allowed)

"""
struct FTestModel{T<:Real} <: TestModel
    x1::AbstractVector{T}
    x2::AbstractVector{T}
    levels::Vector{String}
    function FTestModel(
        X::Matrix{T},
        y::AbstractVector{T},
        levels::Vector{String};
    ) where T
        dummy = X[:,2]
        x1 = y[dummy .== 0]
        x2 = y[dummy .== 1]
        new{T}(x1, x2, levels)
    end
end

function fit(obj::FTestModel)
    f_test(obj.x1, obj.x2, levels=obj.levels)
end
function fit(
    ::Type{FTestModel},
    X::AbstractMatrix{<:Real},
    y::AbstractVector{<:Real},
    levels::Vector{String}
)
    ftest_obj = FTestModel(X, y, levels)
    fit(ftest_obj)
end

"""
    f_test(X, y)

An alias for `fit(FTestModel, X, y, levels)`

The arguments `X` and `y` can be a `Matrix` and a `Vector` or a `Formula` and `DataFrame`.

The argument `levels` can be a `Vector{String}`.
"""
function f_test(
    X::AbstractMatrix{<:Real},
    y::AbstractVector{<:Real}
)
    @info "Using the first and second sample names as 'Sample1' and 'Sample2'"
    fit(FTestModel, X, y, ["Sample1", "Sample2"])
end
function f_test(
    formula::FormulaTerm{Term,Term},
    df::DataFrame
)
    # obtain levels
    f = apply_schema(formula, schema(formula, df))
    levels = DataAPI.levels(df[!,f.rhs.terms[1].sym])
    if eltype(levels) <: Real
        levels = string.(levels)
    end

    # model matrix
    model_frame = ModelFrame(formula, df)
    X = ModelMatrix(model_frame).m
    y = response(model_frame)

    # perform f-test
    fit(FTestModel, X, y, levels)
end