text:
# In this tutorial, we will learn
#
#    -  How to solve a surface-coupled multi-physics problem.
#    -  Construct FE spaces defined in different domains.
#    -  Define interface triangulations and integrate on them.
#
# 1. [Problem Statement](#probStat)
#    - [Strong form](#strongForm)
#    - [Geometry and Discrete model](#geometry)
#    - [Boundary conditions and properties](#conditions)
# 2. [Numerical scheme](#numericalScheme)
#    - [FE spaces](#feSpace)
#    - [Weak form](#weakForm)
#    - [Numerical integration](#integration)
#    - [Algebraic system of equations](#algebraic)
# 3. [Post-processing](#postprocess)
#    - [Visualization](#visualization)
#    - [Quantities of interest](#QOIs)

# ```@raw HTML
# <a name="probStat"></a>
# ```
# ## Problem statement
# ```@raw HTML
# <a name="strongForm"></a>
# ```
# ### Strong form
#  Let $\Gamma_{\rm FS}$ be the interface between a fluid domain $\Omega_{\rm F}$ and a solid domain $\Omega_{\rm S}$. We denote by $\Gamma_{\rm F,D}$ and $\Gamma_{\rm F,N}$ the fluid boundaries with Dirichlet and Neumann conditions, respectively.
# The Fluid-Structure Interaction (FSI) problem reads:
#
# find $u_{\rm F}$, $p_{\rm F}$ and $u_{\rm S}$ such that
# ```math
# \left\lbrace
# \begin{aligned}
# -\nabla\cdot\boldsymbol{\sigma}_{\rm F} = f &\text{ in }\Omega_{\rm F},\\
# \nabla\cdot u_{\rm F} = 0 &\text{ in } \Omega_{\rm F},\\
# -\nabla\cdot\boldsymbol{\sigma}_{\rm S} = s &\text{ in }\Omega_{\rm S},\\
# \end{aligned}
# \right.
# ```
#
# satisfying the Dirichlet and Neumann boundary conditions
# ```math
# \left\lbrace
# \begin{aligned}
# u_{\rm F} = g_{\rm F} &\text{ on } \Gamma_{\rm F,D},\\
# n_{\rm F} \cdot \boldsymbol{\sigma}_{\rm F} = 0 &\text{ on } \Gamma_{\rm F,N},\\
# u_{\rm S} = g_{\rm S} &\text{ on } \Gamma_{\rm S,D},\\
# \end{aligned}
# \right.
# ```
#
# and the kinematic and dynamic conditions at the fluid-solid interface
# ```math
# \left\lbrace
# \begin{aligned}
# u_{\rm F} = u_{\rm S} &\text{ on } \Gamma_{\rm FS},\\
# n_{\rm F} \cdot \boldsymbol{\sigma}_{\rm F}  + n_{\rm S} \cdot \boldsymbol{\sigma}_{\rm S}  = 0 &\text{ on } \Gamma_{\rm FS}.\\
# \end{aligned}
# \right.
# ```
#
# Where $\boldsymbol{\sigma}_{\rm F}(u_{\rm F},p_{\rm F})=2\mu_{\rm F}\boldsymbol{\varepsilon}(u_{\rm F}) - p_{\rm F}\mathbf{I}$ and $\boldsymbol{\sigma}_{\rm S}(u_{\rm S})=2\mu_{\rm S}\boldsymbol{\varepsilon}(u_{\rm S}) +\lambda_{\rm S}tr(\boldsymbol{\varepsilon}(u_{\rm S}))\mathbf{I}$.

# ```@raw HTML
# <a name="geometry"></a>
# ```
# ### Geometry and Discrete model
# In this tutorial we solve the benchmark descrived in [1], consisting on a flow over an elastic flag after a cylinder.
# The computational domain is defined by a channel of size $\Omega \doteq (0,4.5)\times(0,0.41)$, with an embedded cylinder of radius $R=0.05$ and center at $C=(0.2,0.2)$. The associated FE triangulation is denoted by $\mathcal{T}$, the fluid and solid domain and their associated triangulations will be denoted by $\Omega_{\rm F}$, $\Omega_{\rm S}$, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$, respectively.
#
# In order to load the discrete model we first setup Gridap
using Gridap

# The discrete model for the elastic flag problem is generated by loading the `"../models/elasticFlag.json"` file.
model = DiscreteModelFromFile("../models/elasticFlag.json")

# We can inspect the loaded geometry and associated parts by printing to a `vtk` file:
writevtk(model,"model")

# This will produce an output in which we can identify the different parts of the domain, with the associated labels and tags.
#
# | Part | Notation | Label | Tag |
# | :---: | :---:| :---: | :---: |
# | Solid-cylinder wall | $\Gamma_{\rm S,D_{cyl}}$ | "fixed" | 1 |
# | Fluid-solid interface | $\Gamma_{\rm FS}$ | "interface" | 2 |
# | Channel inlet | $\Gamma_{\rm F,D_{in}}$ | "inlet" | 3 |
# | Channel outlet | $\Gamma_{\rm F,N_{out}}$ | "outlet" | 4 |
# | Channel walls | $\Gamma_{\rm F,D_{wall}}$ | "noSlip" | 5 |
# | Fluid-cylinder wall | $\Gamma_{\rm F,D_{cyl}}$ | "cylinder" | 6 |
# | Fluid domain | $\Omega_{\rm F}$ | "fluid" | 7 |
# | Solid domain | $\Omega_{\rm S}$ | "solid" | 8 |
#
# ![](../assets/fsi/tags.png)

# ```@raw HTML
# <a name="conditions"></a>
# ```
# ### External conditions and properties
# #### Boundary conditions
# We apply Dirichlet boundary conditions at the channel inlet, upper and lower boundaries and on the cylinder. A parabolic profile is enforced at the channel inlet, while a no-slip condition is imposed on the remaining Dirichlet boundaries.
# ```math
# \left\lbrace
# \begin{aligned}
# u_{\rm F,in}(x,y) = 1.5U\frac{y(H −y)}{\left(\frac{H}{2}\right)^2}\quad&\textrm{on }\Gamma_{\rm F,D_{in}},\\
# u_{\rm F,0}(x,y) = [0, 0]\quad&\textrm{on }\Gamma_{\rm F,D_{wall}}\cup\Gamma_{\rm F,D_{cyl}},\\
# u_{\rm S,0}(x,y) = [0, 0]\quad&\textrm{on }\Gamma_{\rm S,D_{cyl}}.\\
# \end{aligned}
# \right.
# ```

const U = 1.0
const H = 0.41
uf_in(x) = VectorValue( 1.5 * U * x[2] * ( H - x[2] ) / ( (H/2)^2 ), 0.0 )
uf_0(x) = VectorValue( 0.0, 0.0 )
us_0(x) = VectorValue( 0.0, 0.0 )

# We consider a free tranction condition at the channel outlet
# ```math
# n_{\rm F} \cdot \boldsymbol{\sigma}_{\rm F}  = \mathbf{0}\quad\textrm{on }\Gamma_{\rm F,N}
# ```
hN(x) = VectorValue( 0.0, 0.0 )
p_jump(x) = 0.0

# ```@raw HTML
# <a name="forces"></a>
# ```
# #### External forces
# In this test, the body forces acting on the fluid an solid are zero.
f(x) = VectorValue( 0.0, 0.0 )
s(x) = VectorValue( 0.0, 0.0 )
g(x) = 0.0

# ```@raw HTML
# <a name="properties"></a>
# ```
# #### Material properties
# We use a linear elastic constitutive law for the elastic flag. Given the Young's modulus $E$ and the Poisson ratio $\nu$, we can compute the Lamé constants, $\lambda$ and $\mu$, using the following function:
function lame_parameters(E,ν)
  λ = (E*ν)/((1+ν)*(1-2*ν))
  μ = E/(2*(1+ν))
  (λ, μ)
end

# Then, we get the Lamé parameters for a solid with $E=1.0$ MPa and $\nu=0.33$.
const E_s = 1.0
const ν_s = 0.33
const (λ_s,μ_s) = lame_parameters(E_s,ν_s)

# The Cauchy stress tensor for the solid part is defined by $\sigma_s = 2\mu\varepsilon + \lambda tr(\varepsilon)\mathbf{I}$. Note that we use the trace operator from the `LinearAlgebra` package. Note that this function will be used as a composition (`∘`), using as argument a function whose arguments depend on the coordinates, without the need of passing such coordinates as an argument. That is `σ_s(ε(u)) = σ_s ∘ ε(u)`.
using LinearAlgebra: tr
σₛ(ε) = λ_s*tr(ε)*one(ε) + 2*μ_s*ε

# For the fluid part, we only need to define the viscosity $\mu_f$, which we set to $0.01$.
const μ_f = 0.01

# The Cauchy stress tensor for the fluid part is given by $\sigma_f = \sigma^{dev}_f - p\mathbf{I}$, with $\sigma^{dev}_f=2\mu_f$ the deviatoric part of the stress. Since we use a mixed form with the pressure $p$ as an unknown, the stress law will only involve the deviatoric part.
σ_dev_f(ε) = 2*μ_f*ε

# ```@raw HTML
# <a name="numericalScheme"></a>
# ```
# ## Numerical scheme
# ```@raw HTML
# <a name="feSpace"></a>
# ```
# ### FE Spaces
# In this tutorial we use an *inf-sup* stable velocity-pressure pair, $P_k/P_{k-1}$ elements, with continuous velocities and pressures. We select the same velocity interpolation space for the fluid and the solid parts, defined as follows
# ```math
# V_{\rm X} \doteq \{ v \in [H^1(\Omega_{\rm X})]^d:\ v|_T\in [P_k(T)]^d \text{ for all } T\in\mathcal{T}_{\rm X} \},
# ```
# where $T$ denotes an arbitrary cell of the FE mesh $\mathcal{T}_{\rm X}$, and $P_k(T)$ is the usual Lagrangian FE space of order $k$ defined on a mesh of triangles or tetrahedra. Note that the sub-index $(\cdot)_{\rm X}$ stands for the fluid or solid parts, $(\cdot)_{\rm F}$ or $(\cdot)_{\rm S}$, respectively.
# On the other hand, the space for the pressure is only defined in the fluid domain, $\Omega_{\rm F}$, and is given by
# ```math
# Q \doteq \{ q \in C^0(\Omega):\ q|_T\in P_{k-1}(T) \text{ for all } T\in\mathcal{T}_{\rm F}\}.
# ```

# Before creating the FE spaces, we first need to define the triangulation of each of the sub-domains, which are constructed restricting the global triangulation to the corresponding part. This is done by calling the `Triangulation` (or equivalently `Interior`) function with the desired geometrical part label, i.e. `"solid"` and `"fluid"`, respectively. Here we create the triangulation of the global domain, $\mathcal{T}$, and the solid and fluid triangulations, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$.
Ω = Interior(model)
Ω_s = Interior(model,tags="solid")
Ω_f = Interior(model,tags="fluid")

# We also generate the triangulation and associated outer normal field for the outlet boundary, $\Gamma_{\rm F,N_{out}}$, which will be used to impose a Neumann condition.
Γ_out = BoundaryTriangulation(model,tags="outlet")
n_Γout = get_normal_vector(Γ_out)

# Finally, to impose the interface condition between solid and fluid, we will need the triangulation and normal field of such interface, $\Gamma_{\rm FS}$. The interface triangulation is generated by calling the `InterfaceTriangulation` function specifying the two interfacing domain triangulations. Note that the normal field will point outwards with respect to the first entry.
Γ_fs = InterfaceTriangulation(Ω_f,Ω_s)
n_Γfs = get_normal_vector(Γ_fs)

# In what follows we will assume a second-order veloticty interpolation, i.e. $k=2$
k = 2

# Now we define the reference FE for the velocity and pressure fields. The velocity field reference FE, both for fluid and solid domains, will be defined by a 2-dimensional `VectorValue` type `:Lagrangian` reference FE element of order `k`. The reference FE for the pressure will be given by a scalar value type `:Lagrangian` reference FE element of order `k-1`.
reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},k)
reffeₚ = ReferenceFE(lagrangian,Float64,k-1)

# Having set up all the ingredients, we can create the different FE spaces for the test functions. For the velocity FE spaces we call the `TestFESpace` function with the corresponding discrete model, using the velocity reference FE `reffeᵤ` and conformity `:H1`. Note that we assign different Dirichlet boundary labels for the two different parts, generating the variational spaces with homogeneous Dirichlet boundary conditions, $V_{\rm F,0}$ and $V_{\rm S,0}$ .
Vf = TestFESpace(
  Ω_f,
  reffeᵤ,
  conformity=:H1,
  dirichlet_tags=["inlet", "noSlip", "cylinder"])

Vs = TestFESpace(
  Ω_s,
  reffeᵤ,
  conformity=:H1,
  dirichlet_tags=["fixed"])

# For the pressure test FE space, we use the fluid discrete model, the pressure reference FE `reffeₚ` and `:C0` conformity.
Qf = TestFESpace(
  Ω_f,
  reffeₚ,
  conformity=:C0)

# The trial FE spaces are generated from the test FE spaces, adding the corresponding function for the various Dirichlet boundaries, leading to $U_{\rm F,g_{\rm F}}$, $U_{\rm S,g_{\rm S}}$ and $P_{\rm F}$.
Uf = TrialFESpace(Vf,[uf_in, uf_0, uf_0])
Us = TrialFESpace(Vs,[us_0])
Pf = TrialFESpace(Qf)

# Finally, we glue the test and trial FE spaces together, defining a unique test and trial space for all the fields using the `MultiFieldFESpace` function. That is $Y=[V_{\rm S,0}, V_{\rm F,0}, Q_{\rm F}]^T$ and $X=[U_{\rm S,g_{\rm S}}, U_{\rm F,g_{\rm F}}, P_{\rm F}]^T$
Y = MultiFieldFESpace([Vs,Vf,Qf])
X = MultiFieldFESpace([Us,Uf,Pf])

# ```@raw HTML
# <a name="integration"></a>
# ```
# ### Numerical integration
# Once we have all the triangulations, we can generate the quadrature rules to be applied each domain. This will be generated by calling the `Measure` function, that given a triangulation and an integration degree, it returns the Lebesgue integral measure $d\Omega$.
degree = 2*k
dΩₛ = Measure(Ω_s,degree)
dΩ_f = Measure(Ω_f,degree)

# Idem for the boundary part.
bdegree = 2*k
dΓ_out = Measure(Γ_out,bdegree)

# Idem for the interface part.
idegree = 2*k
dΓ_fs = Measure(Γ_fs,idegree)

# ```@raw HTML
# <a name="weakForm"></a>
# ```
# ### Weak form
# We now introduce the solution and test function vectors as $[\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}]^T$ and $[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}]^T$. The weak form of the coupled FSI problem using the Nitche's method, see for instance [2], reads:
#
# find $[\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}]^T \in X$ such that
# ```math
# a([\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}],[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])=l([\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])\qquad\forall[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}]^T\in Y,
# ```
# where
# ```math
# \begin{aligned}
# &a([\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}],[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])\doteq a_s(\mathbf{u}^h_{\rm s},\mathbf{v}^h_{\rm s}) + a_f((\mathbf{u}^h_{\rm f},p^h_{\rm f}),(\mathbf{v}^h_{\rm f},q^h_{\rm f})) + a_{fs}((\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm f}, p^h_{\rm f}),(\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}))\\
# &l([\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])\doteq l_s(\mathbf{v}^h_{\rm s}) + l_f(\mathbf{v}^h_{\rm f}) + l_{f,\Gamma_N}(\mathbf{v}^h_{\rm f}).
# \end{aligned}
# ```
a((us,uf,p),(vs,vf,q)) = a_s(us,vs) + a_f((uf,p),(vf,q)) + a_fs((us,uf,p),(vs,vf,q))
l((vs,vf,q)) = l_s(vs) + l_f((vf,q)) + l_f_Γn(vf)

# with the following definitions:
#
# - The bilinear form associated with the solid counterpart, $a_s(\mathbf{u}^h_{\rm s},\mathbf{v}^h_{\rm s})$, defined as
# ```math
# a_s(\mathbf{u},\mathbf{v})=\int\varepsilon(\mathbf{v}):\boldsymbol{\sigma}_s(\epsilon(\mathbf{u}))\ d\Omega_{\rm S}.
# ```
a_s(u,v) = ∫( ε(v) ⊙ (σₛ ∘ ε(u)) )dΩₛ

# - The bilinear form associated with the fluid counterpart, $a_f((\mathbf{u}^h_{\rm f},p^h_{\rm f}),(\mathbf{v}^h_{\rm f},q^h_{\rm f}))$, is defined as
# ```math
# a_f((\mathbf{u},p),(\mathbf{v},q))=\int\left(\varepsilon(\mathbf{v}):\boldsymbol{\sigma}^{dev}_f(\mathbf{u}) - (\nabla\cdot\mathbf{v})\ p + q\ (\nabla\cdot\mathbf{u})\right) d\Omega_{\rm F}.
# ```
a_f((u,p),(v,q)) = ∫( ε(v) ⊙ (σ_dev_f ∘ ε(u)) - (∇⋅v)*p + q*(∇⋅u) )dΩ_f

# - The bilinear form associated with the coupling between fluid and solid counterparts, $a_{fs}((\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm f}, p^h_{\rm f}),(\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}))$, is given by:
# ```math
# \begin{aligned}
# a_{fs}((\mathbf{u}_{\rm s},\mathbf{u}_{\rm f}, p),((\mathbf{v}_{\rm s},\mathbf{v}_{\rm f}, q))=&\int\left(\gamma\frac{\mu_f}{h}(\mathbf{v}_{\rm f}-\mathbf{v}_{\rm s})\cdot(\mathbf{u}_{\rm f}-\mathbf{u}_{\rm s})\right.\\
# &-  (\mathbf{v}_{\rm f}-\mathbf{v}_{\rm s})\cdot(\mathbf{n}_{\rm fs}\cdot\boldsymbol{\sigma}^{dev}_f(\mathbf{u}_{\rm f}))
# +   (\mathbf{v}_{\rm f}-\mathbf{v}_{\rm S})\cdot(p_{\rm f}\mathbf{n}_{\rm fs})\\
# &-  \chi(\mathbf{n}_{\rm fs}\cdot\boldsymbol{\sigma}^{dev}_f(\mathbf{v}_{\rm f})\cdot(\mathbf{u}^{\rm f}-\mathbf{u}_{\rm s})
# + \left.(q_{\rm f}\mathbf{n}_{\rm fs})\cdot(\mathbf{u}_{\rm f}-\mathbf{u}_{\rm s})\ \right) d\Gamma_{\rm FS}
# \end{aligned}
# ```
a_fs((us,uf,p),(vs,vf,q)) = ∫( α*(jump_Γ(vf,vs)⋅jump_Γ(uf,us)) +
                               jump_Γ(vf,vs)⋅t_Γfs(1,p,uf,n_Γfs) +
                               t_Γfs(χ,q,vf,n_Γfs)⋅jump_Γ(uf,us) )dΓ_fs
# Where $\chi$ is a parameter that can take values $1.0$ or $-1.0$ and it is used to define the symmetric or antisymmetric version of the method, respectively. To difine this form we used the well known Nitsche's method, which enforces the continuity of fluid and solid velocities as well as the continuity of the normal stresses, see for instance [2].
const γ = 1.0
const χ = -1.0
jump_Γ(wf,ws) = wf.⁺-ws.⁻
t_Γfs(χ,q,w,n) = q.⁺*n.⁺ - χ*n.⁺⋅(σ_dev_f∘ε(w.⁺))

# From the interface triangulation we can obtain the interface elements length, $h$, and the penalty parameter, $\alpha=\gamma\frac{\mu_f}{h}$, used in the Nitsche's terms.
using Gridap.CellData
dim = num_cell_dims(Γ_fs)
h_Γfs = get_array(∫(1)dΓ_fs)
α = CellField( lazy_map(h->γ*μ_f/(h.^dim),h_Γfs), Γ_fs)

# - The linear form associated with the solid counterpart, $l_s(\mathbf{v}^h_{\rm s})$, is defined as
# ```math
# l_s(\mathbf{v})=\int\mathbf{v}\cdot s\ d\Omega_{\rm S}.
# ```
l_s(v) = ∫( v⋅s )dΩₛ

# - The linear form associated with the fluid counterpart, $l_f(\mathbf{v}^h_{\rm f})$ is defined as
# ```math
# l_f(\mathbf{v},q)=\int\mathbf{v}\cdot f + qg\ d\Omega_{\rm F}.
# ```
l_f((v,q)) = ∫( v⋅f + q*g )dΩ_f

# - The linear form associated with the fluid Neumann boundary condition, $l_{f,\Gamma_N}(\mathbf{v}^h_{\rm f})$, is defined as
# ```math
# l_{f,\Gamma_N}(\mathbf{v})=\int\mathbf{v}\cdot h\ d\Gamma_{\rm N_{out}}.
# ```
l_f_Γn(v) = ∫( v⋅hN )dΓ_out

# The final bilinear and linear forms will be given by

# ```@raw HTML
# <a name="algebraic"></a>
# ```
# ### Algebraic System of Equations
# After defining the weak form of the problem and the integration quadrature rules to perform the numerical integration, we are ready to assemble the linear system of equations. In this case, the system will have the following structure:
# ```math
# \begin{bmatrix}
# \mathbf{A}_{u,\rm S}&\mathbf{A}_{u,\rm SF}&0\\
# \mathbf{A}_{u,\rm FS}&\mathbf{A}_{u,\rm F}&\mathbf{A}_{up,\rm F}\\
# 0&\mathbf{A}_{up,\rm F}&\mathbf{A}_{p,\rm F}\\
# \end{bmatrix}\begin{bmatrix}
# \mathbf{U}^h_{\rm S}\\
# \mathbf{U}^h_{\rm F}\\
# \mathbf{P}^h_{\rm F}\\
# \end{bmatrix} = \begin{bmatrix}
# \mathbf{F}^h_{u,\rm S}\\
# \mathbf{F}^h_{u,\rm F}\\
# \mathbf{F}^h_{p,\rm F}\\
# \end{bmatrix}
# ```
# In order to construct the system we define the final FE operator, constructed using the function `AffineFEOperator` passing as arguments the bilinear and linear forms, $a$ and $l$, together with the trial and test FE spaces, $X$ and $Y$.
op = AffineFEOperator(a,l,X,Y)

# Finally, we call `solve` to obtain the solution vector of nodal values $[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$
uhs, uhf, ph = solve(op)

# ```@raw HTML
# <a name="postprocess"></a>
# ```
# ## Post-processing
# ```@raw HTML
# <a name="visualization"></a>
# ```
# ### Visualization
# The solution fields $[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$ are defined over all the domain, extended with zeros on the inactive part. Calling the function `writevtk` passing the global triangulation, we will output the global fields.
writevtk(Ω,"trian", cellfields=["uhs" => uhs, "uhf" => uhf, "ph" => ph])
# ![](../assets/fsi/Global_solution.png)

# However, we can also restrict the fields to the active part by calling the function `restrict` with the field along with the respective active triangulation.
writevtk(Ω_s,"trian_solid",cellfields=["uhs"=>uhs])
writevtk(Ω_f,"trian_fluid",cellfields=["ph"=>ph,"uhf"=>uhf])
# ![](../assets/fsi/Local_solution.png)

# ```@raw HTML
# <a name="QOIs"></a>
# ```
# ### Quantities of Interest
Γ_S = BoundaryTriangulation(model,tags=["cylinder","interface"])
dΓₛ = Measure(Γ_S,bdegree)
n_ΓS = get_normal_vector(Γ_S)
FD, FL = sum( ∫( (n_ΓS⋅σ_dev_f(ε(uhf))) - ph*n_ΓS )*dΓₛ )
println("Drag force: ", FD)
println("Lift force: ", FL)

# ## References
# [1] Turek, S., Hron, J., Madlik, M., Razzaq, M., Wobker, H., & Acker, J. F. (2011).* Numerical simulation and benchmarking of a monolithic multigrid solver for fluid-structure interaction problems with application to hemodynamics*. In Fluid Structure Interaction II (pp. 193-220). Springer, Berlin, Heidelberg.*
#
# [2] Burman, E., and Fernández, M. A. *Stabilized explicit coupling for fluid–structure interaction using Nitsche's method.* Comptes Rendus Mathematique 345.8 (2007): 467-472.


text:
using SyntheticImplicitFeedback
using Base.Test

# define a set of features
features = Feature[]

push!(features, Feature("Age", 1950, () -> rand(1930:2010)))
push!(features, Feature("Geo", "Arizona", () -> rand(["Arizona", "California", "Colorado", "Illinois", "Indiana", "Michigan", "New York", "Utah"])))
push!(features, Feature("Ad", 0, () -> rand(0:4)))

# create a set of rules
rules = Rule[]

push!(rules, Rule(s -> true, 0.001))
push!(rules, Rule(s -> s["Ad"] == 2, 0.01))
push!(rules, Rule(s -> s["Age"] >= 1980 && s["Age"] <= 1989 && s["Geo"] == "New York" && s["Ad"] == 0, 0.30))
push!(rules, Rule(s -> s["Age"] >= 1950 && s["Age"] <= 1959 && s["Geo"] == "New York" && s["Ad"] == 1, 0.30))
push!(rules, Rule(s -> s["Age"] >= 1980 && s["Age"] <= 1989 && s["Geo"] == "Arizona" && s["Ad"] == 1, 0.30))
push!(rules, Rule(s -> s["Age"] >= 1950 && s["Age"] <= 1959 && s["Geo"] == "Arizona" && s["Ad"] == 0, 0.30))

# check if a accumulative CTR is computed correctly
sample1 = Dict("Age" => 1940, "Geo" => "California", "Ad" => 0)
@test accumulate(sample1, rules) == 0.001
@test isa(generate(sample1, rules), Bool)

sample2 = Dict("Age" => 1940, "Geo" => "California", "Ad" => 2)
@test accumulate(sample2, rules) == 0.011
@test isa(generate(sample2, rules), Bool)

sample3 = Dict("Age" => 1953, "Geo" => "New York", "Ad" => 1)
@test accumulate(sample3, rules) == 0.301
@test isa(generate(sample3, rules), Bool)

# generate samples with random pairs of demographics and ad variants
samples = Dict[]

for i in 1:100
    sample = Dict()
    for f in features
        sample[f.name] = f.random()
    end
    push!(samples, sample)
end

feedback = generate(samples, rules)

@test isa(feedback, Array{Bool})
@test length(feedback) == length(samples)

ext: jl


text:
##=================####==============##
## SCRIPT TO PRODUCE FIGURES
##=================####==============##

##=================####==============##

# using CSV
using Plots          # plotting
using LaTeXStrings   # latex labels in plots
using DelimitedFiles # latex labels in plots

##=================####==============##
# PATH TO SAVE IMGS (MODIFY WITH YOUR OWN!)

top_path = joinpath(homedir(), "Code", "Collaborations", "networkEpidemiology_MX", "julia_ABM")

data_path = joinpath(top_path, "data", "EOD")
img_out_path = joinpath(top_path, "figures", "EOD")

##=================####==============##
time_lockdown = 15
HC_capacity = 0.25

files = readdir(data_path)
pop_files = filter(x->occursin("populations", x), files)
time_files = filter(x->occursin("infTime", x), files)

cop_ag_vals = map(x->parse(Float64, match(r"_CA_(\d\.\d)_", x).captures[1]), pop_files)

num_agents = parse(Int, match(r"_N_(\d+)_", pop_files[1]).captures[1])
attack_rate = parse(Float64, match(r"_AR_(\d+\.\d+)_", pop_files[1]).captures[1])

##=================####==============##
# FRACTION OF INFECTED AGENTS AS A FUNTCION OF FRACTION OF LOCKDOWN AGENTS

plt = plot(
xlabel = L"t",
ylabel = L"N_I/N",
guidefontsize = 16,
tickfontsize = 10,
colorbar = :none,
legendtitle = L"p_C",
xlims = (0, 70)
)

max_I = Dict()

# for file in pop_files
for i in 1:length(pop_files)
    # data = readdlm(joinpath(data_path, file), ',')
    # cop_val = parse(Float64, match(r"_CA_(\d\.\d)_", file).captures[1])
    # println(file, "|", cop_val)

    data = readdlm(joinpath(data_path, pop_files[i]), ',')
    cop_val = parse(Float64, match(r"_CA_(\d\.\d)_", pop_files[i]).captures[1])
    println(pop_files[i], "|", cop_val)

    t_vals = data[2:end, 1]
    num_I = data[2:end, 3]

    max_I[cop_val] = maximum(num_I)

    plot!(plt,
        t_vals,
        num_I,
        line_z = cop_ag_vals[i],
        lab = cop_ag_vals[i],
        # line_z = cop_val,
        # lab = cop_val,
        lw = 1.25,
        color = :viridis,
        # ylims = (0.0, 0.2)
        )

end

annotate!(plt, 15, 0.61, text("lockdown", :red, :right, 10))
vline!(plt, [time_lockdown], ls = :dash, color = :red, lab = "")

# annotate!(plt, 50, 0.3, text("HC capacity", :blue, :left, 10))
# hline!(plt, [HC_capacity], ls = :dash, color = :blue, lab = "")

display(plt)

##=================####==============##

savefig(plt, joinpath(img_out_path, "SIR_N_$(num_agents)_AR_$(attack_rate)_CA.png"))

##=================####==============##
# FRACTION OF INFECTED AGENTS AS A FUNTCION OF FRACTION OF LOCKDOWN AGENTS

plt = plot(
xlabel = L"t",
ylabel = L"N_S/N",
guidefontsize = 16,
tickfontsize = 10,
colorbar = :none,
legendtitle = L"p_C",
xlims = (0, 70)
)

max_S = Dict()

# for file in pop_files
for i in 1:length(pop_files)
    # data = readdlm(joinpath(data_path, file), ',')
    # cop_val = parse(Float64, match(r"_CA_(\d\.\d)_", file).captures[1])
    # println(file, "|", cop_val)

    data = readdlm(joinpath(data_path, pop_files[i]), ',')
    cop_val = parse(Float64, match(r"_CA_(\d\.\d)_", pop_files[i]).captures[1])
    println(pop_files[i], "|", cop_val)

    t_vals = data[2:end, 1]
    num_I = data[2:end, 2]

    max_S[cop_val] = num_I[end]

    plot!(plt,
        t_vals,
        num_I,
        line_z = cop_ag_vals[i],
        lab = cop_ag_vals[i],
        # line_z = cop_val,
        # lab = cop_val,
        lw = 1.25,
        color = :viridis,
        # ylims = (0.0, 0.2)
        )

end

annotate!(plt, 15, 0.61, text("lockdown", :red, :right, 10))
vline!(plt, [time_lockdown], ls = :dash, color = :red, lab = "")

# annotate!(plt, 50, 0.3, text("HC capacity", :blue, :left, 10))
# hline!(plt, [HC_capacity], ls = :dash, color = :blue, lab = "")

display(plt)

##=================####==============##

savefig(plt, joinpath(img_out_path, "S_N_$(num_agents)_AR_$(attack_rate)_CA.png"))

##=================####==============##
# FRACTION OF INFECTED AGENTS AS A FUNTCION OF FRACTION OF LOCKDOWN AGENTS

plt = plot(
xlabel = L"\tau_I",
ylabel = L"N(\tau_I)",
guidefontsize = 16,
tickfontsize = 10,
colorbar = :none,
legendtitle = L"p_C",
xlims = (0, 70)
)

max_S = Dict()

# for file in pop_files
for i in 1:length(time_files)
    # data = readdlm(joinpath(data_path, file), ',')
    # cop_val = parse(Float64, match(r"_CA_(\d\.\d)_", file).captures[1])
    # println(file, "|", cop_val)

    data = readdlm(joinpath(data_path, time_files[i]), ',')
    cop_val = parse(Float64, match(r"_CA_(\d\.\d)_", time_files[i]).captures[1])
    println(time_files[i], "|", cop_val)

    t_vals = data[2:end, 1]
    num_I = filter(x -> x > 0, data[2:end, 2])

    max_S[cop_val] = num_I[end]

    histogram!(plt,
        num_I,
        # line_z = cop_ag_vals[i],
        # lab = cop_ag_vals[i],
        fill_z = cop_val,
        lab = cop_val,
        # lw = 1.25,s
        color = :viridis,
        # normalized = true
        xlims = (15, 50)
        )

end

display(plt)

##=================####==============##

savefig(plt, joinpath(img_out_path, "infT_$(num_agents)_AR_$(attack_rate)_CA.png"))

##=================####==============##

# FRACTION OF INFECTED AGENTS AS A FUNTCION OF FRACTION OF LOCKDOWN AGENTS

plt = plot(
xlabel = L"p_C",
ylabel = L"\max(N_\star/N)",
guidefontsize = 16,
tickfontsize = 10,
leg = :bottomright
)

pC = sort(collect(keys(max_I)))

plot!(plt,
    pC,
    [max_I[x] for x in pC],
    m = :o,
    label = L"I"
    )

# plot!(plt,
#     pC,
#     [max_S[x] for x in pC],
#     m = :o,
#     label = L"S"
#     )

# annotate!(plt, 15, 0.61, text("lockdown", :red, :right, 10))
# vline!(plt, [time_lockdown], ls = :dash, color = :red, lab = "")

# annotate!(plt, 50, 0.3, text("HC capacity", :blue, :left, 10))
# hline!(plt, [HC_capacity], ls = :dash, color = :blue, lab = "")

display(plt)

##=================####==============##

savefig(plt, joinpath(img_out_path, "poblaciones_asint_N_$(num_agents)_AR_$(attack_rate)_CA.png"))

##=================####==============##


text:
module test_Variables

import ModiaMath

# Desired:
#   using Test
#
# In order that Test needs not to be defined in the user environment, it is included via ModiaMath:
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using ModiaMath.Test
end

mutable struct Revolute <: ModiaMath.AbstractComponentWithVariables
    _internal::ModiaMath.ComponentInternal
    phi::ModiaMath.RealScalar
    w::ModiaMath.RealScalar
    a::ModiaMath.RealScalar
    tau::ModiaMath.RealScalar
    drive::Bool

    function Revolute(;phi0::Float64=0.0, w0::Float64=0.0, drive::Bool=false)
        this = new(ModiaMath.ComponentInternal(:Revolute))
        phi = ModiaMath.RealScalar(:phi, this, start=phi0, unit="rad",     fixed=true, info="Relative rotation angle",                     numericType=drive ? ModiaMath.WR : ModiaMath.XD_EXP)
        w   = ModiaMath.RealScalar("w", this, start=w0, unit="rad/s",   fixed=true, info="Relative angular velocity",     integral=phi, numericType=drive ? ModiaMath.WR : ModiaMath.XD_EXP,     analysis=ModiaMath.OnlyDynamicAnalysis)
        a   = ModiaMath.RealScalar("a", this, start=0.0, unit="rad/s^2",             info="Relative angular acceleration", integral=w, numericType=drive ? ModiaMath.WR : ModiaMath.DER_XD_EXP, analysis=ModiaMath.OnlyDynamicAnalysis)
        tau = ModiaMath.RealScalar(:tau, this, start=0.0, unit="N*m",                 info="Driving torque",                              numericType=ModiaMath.WR,                                analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
        return this
    end
end

function Base.show(io::IO, rev::Revolute)
    print("Revolute(",
         "\n    phi = ", rev.phi,
         "\n    w   = ", rev.w,
         "\n    a   = ", rev.a,
         "\n    tau = ", rev.tau,
         "\n   )")
end

mutable struct Frame <: ModiaMath.AbstractComponentWithVariables
    _internal::ModiaMath.ComponentInternal
    r::ModiaMath.RealSVector3
    q::ModiaMath.RealSVector{4}

    derq::ModiaMath.RealSVector{4}
    v::ModiaMath.RealSVector3
    w::ModiaMath.RealSVector3

    a::ModiaMath.RealSVector3
    z::ModiaMath.RealSVector3

    f::ModiaMath.RealSVector3
    t::ModiaMath.RealSVector3

    residue_w::ModiaMath.RealSVector3
    residue_f::ModiaMath.RealSVector3
    residue_t::ModiaMath.RealSVector3
    residue_q::ModiaMath.RealScalar

    drive::Bool

    function Frame(;r0=zeros(3), q0=[0,0,0,1], v0=zeros(3), w0=zeros(3), drive::Bool=false)
        this    = new(ModiaMath.ComponentInternal(:Frame))

        r         = ModiaMath.RealSVector3(:r, this, start=r0, unit="m",      fixed=true, info="Relative position",                         numericType=drive ? ModiaMath.WR : ModiaMath.XD_EXP)
        q         = ModiaMath.RealSVector{4}(:q, this, start=q0,              fixed=true, info="Relative quaternion",                       numericType=drive ? ModiaMath.WR : ModiaMath.XD_IMP)

        derq      = ModiaMath.RealSVector{4}(:derq, this,      unit="1/s",                info="der(q)",                        integral=q, numericType=drive ? ModiaMath.WR : ModiaMath.DER_XD_IMP, analysis=ModiaMath.OnlyDynamicAnalysis)
        v         = ModiaMath.RealSVector3(:v, this, start=v0, unit="m/s",    fixed=true, info="Relative velocity",             integral=r, numericType=drive ? ModiaMath.WR : ModiaMath.XD_IMP,     analysis=ModiaMath.OnlyDynamicAnalysis)
        w         = ModiaMath.RealSVector3(:w, this, start=w0, unit="rad/s",  fixed=true, info="Relative angular velocity",                 numericType=drive ? ModiaMath.WR : ModiaMath.XD_IMP,     analysis=ModiaMath.OnlyDynamicAnalysis)

        a         = ModiaMath.RealSVector3(:a, this,           unit="m/s^2",              info="Relative acceleration",         integral=v, numericType=drive ? ModiaMath.WR : ModiaMath.DER_XD_IMP, analysis=ModiaMath.OnlyDynamicAnalysis)
        z         = ModiaMath.RealSVector3(:z, this,           unit="rad/s^2",            info="Relative angular acceleration", integral=w, numericType=drive ? ModiaMath.WR : ModiaMath.DER_XD_IMP, analysis=ModiaMath.OnlyDynamicAnalysis)

        f         = ModiaMath.RealSVector3(:f, this,           unit="N",                  info="Driving force",                             numericType=ModiaMath.WR,                                analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)
        t         = ModiaMath.RealSVector3(:t, this,           unit="N*m",                info="Driving torque",                            numericType=ModiaMath.WR,                                analysis=ModiaMath.QuasiStaticAndDynamicAnalysis)

        residue_w = ModiaMath.RealSVector3(:residue_w, this,                              info="Angular velocity residue",                  numericType=ModiaMath.FD_IMP,                            analysis=ModiaMath.OnlyDynamicAnalysis)
        residue_f = ModiaMath.RealSVector3(:residue_f, this,                              info="Momentum equation residue",                 numericType=ModiaMath.FD_IMP,                            analysis=ModiaMath.OnlyDynamicAnalysis)
        residue_t = ModiaMath.RealSVector3(:residue_t, this,                              info="Angular momentum equation residue",         numericType=ModiaMath.FD_IMP,                            analysis=ModiaMath.OnlyDynamicAnalysis)
        residue_q = ModiaMath.RealScalar(:residue_q, this,                                info="Quaternion constraint residue",             numericType=ModiaMath.FC,                                analysis=ModiaMath.OnlyDynamicAnalysis)

        this.drive = drive
        return this
    end
end

function Base.show(io::IO, frame::Frame)
    print("Revolute(",
         "\n    r = ", frame.r,
         "\n    q = ", frame.q,
         "\n    v = ", frame.v,
         "\n    w = ", frame.w,
         "\n    a = ", frame.a,
         "\n    z = ", frame.z,
         "\n    f = ", frame.f,
         "\n    t = ", frame.t,
         "\n   )")
end

ModiaMath.@component Robot(;phi10=0.0, phi20=0.0, var10=0.0, r0=zeros(3), q0=zeros(4)) begin
    rev1  = Revolute(phi0=phi10)
    rev2  = Revolute(phi0=phi20)
    var1  = ModiaMath.RealScalar(start=var10, numericType=ModiaMath.XA)
    res1  = ModiaMath.RealScalar(numericType=ModiaMath.FD_IMP)
    frame = Frame(r0=r0, q0=q0)
end

ModiaMath.@component Robot2(;phi10=0.0, phi20=0.0, r0=zeros(3), q0=[0.0, 0.0, 0.0, 1.0]) begin
    rev1  = Revolute(phi0=phi10, drive=true)
    rev2  = Revolute(phi0=phi20, drive=true)
    frame = Frame(r0=r0, q0=q0,  drive=true)
end

ModiaMath.@component Robot3(;phi10=0.0, phi20=0.0, phi30=-2.0) begin
    rev1 = Revolute(phi0=phi10, drive=true)
    rev2 = Revolute(phi0=phi20, drive=true)
    rev3 = Revolute(phi0=phi30)
    res1 = ModiaMath.RealScalar(numericType=ModiaMath.FC)
end

@testset "ModiaMath: test_Variables.jl" begin

    @testset "Dynamic analysis" begin # ------------------------------ Dynamic analysis
        r0 = [1.0,2.0,3.0]
        q0 = [0.5,0.5,0.0,sqrt(0.5^2 + 0.5^2)]
        robot = Robot(phi10=1.0, phi20=2.0, var10=3.0, r0=r0, q0=q0)
        robot.rev1.a.value = 10 * 2.22
        robot.rev2.a.value = 10 * 4.44

        println("\n... robot = ", robot)

        println("\n... Print variables of robot")
        m = ModiaMath.ModelVariables(robot)
        ModiaMath.print_ModelVariables(m)

        println("\n... Copy start values to x")
        x       = zeros(5 + 7 + 6)
        x_fixed = fill(false, 5 + 7 + 6)
        ModiaMath.copy_start_to_x!(m, x, x_fixed)
        x0 = [1.0, 0.0, 2.0, 0.0, 1.0, 2.0, 3.0, 0.5, 0.5, 0.0, sqrt(0.5^2 + 0.5^2), 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3.0 ]
        x0_fixed = fill(true, 17)
        append!(x0_fixed, false)
        @test isapprox(x, x0)
        @test x_fixed == x0_fixed

        println("\n... Copy x and der_x to variables")
        x    = [1.11, 2.22 , 3.33, 4.44, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 1.1, 2.2, 3.3, 4.4, 5.5, 6.6 , 5.55]
        derx = [2.22, 22.2, 4.44, 44.4, 1.1, 2.2, 3.3, 4.44, 5.55, 6.66, 7.77, 1.11, 2.22, 3.33, 4.44, 5.55, 6.66, 0.0]
        ModiaMath.copy_x_and_derx_to_variables!(0.5, x, derx, m)
        @test isapprox(x, [robot.rev1.phi.value, robot.rev1.w.value, robot.rev2.phi.value, robot.rev2.w.value,
                            robot.frame.r.value[1], robot.frame.r.value[2], robot.frame.r.value[3],
                            robot.frame.q.value[1], robot.frame.q.value[2], robot.frame.q.value[3], robot.frame.q.value[4],
                            robot.frame.v.value[1], robot.frame.v.value[2], robot.frame.v.value[3],
                            robot.frame.w.value[1], robot.frame.w.value[2], robot.frame.w.value[3],
                            robot.var1.value], rtol=1e-15)
        @test isapprox(derx[2], robot.rev1.a.value, rtol=1e-15)
        @test isapprox(derx[4], robot.rev2.a.value, rtol=1e-15)

        println("\n... Copy variables to residues")
        residues = zeros(5 + 7 + 6)
        ModiaMath.copy_variables_to_residue!(m, x, derx, residues)
        println("residue = ", residues)
        @test isapprox(residues, zeros(5 + 7 + 6), atol=1e-12)
    end

    @testset "Kinematic analysis 1" begin   # ---------------------------- Kinematic analysis 1
        robot2 = Robot2(phi10=1.0, phi20=2.0)
        println("\n... robot2 = ", robot2)

        println("\n... Print variables of robot2")
        m = ModiaMath.ModelVariables(robot2, analysis=ModiaMath.KinematicAnalysis)
        ModiaMath.print_ModelVariables(m)

        println("\n... Copy start values to x")
        x       = [1.11]
        x_fixed = [false]
        ModiaMath.copy_start_to_x!(m, x, x_fixed)
        @test isapprox(x, [0.0])
        @test x_fixed == [true]

        println("\n... Copy x and der_x to variables")
        x    = fill(2.1, 1)
        derx = fill(3.2, 1)
        ModiaMath.copy_x_and_derx_to_variables!(0.5, x, derx, m)
        @test isapprox(x, [m.var[2].value], rtol=1e-15)    # var[2] = _dummy_x

        println("\n... Copy variables to residues")
        residues = zeros(m.nx)
        ModiaMath.copy_variables_to_residue!(m, x, derx, residues)
        @test isapprox(residues, [3.2 - (-2.1)], atol=1e-12)
    end

    @testset "Kinematic analysis 2" begin   # ---------------------------- Kinematic analysis 2
        robot3 = Robot3(phi10=1.0, phi20=2.0, phi30=-2.0)
        println("\n... robot3 = ", robot3)

        println("\n... Print variables of robot3")
        m = ModiaMath.ModelVariables(robot3, analysis=ModiaMath.KinematicAnalysis)
        ModiaMath.print_ModelVariables(m)

        println("\n... Copy start values to x")
        x = zeros(m.nx)
        x_fixed = fill(false, m.nx)
        ModiaMath.copy_start_to_x!(m, x, x_fixed)
        @test isapprox(x, [-2.0])
        @test x_fixed == [true]

        println("\n... Copy x and der_x to variables")
        x    = [1.11]
        derx = [2.22]
        ModiaMath.copy_x_and_derx_to_variables!(0.5, x, derx, m)
        @test isapprox(x, [robot3.rev3.phi.value], rtol=1e-15)

        println("\n... Copy variables to residues")
        residues = zeros(m.nx)
        ModiaMath.copy_variables_to_residue!(m, x, derx, residues)
        @test isapprox(residues, zeros(1), atol=1e-12)
    end

end

end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

const VariableSynchronization = (;[
    Symbol("VARIABLE_SYNCHRONIZATION_AUTO") => Int32(0),
    Symbol("VARIABLE_SYNCHRONIZATION_NONE") => Int32(1),
    Symbol("VARIABLE_SYNCHRONIZATION_ON_WRITE") => Int32(2),
    Symbol("VARIABLE_SYNCHRONIZATION_ON_READ") => Int32(3),
]...)

const VariableAggregation = (;[
    Symbol("VARIABLE_AGGREGATION_NONE") => Int32(0),
    Symbol("VARIABLE_AGGREGATION_SUM") => Int32(1),
    Symbol("VARIABLE_AGGREGATION_MEAN") => Int32(2),
    Symbol("VARIABLE_AGGREGATION_ONLY_FIRST_REPLICA") => Int32(3),
]...)

mutable struct SaveSliceInfoDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SaveSliceInfoDef(; kwargs...)
        obj = new(meta(SaveSliceInfoDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct SaveSliceInfoDef
const __meta_SaveSliceInfoDef = Ref{ProtoMeta}()
function meta(::Type{SaveSliceInfoDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SaveSliceInfoDef)
            __meta_SaveSliceInfoDef[] = target = ProtoMeta(SaveSliceInfoDef)
            pack = Symbol[:full_shape,:var_offset,:var_shape]
            allflds = Pair{Symbol,Union{Type,String}}[:full_name => AbstractString, :full_shape => Base.Vector{Int64}, :var_offset => Base.Vector{Int64}, :var_shape => Base.Vector{Int64}]
            meta(target, SaveSliceInfoDef, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SaveSliceInfoDef[]
    end
end
function Base.getproperty(obj::SaveSliceInfoDef, name::Symbol)
    if name === :full_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :full_shape
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :var_offset
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    elseif name === :var_shape
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int64}
    else
        getfield(obj, name)
    end
end

mutable struct VariableDef <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function VariableDef(; kwargs...)
        obj = new(meta(VariableDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct VariableDef
const __meta_VariableDef = Ref{ProtoMeta}()
function meta(::Type{VariableDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_VariableDef)
            __meta_VariableDef[] = target = ProtoMeta(VariableDef)
            fnum = Int[1,6,2,3,4,5,7,8,9]
            allflds = Pair{Symbol,Union{Type,String}}[:variable_name => AbstractString, :initial_value_name => AbstractString, :initializer_name => AbstractString, :snapshot_name => AbstractString, :save_slice_info_def => SaveSliceInfoDef, :is_resource => Bool, :trainable => Bool, :synchronization => Int32, :aggregation => Int32]
            meta(target, VariableDef, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_VariableDef[]
    end
end
function Base.getproperty(obj::VariableDef, name::Symbol)
    if name === :variable_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :initial_value_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :initializer_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :snapshot_name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :save_slice_info_def
        return (obj.__protobuf_jl_internal_values[name])::SaveSliceInfoDef
    elseif name === :is_resource
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :trainable
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :synchronization
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :aggregation
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

export VariableSynchronization, VariableAggregation, VariableDef, SaveSliceInfoDef


text:
module mzXML

using LightXML, Unitful, Codecs

struct MSscan{T<:AbstractFloat,TI<:Real}
    polarity::Char
    msLevel::Int
    retentionTime::typeof(1.0u"s")
    lowMz::Float64
    highMz::Float64
    basePeakMz::Float64
    totIonCurrent::Float64
    mz::Vector{T}
    I::Vector{TI}
    children::Vector{MSscan{T,TI}}
end

const empty32 = Array{MSscan{Float32,Float32}}[]
const empty64 = Array{MSscan{Float64,Float64}}[]

function index(filename)
    scanpositions = open(filename) do file
        # Find the indexOffset element
        seekend(file)
        p = position(file)
        l = 200
        if p < l
            error("Too short to be mzXML")
        end
        skip(file, -l)
        str = loadrest(file)
        m = match(r"<indexOffset>([0-9].*)</indexOffset>", str)
        if m == nothing
            error("Cannot find indexOffset element")
        end
        length(m.captures) == 1 || error("indexOffset should contain a single number")
        indexpos = parse(Int, m.captures[1])
        length_tail = l - m.offset
        # Read the index
        seek(file, indexpos)
        str = loadrest(file)
        xindex = parse_string(str[1:end-length_tail-1])
        xroot = root(xindex)
        scanpositions = Array{Int}(0)
        for c in child_elements(xroot)
            name(c) == "offset" || error("index could not be parsed")
            push!(scanpositions, parse(Int, content(c)))
        end
        scanpositions
    end
    scanpositions
end

function load(filename; maxlevel = 1)
    xdoc = parse_file(filename)
    xroot = root(xdoc)
    if name(xroot) != "mzXML" || name(xroot) != "indexedmzML"
        error("Not an mzXML file")
    end

    # Find the msRun node
    msRun = find_element(xroot, "msRun")

    props = Dict{Symbol,Any}()
    props[:startTime] = parse_time(attribute(msRun, "startTime"))
    props[:endTime] = parse_time(attribute(msRun, "endTime"))
    el = find_element(msRun, "msInstrument")
    el = find_element(el, "msModel")
    props[:msModel] = attribute(el, "value")

    load_scans(msRun, maxlevel-1), props
end

function load_scans(elm, ndeeper=0)
    # First we discover the data type
    local T
    local TI
    for c in child_elements(elm)
        n = name(c)
        if n != "scan"
            continue
        end
        peak = find_element(c, "peaks")
        TI, T, nochildren = precisiondict[attribute(peak, "precision")]
        break
    end
    # Now load the data
    scans = Array{MSscan{T,TI}}(0)
    load_scans!(scans, elm, ndeeper)
end

function load_scans!(scans, elm, ndeeper)
    for c in child_elements(elm)
        n = name(c)
        if n != "scan"
            continue
        end
        push!(scans, load_scan(c, ndeeper)::eltype(scans))
    end
    scans
end

polaritydict = Dict("+" => '+', "-" => '-')
precisiondict = Dict("32" => (Float32, Float32, empty32), "64" => (Float64, Float64, empty64))

function load_scan(elm, ndeeper)
    polarity = polaritydict[attribute(elm, "polarity")]
    retentionTime = parse_time(attribute(elm, "retentionTime"))
    lMza, hMza = attribute(elm, "lowMz"), attribute(elm, "highMz")
    if lMza != nothing
        lowMz, highMz = parse(Float64, lMza), parse(Float64, hMza)
    else
        lowMz = highMz = NaN
    end
    msLevela = attribute(elm, "msLevel")
    msLevel = msLevela == nothing ? 1 : parse(Int, msLevela)
    basePeakMz = parse(Float64, attribute(elm, "basePeakMz"))
    totIonCurrent = parse(Float64, attribute(elm, "totIonCurrent"))
    npeaks = parse(Int, attribute(elm, "peaksCount"))
    peak = find_element(elm, "peaks")
    data = decode(Base64, content(peak))
    TI, T, nochildren = precisiondict[attribute(peak, "precision")]
    A = reinterpret(TI, data)
    bo = attribute(peak, "byteOrder")
    if bo == "network"
        ntoh!(A)
    else
        error("Don't know what to do with byteOrder $bo")
    end
    po = attribute(peak, "pairOrder")
    if po == nothing
        po = attribute(peak, "contentType")
    end
    po == "m/z-int" || error("Don't know what to do with pairOrder/contentType $po")
    I = A[2:2:end]
    mz = reinterpret(T, A[1:2:end])
    children = ndeeper > 0 ? load_scans!(Array{MSscan{T,TI}}(0), elm, ndeeper-1) : nochildren
    MSscan{T,TI}(polarity, msLevel, retentionTime, lowMz, highMz, basePeakMz, totIonCurrent, mz, I, children)
end

function loadrest(file::IOStream)
    nb = filesize(file) - position(file)
    b = read(file, nb)
    return String(b)
end

function parse_time(tstr)
    if startswith(tstr, "PT") && endswith(tstr, "S")
        return parse(Float64, tstr[3:end-1])*u"s"
    end
    error("Time string $tstr not recognized")
end

function ntoh!(A)
    for i = 1:length(A)
        A[i] = ntoh(A[i])
    end
end

end


text:
## --------  Marginal distributions

"""
	$(SIGNATURES)

Abstract type for marginal distributions of endowments.
"""
abstract type AbstractMarginal{T1} end

"""
	$(SIGNATURES)

Marginal distribution that is Uniform[0, 1].
"""
struct PercentileMarginal{T1 <: Real} <: AbstractMarginal{T1} end

"""
	$(SIGNATURES)

Uniform marginal distribution over fixed bounds.
"""
struct UniformMarginal{T1 <: Real} <: AbstractMarginal{T1}
	lb :: T1
	ub :: T1
end

"""
	$(SIGNATURES)

Normal(mean, std) marginal. Should store the Normal object +++
"""
struct NormalMarginal{T1 <: Real} <: AbstractMarginal{T1}
	mean :: T1
	std :: T1
end

"""
	$(SIGNATURES)

LogNormal(mean, std) marginal with offset `lb`. Should store the Normal object +++
"""
struct LogNormalMarginal{T1 <: Real} <: AbstractMarginal{T1}
	lb :: T1
	mean :: T1
	std :: T1
end

"""
	$(SIGNATURES)

Beta marginal distribution. Characterized by lower bound, upper bound, and Beta parameters alpha and beta.
"""
struct BetaMarginal{T1} <: AbstractMarginal{T1}
	lb :: T1
	ub :: T1
	betaDistr :: Distributions.Beta{T1}
end

"""
	$(SIGNATURES)

Bounded marginal with unkown distribution.
"""
struct BoundedMarginal{T1} <: AbstractMarginal{T1}
	lb :: T1
	ub :: T1
end

"""
	$(SIGNATURES)

Marginal with unknown distribution and no bounds.
"""
struct UnboundedMarginal{T1} <: AbstractMarginal{T1} end

"""
	$(SIGNATURES)

Fallback marginal for cases where no automatic checking of endowments takes place.
"""
struct UnknownMarginal{T1} <: AbstractMarginal{T1} end

## --------  Endowment info

"""
	$(SIGNATURES)

Meta information about one endowment. Contains a name, label, and its marginal distribution.

# Example
```
e = Endowment(:afqt, "Test score", PercentileMarginal{Float64}())
```
"""
struct Endowment{T1}
	name :: Symbol
	label :: String
	marginal :: AbstractMarginal{T1}
end

## ---------  EndowmentDraws

struct EndowmentDraws
	meta :: Vector{Endowment}
	draws :: DataFrame
end

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


text:
@testset "normal" begin

    for D in [1, 5]
        let rng = MersenneTwister(123456), N = 10_000_000, D = 5, tol = 1e-2

            # Test dimensionality.
            @test dim(DiagonalStandardNormal(1)) == 1
            @test dim(DiagonalStandardNormal(D)) == D

            # Allocating rng.
            @test abs(mean(rand(rng, DiagonalStandardNormal(1), N))) < tol
            @test abs(mean(rand(rng, DiagonalStandardNormal(D), N))) < tol
            @test abs(std(rand(rng, DiagonalStandardNormal(1), N)) - 1) < tol
            @test abs(std(rand(rng, DiagonalStandardNormal(D), N)) - 1) < tol

            # Test logpdf.
            let x = randn(rng, D)
                d = DiagonalStandardNormal(D)
                d′ = MvNormal(PDMats.ScalMat(D, 1.0))
                @test logpdf(d, x) ≈ logpdf(d′, x)
                for N′ in [1, 7]
                    X = randn(rng, D, N′)
                    @test logpdf(d, X) ≈ logpdf(d′, X)
                end
            end

            # Check that gradients propagate properly through second argument.
            for N′ in [1, 6]
                x = randn(rng, D)
                X = randn(rng, D, N′)
                d = DiagonalStandardNormal(D)
                @test ∇(x->logpdf(d, x))(x)[1] ≈ -x
                @test ∇(X->sum(logpdf(d, X)))(X)[1] ≈ -X
            end
        end
    end
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module RandTest
using ParallelAccelerator

#ParallelAccelerator.DomainIR.set_debug_level(4)
#ParallelAccelerator.ParallelIR.set_debug_level(4)
#ParallelAccelerator.CGen.set_debug_level(4)
#ParallelAccelerator.set_debug_level(4)

@acc function simple_rand(size1::Int64, size2::Int64)
    A = rand(size1,size2)
    C = A.*2.0
    return C
end

@acc function simple_randn(size1::Int64, size2::Int64)
    A = randn(size1,size2)
    C = A.*2.0
    return C
end

@acc function tuple_rand()
    y = [1.0 2.0 3.0 ;1.1 2.1 3.1]
    y += 0.1*rand(size(y))
    return y
end

@acc function tuple_randn()
    y = [1.0 2.0 3.0 ;1.1 2.1 3.1]
    y += 0.1*randn(size(y))
    return y
end

function test1()
    return simple_rand(2,3)
end

function test2()
    return simple_randn(2,3)
end

function test3()
    return tuple_rand()
end

function test4()
    return tuple_randn()
end

end

using Base.Test
println("Testing rand()...")
@test all(RandTest.test1() .<= ones(2,3).*2.0) && all(RandTest.test1() .>= zeros(2,3))
@test size(RandTest.test2())==(2,3)
@test size(RandTest.test3())==(2,3)
@test size(RandTest.test4())==(2,3)
println("Done testing rand()...")


text:
using Test;

function run(fname, version=1)
    f = open(fname, "r")
    lines = readlines(f)
    close(f)

    result = 0
    for l in lines
        if version == 1
            result += compute(l)
        elseif version == 2
            result += compute2(l)
        end
    end

    println("[Part ", version, "] Result = ", result)
    return result
end

function add(a, b)
    return a + b
end

function multiply(a, b)
    return a * b
end

funcs = Dict()
funcs["+"] = add
funcs["*"] = multiply

# WARNING: This is special math
function compute(s, version=1)
    if occursin('(', s)
        # if you find a pair of parenthesis, compute what's inside recursively
        start_idx = findlast(isequal('('), s)
        end_idx = findnext(isequal(')'), s, start_idx)

        pre = s[1:start_idx - 1]
        mid = compute(s[start_idx + 1:end_idx - 1], version)
        post = s[end_idx + 1:length(s)]
        return compute(pre * string(mid) * post, version)
    else
        parts = split(s)
        while length(parts) > 1
            start = 1
            if version == 2 && in("+", parts)
                # "Advanced Math": plus is resolved before multiply
                start = findfirst(isequal("+"), parts) - 1
            end
            val = resolve(parts[start:start + 2])
            parts = [parts[1:start - 1]..., string(val), parts[start + 3:length(parts)]...]
        end

        total = parse(Int, parts[1])
        return total
    end
end

function resolve(list)
    # println(list)
    # given a list of 3 items (value1, operation, value2), reduces to 1 item (output value)
    v1 = parse(Int, list[1])
    v2 = parse(Int, list[3])
    op = add
    if list[2] == "*"
        op = multiply
    end

    return op(v1, v2)
end

function compute2(s)
    return compute(s, 2)
end

# @test compute("1 + 2 * 3 + 4 * 5 + 6") == 71
# @test compute("2 * 3 + (4 * 5)") == 26
# @test compute("5 + (8 * 3 + 9 + 3 * 4 * 3)") == 437
# @test compute("5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))") == 12240
# @test compute("((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2") == 13632

# run("18.txt")

@test resolve(["1", "+", "3"]) == 4
@test resolve(["1", "*", "3"]) == 3

@test compute2("1 + 2 * 3 + 4 * 5 + 6") == 231
@test compute2("1 + (2 * 3) + (4 * (5 + 6))") == 51
@test compute2("2 * 3 + (4 * 5)") == 46
@test compute2("5 + (8 * 3 + 9 + 3 * 4 * 3)") == 1445
@test compute2("5 * 9 * (7 * 3 * 3 + 9 * 3 + (8 + 6 * 4))") == 669060
@test compute2("((2 + 4 * 9) * (6 + 9 * 8 + 6) + 6) + 2 + 4 * 2") == 23340

run("18.txt", 2)


text:
# This file is a part of ValueShapes.jl, licensed under the MIT License (MIT).

using ValueShapes
using Test

@testset "valueaccessor" begin
    acc = @inferred ValueAccessor(ArrayShape{Real}(2,3), 2)
    @test @inferred(valshape(acc)) == ArrayShape{Real,2}((2, 3))
    data = [1, 2, 3, 4, 5, 6, 7, 8, 9]
    @test @inferred(data[acc]) == [3 5 7; 4 6 8]
end


text:
module FStep

using Images
import JLD
using XGBoost
using MLBase
using Logging
@Logging.configure(level=DEBUG)

const debug = false

const data_dir = "../data"
const trainresize_dir = joinpath(data_dir, "trainResized")
const testresize_dir = joinpath(data_dir, "testResized")
const train_lb = joinpath(data_dir, "trainLabels.csv")
const sample_label_fl = joinpath(data_dir,"sample_label_dict.jld")
const label_sample_fl = joinpath(data_dir,"label_sample_dict.jld")

const train_fl = joinpath(data_dir, "train.jld")
const val_fl = joinpath(data_dir, "val.jld")
const test_fl = joinpath(data_dir, "test.jld")

include("preprocess.jl")
include("model.jl")
include("eval.jl")
include("dataflow.jl")
include("utils.jl")

end


text:
function _interpret(engine, wordstr)
    isnothing(wordstr) || wordstr == "" && return true
    word = codeof(engine, wordstr)
    if engine.mode == MODE_COMPILE && (!(word isa ExecutionToken) || !word.immediate)
        return _compile(engine, wordstr)
    else
        if isnothing(word)
            print(engine.out, "$(wordstr)?\n")
            return false
        end
        execute_codeword(engine, word) # No branching during interpretation
        return true
    end
end

# Execute the word either by dispatching to its Julia function if its a native word,
# or by recursively executing its compiled definition. (The "inner interpreter")
function execute_codeword(engine, word::ExecutionToken, parent = nothing, myidx = 0)
    code = word.code
    if code === op_swap # Fast, type.stabilized routes for a few words
        return code(engine, parent, myidx)
    elseif code === op_dup
        return code(engine, parent, myidx)
    elseif code === op_drop
        return code(engine, parent, myidx)
    elseif code isa Function
        return code(engine, parent, myidx)
    elseif code isa Number
        push!(engine.stack, code)
        return 1
    else
        idx = 1 # DOCOL
        codelength = length(code)
        while idx > 0 && idx <= codelength
            idx += execute_codeword(engine, code[idx], word, idx)
        end
        return 1
    end
end

function execute_codeword(engine, num::Number, parent = nothing, myidx = 0)
    push!(engine.stack, num)
    return 1
end

# read a space-separated word from io and return (word, end_of_line)
function word(io)
    buf = Char[]
    chr = ' '
    while isspace(chr) && !eof(io)
        chr = read(io, Char)
    end
    while true
        if isspace(chr)
            return (String(buf), chr == '\n' || chr == '\r' || eof(io))
        end
        push!(buf, chr)
        chr = eof(io) ? '\n' : read(io, Char)
    end
end

function interpret(machine, sentence)
    oldinput = machine.input
    machine.input = IOBuffer(sentence)
    try
        ok = true
        while !eof(machine.input) && ok
            (w, eol) = word(machine.input)
            ok = _interpret(machine, w)
        end
        return ok
    finally
        machine.input = oldinput
    end
end

function repl(engine = interpreter(); silent = false)
    if !silent
        print(engine.out, "ChainForth.jl v\"$VERSION\":\n")
    end
    while !eof(engine.input)
        try
            sentence = readline(stdin)
            if sentence == "exit"
                return
            end
            ok = interpret(engine, sentence)
            if ok && !silent
                println(engine.out, " ok")
            end
        catch e
            if e isa ArgumentError && e.msg == "array must be non-empty"
                println(engine.out, " Stack underflow")
            elseif e isa InterruptException
                return nothing
            else
                println(engine.out)
                if !silent
                    for (exc, bt) in Base.catch_stack()
                        showerror(engine.out, exc, bt)
                        println(engine.out)
                    end
                end
            end
        end
    end
end


text:
abstract type ALogicOperator <: AOperator end
abstract type AAnd <: ALogicOperator end
struct And <: AAnd
end
∧(x::Node, y::Node) = Node("∧", AAnd, (x, y))
export ∧

abstract type AForAll <: ALogicOperator end
ForAll(variable::Node, set::Node, statement::Node) = Node("∀", AForAll, (variable, set, statement))

abstract type AOr <: ALogicOperator end
struct Or <: AOr
end
∨(x::Node, y::Node) = Node("∨", AOr, (x, y))
export ∨

abstract type AExists <: ALogicOperator end
Exists(variable::Node, set::Node, statement::Node) = Node("∃", AOr, (variable, set, statement))
export Exists

abstract type ANot <: ALogicOperator end
struct Not <: ANot
end
Base.:(~)(x::Node) = Node("~", ANot, (x,))

abstract type AImplies <: ALogicOperator end
struct Implies <: AImplies
end
⟹(x::Node, y::Node) = Node("⟹", AImplies, (x, y))
export ⟹

abstract type AEquivalent <: ALogicOperator end
struct Equivalent <: AEquivalent
end
≡(x::Node, y::Node) = Node("≡", AEquivalent, (x, y))
export ≡

abstract type ASubset <: ALogicOperator end
struct Subset <: ASubset
end
⊂(x::Node, y::Node) = Node("⊂", ASubset, (x, y))
export ⊂

abstract type AElementOf <: ALogicOperator end
struct ElementOf <: AElementOf
end
∈(x::Node, y::Node) = Node("∈", AElementOf, (x, y))
export ∈


text:
using Turing, ContinuousBenchmarks, DelimitedFiles

odeDataRaw = readdlm(joinpath(ContinuousBenchmarks.DATA_DIR, "toy-data", "ode.csv"), ',')

t = Vector{Float64}(odeDataRaw[2:end, 1])
x0 = Vector{Float64}(odeDataRaw[2:end, 2])
x1 = Vector{Float64}(odeDataRaw[2:end, 3])

function sho(t, y, theta, x_r, x_i)
    dydt = Vector{Real}(undef, 2)
    dydt[1] = y[2]
    dydt[2] = -y[1] - theta[1] * y[2]

    return dydt
end


text:
schedulcolumn=6
jobname=Array{Any}(["スケジュール1","スケジュール2","スケジュール3","スケジュール4","スケジュール5","スケジュール6","スケジュール7","スケジュール8","スケジュール9","スケジュール10"])
runtime=Array{Any}([30,60,15,30,60,45,30,15,60,30])
Specifiedtime=Array{Any}([true,true,true,true,true,false,false,false,false,false])
JobStartTime=Array{Any}(["0:30","0:45","0:15","0:00","1:00",missing,missing,missing,missing,missing])
JobEndTime=Array{Any}(["1:00","1:45","0:30","0:30","2:00",missing,missing,missing,missing,missing])
#   0:00
col1=Array{Any}([missing,missing,missing,1,missing,missing,missing,missing,missing,missing])
#    0:15
col2=Array{Any}([missing,missing,1,1,missing,missing,missing,missing,missing,missing])
#    0:30
col3=Array{Any}([1,missing,missing,missing,missing,missing,missing,missing,missing,missing])
#    0:45
col4=Array{Any}([1,1,missing,missing,missing,missing,missing,missing,missing,missing])
#    1:00
col5=Array{Any}([missing,1,missing,missing,1,missing,missing,missing,missing,missing])
#    1:15
col6=Array{Any}([missing,1,missing,missing,1,missing,missing,missing,missing,missing])
#    1:30
col7=Array{Any}([missing,1,missing,missing,1,missing,missing,missing,missing,missing])
#    1:45
col8=Array{Any}([missing,missing,missing,missing,1,missing,missing,missing,missing,missing])
#    2:00
col9=Array{Any}([missing,missing,missing,missing,missing,missing,missing,missing,missing,missing])
#    DataFrame作成
schedulemaster=DataFrame(jobname=jobname,runtime=runtime,Specifiedtime=Specifiedtime,JobStartTime=JobStartTime,JobEndTime=JobEndTime,col1=col1,col2=col2,col3=col3,col4=col4,col5=col5,col6=col6,col7=col7,col8=col8,col9=col9)
rename!(schedulemaster,[:jobname,:runtime,:Specifiedtime,:JobStartTime,:JobEndTime,Symbol("0:00"),Symbol("0:15"),Symbol("0:30"),Symbol("0:45"),Symbol("1:00"),Symbol("1:15"),Symbol("1:30"),Symbol("1:45"),Symbol("2:00")])

InputFilePath=joinpath(@__DIR__, "schedule.xlsx")
OutputFilePath=joinpath(@__DIR__, "scheduleoutput.xlsx")

scheduleplan,robotn,run_unit_time,jobn,timen,schedule,blocktime_start,blocktime_end,blocktime_dow=readprerequisite(InputFilePath,"parameters","schedule")

rename!(schedulemaster,names(scheduleplan))


text:
using Pkg;
Pkg.activate(@__DIR__)
Pkg.instantiate()

using MLJ

using HTTP
using CSV
import DataFrames: DataFrame, select!, Not
req1 = HTTP.get("http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.data")
req2 = HTTP.get("http://archive.ics.uci.edu/ml/machine-learning-databases/horse-colic/horse-colic.test")
header = ["surgery", "age", "hospital_number",
    "rectal_temperature", "pulse",
    "respiratory_rate", "temperature_extremities",
    "peripheral_pulse", "mucous_membranes",
    "capillary_refill_time", "pain",
    "peristalsis", "abdominal_distension",
    "nasogastric_tube", "nasogastric_reflux",
    "nasogastric_reflux_ph", "feces", "abdomen",
    "packed_cell_volume", "total_protein",
    "abdomcentesis_appearance", "abdomcentesis_total_protein",
    "outcome", "surgical_lesion", "lesion_1", "lesion_2", "lesion_3",
    "cp_data"]
csv_opts = (header=header, delim=' ', missingstring="?",
            ignorerepeated=true)
data_train = CSV.read(req1.body; csv_opts...)
data_test  = CSV.read(req2.body; csv_opts...)
@show size(data_train)
@show size(data_test)

unwanted = [:lesion_1, :lesion_2, :lesion_3]
data = vcat(data_train, data_test)
select!(data, Not(unwanted));

train = 1:nrows(data_train)
test = last(train) .+ (1:nrows(data_test));

datac = coerce(data, autotype(data));

sch0 = schema(data)
sch = schema(datac)

old_scitype_given_name = Dict(
sch0.names[j] => sch0.scitypes[j] for j in eachindex(sch0.names))

length(unique(datac.hospital_number))

datac = select!(datac, Not(:hospital_number));

datac = coerce(datac, autotype(datac, rules=(:discrete_to_continuous,)));

missing_outcome = ismissing.(datac.outcome)
idx_missing_outcome = missing_outcome |> findall

train = setdiff!(train |> collect, idx_missing_outcome)
test = setdiff!(test |> collect, idx_missing_outcome)
datac = datac[.!missing_outcome, :];

for name in names(datac)
    col = datac[:, name]
    ratio_missing = sum(ismissing.(col)) / nrows(datac) * 100
    println(rpad(name, 30), round(ratio_missing, sigdigits=3))
end

unwanted = [:peripheral_pulse, :nasogastric_tube, :nasogastric_reflux,
        :nasogastric_reflux_ph, :feces, :abdomen, :abdomcentesis_appearance, :abdomcentesis_total_protein]
select!(datac, Not(unwanted));

@load FillImputer
filler = machine(FillImputer(), datac)
fit!(filler)
datac = transform(filler, datac)

cat_fields = filter(schema(datac).names) do field
    datac[:, field] isa CategoricalArray
end

for f in cat_fields
    datac[!, f] = get.(datac[:, f])
end

datac.pulse = coerce(datac.pulse, Count)
datac.respiratory_rate = coerce(datac.pulse, Count)

sch1 = schema(datac)

CSV.write("horse.csv", datac)


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

"""Mirror callback function

Function sets `+refs/*:refs/*` refspecs and `mirror` flag for remote reference.
"""
function mirror_callback(remote::Ptr{Ptr{Cvoid}}, repo_ptr::Ptr{Cvoid},
                         name::Cstring, url::Cstring, payload::Ptr{Cvoid})
    # Create the remote with a mirroring url
    fetch_spec = "+refs/*:refs/*"
    err = ccall((:git_remote_create_with_fetchspec, :libgit2), Cint,
                (Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Cstring, Cstring),
                remote, repo_ptr, name, url, fetch_spec)
    err != 0 && return Cint(err)

    # And set the configuration option to true for the push command
    config = GitConfig(GitRepo(repo_ptr,false))
    name_str = unsafe_string(name)
    err= try set!(config, "remote.$name_str.mirror", true)
         catch; -1
         finally close(config)
         end
    err != 0 && return Cint(err)
    return Cint(0)
end

"""
    LibGit2.is_passphrase_required(private_key) -> Bool

Return `true` if the `private_key` file requires a passphrase, `false` otherwise.
"""
function is_passphrase_required(private_key::AbstractString)
    !isfile(private_key) && return false

    # In encrypted private keys, the second line is "Proc-Type: 4,ENCRYPTED"
    return open(private_key) do f
        readline(f)
        readline(f) == "Proc-Type: 4,ENCRYPTED"
    end
end

function user_abort()
    # Note: Potentially it could be better to just throw a Julia error.
    ccall((:giterr_set_str, :libgit2), Cvoid,
          (Cint, Cstring), Cint(Error.Callback),
          "Aborting, user cancelled credential request.")
    return Cint(Error.EUSER)
end

function prompt_limit()
    ccall((:giterr_set_str, :libgit2), Cvoid,
          (Cint, Cstring), Cint(Error.Callback),
          "Aborting, maximum number of prompts reached.")
    return Cint(Error.EAUTH)
end

function exhausted_abort()
    ccall((:giterr_set_str, :libgit2), Cvoid,
          (Cint, Cstring), Cint(Error.Callback),
          "All authentication methods have failed.")
    return Cint(Error.EAUTH)
end

function authenticate_ssh(libgit2credptr::Ptr{Ptr{Cvoid}}, p::CredentialPayload, username_ptr)
    cred = p.credential::SSHCredential
    revised = false

    # Use a filled credential as-is on the first pass. Reset password on sucessive calls.
    if p.first_pass && isfilled(cred)
        revised = true
    elseif !p.first_pass
        cred.pass = ""
    end

    # first try ssh-agent if credentials support its usage
    if p.use_ssh_agent && username_ptr != Cstring(C_NULL) && (!revised || !isfilled(cred))
        err = ccall((:git_cred_ssh_key_from_agent, :libgit2), Cint,
                    (Ptr{Ptr{Cvoid}}, Cstring), libgit2credptr, username_ptr)

        p.use_ssh_agent = false  # use ssh-agent only one time
        err == 0 && return Cint(0)
    end

    if p.use_env && (!revised || !isfilled(cred))
        if isempty(cred.user) && username_ptr != Cstring(C_NULL)
            cred.user = unsafe_string(username_ptr)
        end

        cred.prvkey = Base.get(ENV, "SSH_KEY_PATH") do
            default = joinpath(homedir(), ".ssh", "id_rsa")
            if isempty(cred.prvkey) && isfile(default)
                default
            else
                cred.prvkey
            end
        end

        cred.pubkey = Base.get(ENV, "SSH_PUB_KEY_PATH") do
            default = cred.prvkey * ".pub"
            if isempty(cred.pubkey) && isfile(default)
                default
            else
                cred.pubkey
            end
        end

        cred.pass = Base.get(ENV, "SSH_KEY_PASS", cred.pass)

        revised = true
        p.use_env = false
    end

    if p.remaining_prompts > 0 && (!revised || !isfilled(cred))
        if isempty(cred.user) || username_ptr == Cstring(C_NULL)
            url = git_url(scheme=p.scheme, host=p.host)
            response = Base.prompt("Username for '$url'", default=cred.user)
            response === nothing && return user_abort()
            cred.user = response
        end

        url = git_url(scheme=p.scheme, host=p.host, username=cred.user)

        # For SSH we need a private key location
        last_private_key = cred.prvkey
        if !isfile(cred.prvkey) || !revised || !haskey(ENV, "SSH_KEY_PATH")
            response = Base.prompt("Private key location for '$url'", default=cred.prvkey)
            response === nothing && return user_abort()
            cred.prvkey = expanduser(response)

            # Only update the public key if the private key changed
            if cred.prvkey != last_private_key
                cred.pubkey = cred.prvkey * ".pub"
            end
        end

        # For SSH we need a public key location. Avoid asking about the public key as
        # typically this will just annoy users.
        stale = !p.first_pass && cred.prvkey == last_private_key && cred.pubkey != cred.prvkey * ".pub"
        if isfile(cred.prvkey) && (stale || !isfile(cred.pubkey))
            response = Base.prompt("Public key location for '$url'", default=cred.pubkey)
            response === nothing && return user_abort()
            cred.pubkey = expanduser(response)
        end

        # Ask for a passphrase when the private key exists and requires a passphrase
        if isempty(cred.pass) && is_passphrase_required(cred.prvkey)
            if Sys.iswindows()
                response = Base.winprompt(
                    "Your SSH Key requires a password, please enter it now:",
                    "Passphrase required", cred.prvkey; prompt_username=false)
                response === nothing && return user_abort()
                cred.pass = response[2]
            else
                response = Base.prompt("Passphrase for $(cred.prvkey)", password=true)
                response === nothing && return user_abort()
                cred.pass = response
                isempty(cred.pass) && return user_abort()  # Ambiguous if EOF or newline
            end
        end

        revised = true

        p.remaining_prompts -= 1
        p.remaining_prompts <= 0 && return prompt_limit()
    end

    if !revised
        return exhausted_abort()
    end

    return ccall((:git_cred_ssh_key_new, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Cstring, Cstring, Cstring, Cstring),
                 libgit2credptr, cred.user, cred.pubkey, cred.prvkey, cred.pass)
end

function authenticate_userpass(libgit2credptr::Ptr{Ptr{Cvoid}}, p::CredentialPayload)
    cred = p.credential::UserPasswordCredential
    revised = false

    # Use a filled credential as-is on the first pass. Reset password on sucessive calls.
    if p.first_pass && isfilled(cred)
        revised = true
    elseif !p.first_pass
        cred.pass = ""
    end

    if p.use_git_helpers && (!revised || !isfilled(cred))
        git_cred = GitCredential(p.config, p.url)

        # Use `deepcopy` to ensure zeroing the `git_cred` doesn't also zero the `cred`s copy
        cred.user = deepcopy(coalesce(git_cred.username, ""))
        cred.pass = deepcopy(coalesce(git_cred.password, ""))
        securezero!(git_cred)
        revised = true

        p.use_git_helpers = false
    end

    if p.remaining_prompts > 0 && (!revised || !isfilled(cred))
        url = git_url(scheme=p.scheme, host=p.host)
        username = isempty(cred.user) ? p.username : cred.user
        if Sys.iswindows()
            response = Base.winprompt(
                "Please enter your credentials for '$url'", "Credentials required",
                username; prompt_username=true)
            response === nothing && return user_abort()
            cred.user, cred.pass = response
        else
            response = Base.prompt("Username for '$url'", default=username)
            response === nothing && return user_abort()
            cred.user = response

            url = git_url(scheme=p.scheme, host=p.host, username=cred.user)
            response = Base.prompt("Password for '$url'", password=true)
            response === nothing && return user_abort()
            cred.pass = response
            isempty(cred.pass) && return user_abort()  # Ambiguous if EOF or newline
        end

        revised = true

        p.remaining_prompts -= 1
        p.remaining_prompts <= 0 && return prompt_limit()
    end

    if !revised
        return exhausted_abort()
    end

    return ccall((:git_cred_userpass_plaintext_new, :libgit2), Cint,
                 (Ptr{Ptr{Cvoid}}, Cstring, Cstring),
                 libgit2credptr, cred.user, cred.pass)
end

"""
    credential_callback(...) -> Cint

A LibGit2 credential callback function which provides different credential acquisition
functionality w.r.t. a connection protocol. The `payload_ptr` is required to contain a
`LibGit2.CredentialPayload` object which will keep track of state and settings.

The `allowed_types` contains a bitmask of `LibGit2.Consts.GIT_CREDTYPE` values specifying
which authentication methods should be attempted.

Credential authentication is done in the following order (if supported):
- SSH agent
- SSH private/public key pair
- Username/password plain text

If a user is presented with a credential prompt they can abort the prompt by typing `^D`
(pressing the control key together with the `d` key).

**Note**: Due to the specifics of the `libgit2` authentication procedure, when
authentication fails, this function is called again without any indication whether
authentication was successful or not. To avoid an infinite loop from repeatedly
using the same faulty credentials, we will keep track of state using the payload.

For addition details see the LibGit2 guide on
[authenticating against a server](https://libgit2.github.com/docs/guides/authentication/).
"""
function credentials_callback(libgit2credptr::Ptr{Ptr{Cvoid}}, url_ptr::Cstring,
                              username_ptr::Cstring,
                              allowed_types::Cuint, payload_ptr::Ptr{Cvoid})
    err = Cint(0)

    # get `CredentialPayload` object from payload pointer
    @assert payload_ptr != C_NULL
    p = unsafe_pointer_to_objref(payload_ptr)::CredentialPayload

    # Parse URL only during the first call to this function. Future calls will use the
    # information cached inside the payload.
    if isempty(p.url)
        p.url = unsafe_string(url_ptr)
        m = match(URL_REGEX, p.url)

        p.scheme = coalesce(m[:scheme], "")
        p.username = coalesce(m[:user], "")
        p.host = m[:host]

        # When an explicit credential is supplied we will make sure to use the given
        # credential during the first callback by modifying the allowed types. The
        # modification only is in effect for the first callback since `allowed_types` cannot
        # be mutated.
        if p.explicit !== nothing
            cred = p.explicit

            # Copy explicit credentials to avoid mutating approved credentials.
            p.credential = deepcopy(cred)

            if isa(cred, SSHCredential)
                allowed_types &= Cuint(Consts.CREDTYPE_SSH_KEY)
            elseif isa(cred, UserPasswordCredential)
                allowed_types &= Cuint(Consts.CREDTYPE_USERPASS_PLAINTEXT)
            else
                allowed_types &= Cuint(0)  # Unhandled credential type
            end
        elseif p.cache !== nothing
            cred_id = credential_identifier(p.scheme, p.host)

            # Perform a deepcopy as we do not want to mutate approved cached credentials
            if haskey(p.cache, cred_id)
                p.credential = deepcopy(p.cache[cred_id])
            end
        end

        p.first_pass = true
    else
        p.first_pass = false
    end

    # use ssh key or ssh-agent
    if isset(allowed_types, Cuint(Consts.CREDTYPE_SSH_KEY))
        if p.credential === nothing || !isa(p.credential, SSHCredential)
            p.credential = SSHCredential(p.username)
        end
        err = authenticate_ssh(libgit2credptr, p, username_ptr)
        err == 0 && return err
    end

    if isset(allowed_types, Cuint(Consts.CREDTYPE_USERPASS_PLAINTEXT))
        if p.credential === nothing || !isa(p.credential, UserPasswordCredential)
            p.credential = UserPasswordCredential(p.username)
        end
        err = authenticate_userpass(libgit2credptr, p)
        err == 0 && return err
    end

    # No authentication method we support succeeded. The most likely cause is
    # that explicit credentials were passed in, but said credentials are incompatible
    # with the requested authentication method.
    if err == 0
        if p.explicit !== nothing
            ccall((:giterr_set_str, :libgit2), Cvoid, (Cint, Cstring), Cint(Error.Callback),
                  "The explicitly provided credential is incompatible with the requested " *
                  "authentication methods.")
        end
        err = Cint(Error.EAUTH)
    end
    return err
end

function fetchhead_foreach_callback(ref_name::Cstring, remote_url::Cstring,
                        oid_ptr::Ptr{GitHash}, is_merge::Cuint, payload::Ptr{Cvoid})
    fhead_vec = unsafe_pointer_to_objref(payload)::Vector{FetchHead}
    push!(fhead_vec, FetchHead(unsafe_string(ref_name), unsafe_string(remote_url),
        unsafe_load(oid_ptr), is_merge == 1))
    return Cint(0)
end

"C function pointer for `mirror_callback`"
mirror_cb() = cfunction(mirror_callback, Cint, Tuple{Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Cstring, Cstring, Ptr{Cvoid}})
"C function pointer for `credentials_callback`"
credentials_cb() = cfunction(credentials_callback, Cint, Tuple{Ptr{Ptr{Cvoid}}, Cstring, Cstring, Cuint, Ptr{Cvoid}})
"C function pointer for `fetchhead_foreach_callback`"
fetchhead_foreach_cb() = cfunction(fetchhead_foreach_callback, Cint, Tuple{Cstring, Cstring, Ptr{GitHash}, Cuint, Ptr{Cvoid}})

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: YingboMa/julia

max_forks_repo_name: YingboMa/julia


text:
# Trajectory-based speech parameter mapping for voice conversion
# based on the maximum likelihood criterion.
type TrajectoryGMMMap <: TrajectoryConverter
    gmmmap::GMMMap

    W::SparseMatrixCSC{Float64, Int}     # weight matrix
    Eʸ::Vector{Float64}                  # vectorized version of eq. (40)
    Dʸ::Array{Float64, 3}                # diagonal components of eq. (41)
    Dʸ⁻¹::SparseMatrixCSC{Float64, Int}  # Dʸ^-1 in eq. (41)

    function TrajectoryGMMMap(g::GMMMap, T::Int)
        D = dim(g)>>1 # dimension of static feature vector

        W = constructW(D, T)

        # the number of mixtures
        M = ncomponents(g)

        Σˣʸ = g.params.Σˣʸ
        Σʸʸ = g.params.Σʸʸ
        ΣʸˣΣˣˣ⁻¹ = g.params.ΣʸˣΣˣˣ⁻¹

        # pre-computations
        Dʸ = Array(Float64, 2D, 2D, M)
        for m=1:M
            Dʸ[:,:,m] = Σʸʸ[:,:,m] - ΣʸˣΣˣˣ⁻¹[:,:,m] * Σˣʸ[:,:,m]
            Dʸ[:,:,m] = Dʸ[:,:,m]^-1
        end

        new(g, W, zeros(0), Dʸ, spzeros(0, 0))
    end
end

Base.length(t::TrajectoryGMMMap) = div(size(t.W, 2), div(dim(t), 2))
dim(t::TrajectoryGMMMap) = dim(t.gmmmap)
ncomponents(t::TrajectoryGMMMap) = ncomponents(t.gmmmap)
Base.size(g::TrajectoryGMMMap) = (dim(g), length(g))

function compute_wt(t::Int, D::Int, T::Int)
    @assert t > 0
    w⁰ = spzeros(D, D*T)
    w¹ = spzeros(D, D*T)
    w⁰[:, (t-1)*D+1:t*D] = speye(D)

    if t >= 2
        w¹[:, (t-2)*D+1:(t-1)*D] = -0.5*speye(D)
    end
    if t < T
        w¹[:, t*D+1:(t+1)*D] = 0.5*speye(D)
    end

    [w⁰; w¹]
end

function constructW(D::Int, T::Int)
    W = spzeros(2D*T, D*T)
    for t=1:T
        W[2D*(t-1)+1:2D*t,:] = compute_wt(t, D, T)
    end
    W
end

# Mapping source spectral feature x to target spectral feature y
# so that maximize the likelihood of y given x.
function fvconvert(tgmm::TrajectoryGMMMap, X::Matrix{Float64})
    # input feature vector must contain delta feature
    D, T = size(X, 1)>>1, size(X, 2)
    2D == dim(tgmm) || throw(DimensionMismatch("Inconsistent dimentions."))

    if T != length(tgmm)
        tgmm.W = constructW(D, T)
    end

    # aliases
    g = tgmm.gmmmap
    μʸ = g.params.μʸ
    μˣ = g.params.μˣ
    ΣʸˣΣˣˣ⁻¹ = g.params.ΣʸˣΣˣˣ⁻¹
    W = tgmm.W

    # A suboptimum mixture sequence  eq. (37)
    m̂ = predict(g.px, X)

    # Compute Eʸ eq.(40)
    Eʸ = Array(Float64, 2D, T)
    for t=1:T
        m = m̂[t]::Int
        Eʸ[:,t] = μʸ[:,m] + ΣʸˣΣˣˣ⁻¹[:,:,m] * (X[:,t] - μˣ[:,m])
    end
    Eʸ = vec(Eʸ)
    tgmm.Eʸ = Eʸ # keep Eʸ for GV optimization
    @assert size(Eʸ) == (2D*T,)

    # Compute D^-1 eq.(41)
    Dʸ⁻¹ = blkdiag([sparse(tgmm.Dʸ[:,:,m̂[t]]) for t=1:T]...)
    tgmm.Dʸ⁻¹ = Dʸ⁻¹ # Keep Dʸ⁻¹ for GV

    @assert size(Dʸ⁻¹) == (2D*T, 2D*T)
    @assert issparse(Dʸ⁻¹)

    # Compute target static feature vector
    # eq. (39)
    WᵀDʸ⁻¹ = W' * Dʸ⁻¹
    @assert issparse(WᵀDʸ⁻¹)
    y = (WᵀDʸ⁻¹ * W) \ (WᵀDʸ⁻¹ * Eʸ)
    @assert size(y) == (D*T,)

    # Finally we get static feature vector
    reshape(y, D, T)
end

# Trajectory-based speech parameter mapping considering global variance
# based on the maximum likelihood criterion.
type TrajectoryGVGMMMap <: TrajectoryConverter
    tgmm::TrajectoryGMMMap
    μᵛ::Vector{Float64}
    Σᵛᵛ::Matrix{Float64}
    pᵥ::Matrix{Float64}

    function TrajectoryGVGMMMap(tgmm::TrajectoryGMMMap,
                                μᵛ::Vector{Float64}, # mean of GV
                                Σᵛᵛ::Matrix{Float64} # covars of GV
        )
        @assert sum(μᵛ .< 0) == 0
        new(tgmm, μᵛ, Σᵛᵛ, inv(Σᵛᵛ))
    end
end

Base.length(t::TrajectoryGVGMMMap) = Base.length(t.tgmm)
dim(t::TrajectoryGVGMMMap) = dim(t.tgmm)
ncomponents(t::TrajectoryGVGMMMap) = ncomponents(t.tgmm)
Base.size(g::TrajectoryGVGMMMap) = Base.size(t.tgmm)

# Mapping source spectral feature x to target spectral feature y
# so that maximize the likelihood of y given x with considering
# global variance.
# Note that step size `α` should be carefully chosen.
function fvconvert(tgv::TrajectoryGVGMMMap, X::Matrix{Float64};
                   epochs::Int=100,
                   α::Float64=1.0e-5, # step-size
                   verbose::Bool=false
                   )
    # Initialize target static features without considering GV
    y⁰ = fvconvert(tgv.tgmm, X)
    D, T = size(y⁰)

    # Better initial value based on eq. (58)
    y⁰ = sqrt(tgv.μᵛ ./ var(y⁰[1:D,:], 2)) .* (y⁰ .- mean(y⁰, 2)) .+ mean(y⁰, 2)

    ω = 1.0/(2T)

    # aliases
    Eʸ = tgv.tgmm.Eʸ
    W = tgv.tgmm.W
    Dʸ⁻¹ = tgv.tgmm.Dʸ⁻¹
    WᵀDʸ⁻¹ = W' * Dʸ⁻¹

    # update y based on gradient decent
    yⁱ = y⁰
    for epoch=1:epochs
        Δyⁱ = ω*(-WᵀDʸ⁻¹ * W * vec(yⁱ) + WᵀDʸ⁻¹ * Eʸ) + vec(gvgrad(tgv, yⁱ))
        if verbose
            println("Epoch #$(epoch): norm $(norm(Δyⁱ))")
        end
        @assert !any(isnan(Δyⁱ))
        Δyⁱ = reshape(Δyⁱ, D, T)
        # Eq. (52)
        yⁱ = yⁱ + α * Δyⁱ
    end

    yⁱ
end

# gvgrad computes gradient of the likelihood with regard to GV.
function gvgrad(tgv::TrajectoryGVGMMMap, y::Matrix{Float64})
    D, T = size(y)

    gv = var(y, 2) # global variance over time
    @assert size(gv) == (D, 1)
    μʸ = mean(y, 2)
    @assert size(μʸ) == (D, 1)

    v = Array(Float64, D, T)
    for t=1:T
        @inbounds v[:,t] = -2.0/T*(tgv.pᵥ'*(gv - tgv.μᵛ)) .* (y[:,t] - μʸ)
    end

    v
end


text:
"""
`module  JAC.Semiempirical`  ... a submodel of JAC that contains all methods to set-up and process (simple)  semiempirical estimations of
                                 atomic properties that cannot be calculated so easily in a rigorious manner; it is using JAC,
                                 JAC.ManyElectron, JAC.ImpactIonization.
"""
module Semiempirical

    using JAC, JAC.ManyElectron, JAC.ImpactIonization

    """
    `@enum   AtomicCrossSection`  ... defines a enumeration for (supported) cross section estimation.
    """
    @enum   AtomicCrossSection    NoCrossSection    ImpactIonizationCS

    """
    `JAC.AtomicCrossSection(sa::String)`  ... constructor for a given String.
    """
    function AtomicCrossSection(sa::String)
        if       sa == "No cross section"      wa = NoCrossSection
        elseif   sa == "impact-ionization"     wa = ImpactIonizationCS
        else
            error("stop a")
        end

        return( wa )
    end

    # `Base.show(io::IO, cs::AtomicCrossSection)`  ... prepares a proper printout of the variable cs::AtomicCrossSection.
    function Base.show(io::IO, cs::AtomicCrossSection)
        print(io, string(cs) )
    end

    # `Base.string(cs::AtomicCrossSection)`  ... provides a proper printout of the variable cs::AtomicCrossSection.
    function Base.string(cs::AtomicCrossSection)
        if       cs == NoCrossSection       return("No cross section")
        elseif   cs == ImpactIonizationCS   return("impact-ionization cross section")
        else     error("stop a")
        end
    end

    """
    `struct  Estimation`  ... defines a type for defining  (a framework for) semiempirical estimates of different properties that cannot
                              be calculated rigoriously in JAC.

        + calcCrossSections    ::Bool                        ... True, if atomic cross sections of some kind are to be calculated
        + crossSection         ::AtomicCrossSection          ... Selected cross sections type.
        + config               ::Configuration               ... A non-relativistic configurations for which computations are to be carried out.
        + crossSectionSettings ::Union{JAC.ImpactIonization.Settings} ... Provides the settings for the selected cross section estimates.
    """
    struct  Estimation
        calcCrossSections      ::Bool
        crossSection           ::AtomicCrossSection
        config                 ::Configuration
        crossSectionSettings   ::Union{JAC.ImpactIonization.Settings}
    end

    """
    `JAC.Semiempirical.Estimation()`  ... constructor for an 'empty' instance::Semiempirical.Estimation.
    """
    function Estimation()
        Estimation(false, NoCrossSection, Configuration(), JAC.ImpactIonization.Settings() )
    end

    # `Base.show(io::IO, estimation::Semiempirical.Estimation)`  ... prepares a proper printout of the variable estimation::Semiempirical.Estimation.
    function Base.show(io::IO, estimation::Semiempirical.Estimation)
        println(io, "calcCrossSections:        $(estimation.calcCrossSections)  ")
        println(io, "crossSection:             $(estimation.crossSection)  ")
        println(io, "config:                   $(estimation.config)  ")
        println(io, "crossSectionSettings:     $(estimation.crossSectionSettings)  ")
    end

    """
    `JAC.Semiempirical.Estimation("interactive")`  ... constructor to generate a new instance of Semiempirical.Estimation interactively by
         replying to some detailed dialog.  **Not yet fully implemented !**
    """
    function Semiempirical.Estimation(sa::String)
        sa != "interactive"   &&    error("Unsupported keystring = $sa.")
    end

end # module


text:
import Base.show_expr_type, Base.show_linenumber, Base.show_call
# Re-define the function. Otherwise, the show of an expression may be too verbose
function show_expr_type(io::IO, ty)
    return
end

function show_linenumber(io::IO, file, line)
    return
end

max_stars_repo_name: IntelLabs/Sparso

max_issues_repo_name: IntelLabs/Sparso

max_forks_repo_name: IntelLabs/Sparso


text:
# lists table-specific iterators
iters(session::AccumuloSession, tablename::AbstractString) = listIterators(client(session), handle(session), utf8(tablename))

# get a single iterator with its settings
iter(session::AccumuloSession, tablename::AbstractString, itername::AbstractString, scope::Integer) = getIteratorSetting(client(session), handle(session), utf8(tablename), utf8(itername), Int32(scope))

# create an iterator type
function iter(name::AbstractString, class::AbstractString, priority::Integer, properties::Dict=Dict())
    iter_properties = Dict{UTF8String,UTF8String}([utf8(string(n))=>utf8(string(v)) for (n,v) in properties]...)

    thriftbuild(IteratorSetting, Dict(:priority => Int32(priority),
        :name => utf8(name),
        :iteratorClass => utf8(class),
        :properties => iter_properties))
end

# attach an iterator to a table
function add_iter(session::AccumuloSession, tablename::AbstractString, iter_setting::IteratorSetting, scopes::SET; check::Bool=true)
    iter_scopes = Set([Int32(scope) for scope in scopes])

    if check
        checkIteratorConflicts(client(session), handle(session), utf8(tablename), iter_setting, iter_scopes)
    end
    attachIterator(client(session), handle(session), utf8(tablename), iter_setting, iter_scopes)
end

# detach the iterator from the table
remove_iter(session::AccumuloSession, tablename::AbstractString, iter_setting::IteratorSetting, scopes::SET) = remove_iter(session, tablename, iter_setting.name, scopes)
function remove_iter(session::AccumuloSession, tablename::AbstractString, name::AbstractString, scopes::SET)
    iter_scopes = Set([Int32(scope) for scope in scopes])
    removeIterator(client(session), handle(session), utf8(tablename), iter_scopes)
end

# check if there is any conflicts with the iterator name and priority
function check_iter(session::AccumuloSession, tablename::AbstractString, iter_setting::IteratorSetting, scopes::SET)
    iter_scopes = Set([Int32(scope) for scope in scopes])
    checkIteratorConflicts(client(session), handle(session), utf8(tablename), iter_setting, iter_scopes)
end