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. With the macro `@law` we are able to define a function whose arguments depend on the coordinates, without the need of passing such coordinates as an argument.
using LinearAlgebra: tr
@law σ_s(ε) = λ_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.
@law σ_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 discrete model of each of the sub-domains, which are constructed restricting the global discrete model to the corresponding part. This is done by calling the `DiscreteModel` function with the desired geometrical part label, i.e. `"solid"` and `"fluid"`, respectively.
model_solid = DiscreteModel(model,"solid")
model_fluid = DiscreteModel(model,"fluid")

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

# 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 a 2-dimensional `VectorValue` type, a `:Lagrangian` reference FE element of order `k` 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(
    model=model_fluid,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=k,
    conformity =:H1,
    dirichlet_tags=["inlet", "noSlip", "cylinder"])

Vs = TestFESpace(
    model=model_solid,
    valuetype=VectorValue{2,Float64},
    reffe=:Lagrangian,
    order=k,
    conformity =:H1,
    dirichlet_tags=["fixed"])

# For the pressure test FE space, we use the fluid discrete model, a scalar value type, a `:Lagrangian` reference FE element of order `k-1` and `:C0` conformity.
Qf = TestFESpace(
  model=model_fluid,
  valuetype=Float64,
  order=k-1,
  reffe=:Lagrangian,
  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="weakForm"></a>
# ```
# ### Weak form
# We now introduce the solution and test function vectors as $\mathbf{x}^h=[\mathbf{u}^h_{\rm S},\mathbf{u}^h_{\rm F}, p^h_{\rm F}]^T$ and $\mathbf{y}^h=[\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{x}^h \in X$ such that
# ```math
# a_s(\mathbf{x}^h,\mathbf{y}^h) + a_f(\mathbf{x}^h,\mathbf{y}^h) + a_{fs}(\mathbf{x}^h,\mathbf{y}^h) = l_s(\mathbf{y}^h) + l_f(\mathbf{y}^h) + l_{f,\Gamma_N}(\mathbf{y}^h)\qquad\forall\mathbf{y}^h\in Y,
# ```
# with the following definitions:

# - $a_s(\mathbf{x}^h,\mathbf{y}^h)$ is the bilinear form associated with the solid counterpart, defined as
# ```math
# a_s(\mathbf{x}^h,\mathbf{y}^h)=(\varepsilon(\mathbf{v}^h_{\rm S}),\boldsymbol{\sigma}_s(\mathbf{u}^h_{\rm S})).
# ```
function a_s(x,y)
  us,uf,p = x
  vs,vf,q = y
	ε(vs) ⊙ σ_s(ε(us))
end

# - $a_f(\mathbf{x}^h,\mathbf{y}^h)$ is the bilinear form associated with the fluid counterpart, defined as
# ```math
# a_f(\mathbf{x}^h,\mathbf{y}^h)=(\varepsilon(\mathbf{v}^h_{\rm F}),\boldsymbol{\sigma}^{dev}_f(\mathbf{u}^h_{\rm F})) - (\nabla\cdot\mathbf{v}^h_{\rm F},p_{\rm F}) + (q_{\rm F}, \nabla\cdot\mathbf{u}^h_{\rm F}).
# ```
function a_f(x,y)
  us,uf,p = x
  vs,vf,q = y
  (ε(vf) ⊙ σ_dev_f(ε(uf))) - (∇⋅vf)*p + q*(∇⋅uf)
end

# - $a_{fs}(\mathbf{x}^h,\mathbf{y}^h)$ is the bilinear form associated with the coupling between fluid and solid counterparts. To difine this form we use 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]. The final expression for this term reads:
# ```math
# \begin{aligned}
# a_{fs}(\mathbf{x}^h,\mathbf{y}^h)=&\langle\gamma\frac{\mu_f}{h}(\mathbf{v}^h_{\rm F}-\mathbf{v}^h_{\rm S}),(\mathbf{u}^h_{\rm F}-\mathbf{u}^h_{\rm S})\rangle_{\Gamma_{\rm FS}}\\
# &-  \langle(\mathbf{v}^h_{\rm F}-\mathbf{v}^h_{\rm S}),\mathbf{n}_{\rm FS}\cdot\boldsymbol{\sigma}^{dev}_f(\mathbf{u}^h_{\rm F})\rangle_{\Gamma_{\rm FS}}
# +  \langle(\mathbf{v}^h_{\rm F}-\mathbf{v}^h_{\rm S}),p^h_{\rm F}\mathbf{n}_{\rm FS}\rangle_{\Gamma_{\rm FS}}\\
# &-  \chi\langle\mathbf{n}_{\rm FS}\cdot\boldsymbol{\sigma}^{dev}_f(\mathbf{v}^h_{\rm F}),(\mathbf{u}^h_{\rm F}-\mathbf{u}^h_{\rm S})\rangle_{\Gamma_{\rm FS}}
# +  \langle q^h_{\rm F}\mathbf{n}_{\rm FS},(\mathbf{u}^h_{\rm F}-\mathbf{u}^h_{\rm S})\rangle_{\Gamma_{\rm FS}}\\
# \end{aligned}
# ```
# 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.
const γ = 1.0
const χ = -1.0
function a_fs(x,y)
  us_Γ, uf_Γ, p_Γ = x
  vs_Γ, vf_Γ, q_Γ = y
  us, uf, p = us_Γ.⁻, uf_Γ.⁺, p_Γ.⁺
  vs, vf, q = vs_Γ.⁻, vf_Γ.⁺, q_Γ.⁺
  α*(vf-vs)⋅(uf-us) + (vf-vs)⋅(p*n_Γfs-n_Γfs⋅σ_dev_f(ε(uf))) + (q*n_Γfs-χ*n_Γfs⋅σ_dev_f(ε(vf)))⋅(uf-us)
end

# - $l_s(\mathbf{y}^h)$ is the linear form associated with the solid counterpart, defined as
# ```math
# l_s(\mathbf{y}^h)=(\mathbf{v}^h_{\rm S},s).
# ```
function l_s(y)
  vs,vf,q = y
  vs⋅s
end

# - $l_f(\mathbf{y}^h)$ is the linear form associated with the fluid counterpart, defined as
# ```math
# l_f(\mathbf{y}^h)=(\mathbf{v}^h_{\rm F},f).
# ```
function l_f(y)
  vs,vf,q = y
  vf⋅f + q*g
end

# - $l_{f,\Gamma_N}(\mathbf{y}^h)$ is the linear form associated with the fluid Neumann boundary condition, defined as
# ```math
# l_{f,\Gamma_N}(\mathbf{y}^h)=\langle\mathbf{v}^h_{\rm F},h\rangle_{\Gamma_N}.
# ```
function l_f_Γn(y)
  vs,vf,q = y
  vf⋅hN
end

# ```@raw HTML
# <a name="integration"></a>
# ```
# ### Numerical integration
# To define the quadrature rules used in the numerical integration of the different terms, we first need to generate the domain triangulation. Here we create the triangulation of the global domain, $\mathcal{T}$.
trian = Triangulation(model)

# The solid and fluid triangulations, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$, are constructed from the discrete models restricted to the respective parts.
trian_solid = Triangulation(model_solid)
trian_fluid = Triangulation(model_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.
trian_Γout = BoundaryTriangulation(model,"outlet")
n_Γout = get_normal_vector(trian_Γ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 models. Note that the normal field will point outwards with respect to the first entry.
trian_Γfs = InterfaceTriangulation(model_fluid,model_solid)
n_Γfs = get_normal_vector(trian_Γfs)

# 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.Arrays
using LinearAlgebra: norm
xe = get_cell_coordinates(trian_Γfs)
he = apply(x->norm(x[2]-x[1]),xe)
α = apply(h->γ*μ_f/h, he)

# Once we have all the triangulations, we can generate the quadrature rules to be applied each domain.
degree = 2*k
quad_solid = CellQuadrature(trian_solid,degree)
quad_fluid = CellQuadrature(trian_fluid,degree)

# Idem for the boundary part.
bdegree = 2*k
quad_Γout = CellQuadrature(trian_Γout,bdegree)

# Idem for the interface part.
idegree = 2*k
quad_Γfs = CellQuadrature(trian_Γfs,idegree)

# ```@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 first define the diferent discrete terms using the functions `AffineFETerm`, that assemble contributions on the left-hand side and the right-hand side of the system, for the fluid and solid interior terms evaluated in the respective triangulations, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$. The fluid Neumann boundary condition is assembled using the function `FESource`, which only affects the right-hand side of the system, evaluating the corresponding form on the triangulation of the outlet boundary, $\Gamma_{\rm F,N_{out}}$. The coupling terms are evaluated at the interface $\Gamma_{\rm FS}$ using the function `LinearFETerm`, assembling only to the left-hand side of the system.
t_Ω_s= AffineFETerm(a_s,l_s,trian_solid,quad_solid)
t_Ω_f = AffineFETerm(a_f,l_f,trian_fluid,quad_fluid)
t_Γfs = LinearFETerm(a_fs,trian_Γfs,quad_Γfs)
t_Γn_f = FESource(l_f_Γn,trian_Γout,quad_Γout)

# The final FE operator is constructed using the function `AffineFEOperator` and passing as arguments the trial and test FE spaces, $X$ and $Y$, and all the different terms previously defined.
op = AffineFEOperator(X,Y,t_Ω_s,t_Ω_f,t_Γn_f,t_Γfs)

# 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,"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.
uhs_solid = restrict(uhs, trian_solid)
uhf_fluid = restrict(uhf, trian_fluid)
ph_fluid = restrict(ph, trian_fluid)
writevtk(trian_solid,"trian_solid",cellfields=["uhs"=>uhs_solid])
writevtk(trian_fluid,"trian_fluid",cellfields=["ph"=>ph_fluid,"uhf"=>uhf_fluid])
# ![](../assets/fsi/Local_solution.png)

# ```@raw HTML
# <a name="QOIs"></a>
# ```
# ### Quantities of Interest
trian_ΓS = BoundaryTriangulation(model,["cylinder","interface"])
quad_ΓS = CellQuadrature(trian_ΓS,bdegree)
n_ΓS = get_normal_vector(trian_ΓS)
uh_ΓS = restrict(uhf_fluid,trian_ΓS)
ph_ΓS = restrict(ph_fluid,trian_ΓS)
FD, FL = sum( integrate( (n_ΓS⋅σ_dev_f(ε(uh_ΓS))) - ph_ΓS*n_ΓS, trian_ΓS, quad_ΓS ) )
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.

max_stars_repo_name: andre-massing/Tutorials

max_issues_repo_name: andre-massing/Tutorials

max_forks_repo_name: andre-massing/Tutorials


text:
"""
    actual_time_step(max_time_step, time_interval, time_step_options)

Actual size and number of time steps that should be taken if the maximum
is `max_time_step`. No more than `time_step_options.update_period` steps should
be taken, and they should fit in `time_interval`.

# Examples

```jldoctest
julia> using UltraDark: actual_time_step, TimeStepOptions

julia> actual_time_step(0.11, 1, TimeStepOptions())
(0.1, 10)
```
"""
function actual_time_step(max_time_step, time_interval, time_step_options)::Tuple{Float64, Integer}
    max_time_step *= time_step_options.multiplier
    if max_time_step * time_step_options.update_period > time_interval
        num_steps = ceil(time_interval / max_time_step)
        time_interval / num_steps, num_steps
    else
        max_time_step, time_step_options.update_period
    end
end

"""
    max_time_step(grids, a)

Calculate an upper bound on the time step
"""
function max_time_step(grids, a, time_step_options)
    max_time_step_gravity = 2π / maximum(abs.(grids.Φx))
    max_time_step_pressure = 2π * 2 / maximum(grids.k)^2 * a^2  # TODO: cache k_max

    @assert isfinite(max_time_step_gravity)
    @assert isfinite(max_time_step_pressure)

    time_step = min(max_time_step_gravity, max_time_step_pressure)

    time_step
end

"""
    max_time_step(grids, a)

Calculate an upper bound on the time step for `PencilGrids`
"""
function max_time_step(grids::PencilGrids, a, time_step_options)
    # Find maximum on local grid
    local_max_time_step_gravity = 2π / maximum(abs.(grids.Φx))
    local_max_time_step_pressure = 2π * 2 / maximum(grids.k)^2 * a^2  # TODO: cache k_max

    # Maximize over other grids
    max_time_step_gravity = MPI.Allreduce(local_max_time_step_gravity, MPI.MIN, grids.MPI_COMM)
    max_time_step_pressure = MPI.Allreduce(local_max_time_step_pressure, MPI.MIN, grids.MPI_COMM)

    @assert isfinite(max_time_step_gravity)
    @assert isfinite(max_time_step_pressure)

    time_step = min(max_time_step_gravity, max_time_step_pressure)

    time_step
end

ext: jl


text:
function IpoptSolve(nlp;print_level=5)
    tic()
    #solver = IpoptSolver(print_level=5, max_iter=3000, bound_relax_factor=0.0, nlp_scaling_method="none", acceptable_iter=999999)
    solver = IpoptSolver(print_level=print_level, max_iter=3000, bound_relax_factor=0.0, nlp_scaling_method="none", acceptable_iter=999999)

    #,mehrotra_algorithm="yes") #, tol_dual_abs=1e-6)
    #solver = IpoptSolver(print_level=5, tol=1e-8)
    mp = NLPModelsJuMP.NLPtoMPB(nlp, solver)
    MathProgBase.optimize!(mp)
    @show norm(mp.inner.mult_g, Inf)
    #y = MathProgBase.getdual(mp)
    solver = MathProgBase.getrawsolver(mp)
    #finalize(nlp_raw2)
    toc();
end


text:
using JuSwarm, StaticArrays, Test
import Optim

function rosenbrockfunc(x::AbstractVector)
    sum = 0.0
    for i in 1:length(x) - 1
        sum += 100*(x[i + 1] - x[i]^2)^2 + (x[i] - 1)^2
    end
    return sum
end

N = 2
opts = PSOOptions(N; SwarmSize = 1000,
                     HybridOptimizer = Optim.LBFGS(),
                     Display = "none",
                     UseParallel = true)

sol = psoptimize(rosenbrockfunc,opts)

@test sol.fbest <= sol.fbestPSO

@test sol.fbestPSO <= 1.0e-6

@test sol.fbest <= 1.0e-10

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
using AbstractTrees

# TODO these tests are pretty silly in that they now pretty much just
# test code in AbstractTrees.jl. They are still sort of useful though
# in that they test our implementation of `children` indirectly, and I
# an just loath to remove tests in general, so I'm going to leave them
# here for now

const ex = parsehtml("""

abc

                     """)

let res = Any[]
    for node in StatelessBFS(ex.root)
        push!(res, node)
    end
    @assert tag(res[3]) == :body
    @assert tag(res[4]) == :p
    @assert text(last(res)) == "b"
end

let res = Any[]
    for node in PreOrderDFS(ex.root)
        push!(res, node)
    end
    @assert tag(res[3]) == :body
    @assert tag(res[4]) == :p
    @assert text(last(res)) == "c"
end

let res = Any[]
    for node in PostOrderDFS(ex.root)
        push!(res, node)
    end
    @assert tag(res[1]) == :head
    @assert text(res[2]) == "a"
    @assert tag(res[4]) == :strong
    @assert tag(last(res)) == :HTML
end

isp(node::HTMLNode) = isa(node, HTMLElement) && tag(node) == :p

for itr in [PostOrderDFS, PreOrderDFS, StatelessBFS]
    @assert mapreduce(isp,+,itr(ex.root)) == 1
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct GenericResourceExpanded <: SwaggerModel
    plan::Any # spec type: Union{ Nothing, Plan } # spec name: plan
    properties::Any # spec type: Union{ Nothing, Any } # spec name: properties
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    managedBy::Any # spec type: Union{ Nothing, String } # spec name: managedBy
    sku::Any # spec type: Union{ Nothing, Sku } # spec name: sku
    identity::Any # spec type: Union{ Nothing, Identity } # spec name: identity
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    location::Any # spec type: Union{ Nothing, String } # spec name: location
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

    function GenericResourceExpanded(;plan=nothing, properties=nothing, kind=nothing, managedBy=nothing, sku=nothing, identity=nothing, id=nothing, name=nothing, type=nothing, location=nothing, tags=nothing)
        o = new()
        validate_property(GenericResourceExpanded, Symbol("plan"), plan)
        setfield!(o, Symbol("plan"), plan)
        validate_property(GenericResourceExpanded, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        validate_property(GenericResourceExpanded, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(GenericResourceExpanded, Symbol("managedBy"), managedBy)
        setfield!(o, Symbol("managedBy"), managedBy)
        validate_property(GenericResourceExpanded, Symbol("sku"), sku)
        setfield!(o, Symbol("sku"), sku)
        validate_property(GenericResourceExpanded, Symbol("identity"), identity)
        setfield!(o, Symbol("identity"), identity)
        validate_property(GenericResourceExpanded, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(GenericResourceExpanded, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(GenericResourceExpanded, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(GenericResourceExpanded, Symbol("location"), location)
        setfield!(o, Symbol("location"), location)
        validate_property(GenericResourceExpanded, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        o
    end
end # type GenericResourceExpanded

const _property_map_GenericResourceExpanded = Dict{Symbol,Symbol}(Symbol("plan")=>Symbol("plan"), Symbol("properties")=>Symbol("properties"), Symbol("kind")=>Symbol("kind"), Symbol("managedBy")=>Symbol("managedBy"), Symbol("sku")=>Symbol("sku"), Symbol("identity")=>Symbol("identity"), Symbol("id")=>Symbol("id"), Symbol("name")=>Symbol("name"), Symbol("type")=>Symbol("type"), Symbol("location")=>Symbol("location"), Symbol("tags")=>Symbol("tags"))
const _property_types_GenericResourceExpanded = Dict{Symbol,String}(Symbol("plan")=>"Plan", Symbol("properties")=>"Any", Symbol("kind")=>"String", Symbol("managedBy")=>"String", Symbol("sku")=>"Sku", Symbol("identity")=>"Identity", Symbol("id")=>"String", Symbol("name")=>"String", Symbol("type")=>"String", Symbol("location")=>"String", Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ GenericResourceExpanded }) = collect(keys(_property_map_GenericResourceExpanded))
Swagger.property_type(::Type{ GenericResourceExpanded }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_GenericResourceExpanded[name]))}
Swagger.field_name(::Type{ GenericResourceExpanded }, property_name::Symbol) =  _property_map_GenericResourceExpanded[property_name]

function check_required(o::GenericResourceExpanded)
    true
end

function validate_property(::Type{ GenericResourceExpanded }, name::Symbol, val)
    if name === Symbol("kind")
    end
end


text:
@testset "Polygon" begin
    a = SVector{2,Float32}(0, 0)
    b = SVector{2,Float32}(0, 1)
    c = SVector{2,Float32}(1, 1)
    d = SVector{2,Float32}(1, 0)
    p = Polygon([a,b,c,d,a])

    @test boundingbox(p) == BoundingBox(0,0,-Inf,1,1,Inf)

    @test convert2d(p) == p
    @test convert3d(p) == Polygon(map(x->SVector{3,Float32}(x[1], x[2], 0.0), [a,b,c,d,a]))
    @test convert3d(p, 1) == Polygon(map(x->SVector{3,Float32}(x[1], x[2], 1.0), [a,b,c,d,a]))
    @test convert3d(Polygon(map(x->SVector{3,Float32}(x[1], x[2], 1.0), [a,b,c,d,a])), 2) == Polygon{3}([Float32[0.0, 0.0, 2.0], Float32[0.0, 1.0, 2.0], Float32[1.0, 1.0, 2.0], Float32[1.0, 0.0, 2.0], Float32[0.0, 0.0, 2.0]])

    @testset "Point containment and distance" begin
        ls_outer = LineString([SVector(0.0,0.0), SVector(0.0,1.0), SVector(1.0,1.0), SVector(1.0,0.0), SVector(0.0,0.0)])
        ls_inner = LineString([SVector(0.4,0.4), SVector(0.6,0.4), SVector(0.6,0.6), SVector(0.4,0.6), SVector(0.4,0.4)])

        poly1 = Polygon(ls_outer)
        poly2 = Polygon(ls_outer, [ls_inner])

        p_centre = SVector(0.5, 0.5)
        ps_inner_ring = SVector{2,Float64}[[0.3,0.3],[0.4,0.3],[0.5,0.3],[0.6,0.3],[0.7,0.3],[0.7,0.4],[0.7,0.5],[0.7,0.6],[0.7,0.7],[0.6,0.7],[0.5,0.7],[0.4,0.7],[0.3,0.7],[0.3,0.6],[0.3,0.5],[0.3,0.4]]
        ps_outer_ring = SVector{2,Float64}[[-1.0,-1.0],[0.0,-1.0],[0.3,-1.0],[0.4,-1.0],[0.5,-1.0],[0.6,-1.0],[0.7,-1.0],[1.0,-1.0],
                                           [2.0,-1.0],[2.0,0.0],[2.0,0.3],[2.0,0.4],[2.0,0.5],[2.0,0.6],[2.0,0.7],[2.0,1.0],
                                           [2.0,2.0],[1.0,2.0],[0.7,2.0],[0.6,2.0],[0.5,2.0],[0.4,2.0],[0.3,2.0],[0.0,2.0],
                                           [-1.0,2.0],[-1.0,1.0],[-1.0,0.7],[-1.0,0.6],[-1.0,0.5],[-1.0,0.4],[-1.0,0.3],[-1.0,0.0]]

        dists_outer = [sqrt(2), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                       sqrt(2), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                       sqrt(2), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,
                       sqrt(2), 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

        @test p_centre ∈ poly1
        @test p_centre ∉ poly2
        @test distance(p_centre, poly1) == 0.0
        @test distance(p_centre, poly2) ≈ 0.1

        for p in ps_inner_ring
            @test p ∈ poly1
            @test p ∈ poly2
            @test distance(p, poly1) == 0.0
            @test distance(p, poly2) == 0.0
        end

        for (i, p) in enumerate(ps_outer_ring)
            @test p ∉ poly1
            @test p ∉ poly2
            @test distance(p, poly1) == dists_outer[i]
            @test distance(p, poly2) == dists_outer[i]
        end

        poly3 = Polygon(ls_outer) # New instance of poly1
        @test isequal(poly1, poly3)
        @test !isequal(poly1, poly2)
        @test hash(poly1) == hash(poly3)
    end
end


text:
using MarketData
    c1  = convert(TimestampArray, cl[1:5])
    c2  = c1
    cs  = merge(c1, c2)
    csu = merge(c1, c2, sorted=false)

facts("merge works") do

    context("merge and sort") do
        @fact length(cs) --> 10
        @fact cs[2].timestamp --> Date(2000,1,3)
    end

    context("merge and don't sort") do
        @fact length(csu) --> 10
        @fact csu[2].timestamp --> Date(2000,1,4)
    end
end


text:
using PDMats, Distributions
writing_output = false

####################################################################
# Testing Adaptive Φ Solution
####################################################################

file   = JLD2.jldopen("$path/../../reference/helpers_input.jld2", "r")
cloud  = old_to_new_cloud(DSGE.Cloud(read(file, "cloud")))
i_smc  = read(file, "i")
j_smc  = read(file, "j")
ϕ_prop = read(file, "phi_prop")
ϕ_n1   = read(file, "phi_n1")
proposed_fixed_schedule = read(file, "proposed_fixed_schedule")
tempering_target        = read(file, "tempering_target")
resampled_last_period   = read(file, "resampled_last_period")
close(file)

@everywhere Random.seed!(42)
test_ϕ_n, test_resampled_last_period, test_j, test_ϕ_prop = SMC.solve_adaptive_ϕ(cloud,
                                                                proposed_fixed_schedule,
                                                                i_smc, j_smc, ϕ_prop, ϕ_n1,
                                                                tempering_target,
                                                                resampled_last_period)
if writing_output
    jldopen("$path/../../reference/helpers_output.jld2", "w") do file
        write(file, "phi_n", test_ϕ_n)
        write(file, "resampled_last_period", test_resampled_last_period)
        write(file, "j", test_j)
        write(file, "phi_prop", test_ϕ_prop)
    end
    jldopen("$path/../../reference/helpers_output.jld2", "w") do file
        write(file, "phi_n", test_ϕ_n)
        write(file, "resampled_last_period", test_resampled_last_period)
        write(file, "j", test_j)
        write(file, "phi_prop", test_ϕ_prop)
    end
end

file = JLD2.jldopen("$path/../../reference/helpers_output.jld2", "r")
saved_ϕ_n = read(file, "phi_n")
saved_resampled_last_period = read(file, "resampled_last_period")
saved_j = read(file, "j")
saved_ϕ_prop = read(file, "phi_prop")
close(file)

####################################################################
@testset "Solve Adaptive Φ" begin
    @test test_ϕ_n ≈ saved_ϕ_n
    @test test_resampled_last_period == saved_resampled_last_period
    @test test_j == saved_j
    @test test_ϕ_prop == saved_ϕ_prop
end

####################################################################
# Testing MvNormal Mixture Draw Function
####################################################################
#=
file = JLD2.jldopen("reference/mutation_inputs.jld2")
saved_particles = read(file, "particles")
blocks_free     = read(file, "blocks_free")
blocks_all      = read(file, "blocks_all")
d               = read(file, "d")
α         = read(file, "α")
c         = read(file, "c")
close(file)
para = saved_particles[1].value
for (block_f, block_a) in zip(blocks_free, blocks_all)
    para_subset = para[block_a]
    d_subset    = MvNormal(d.μ[block_f], d.Σ.mat[block_f, block_f])
    test_θ_new, test_new_mix_density, test_old_mix_density = DSGE.mvnormal_mixture_draw(para_subset,
                                                                                d_subset; c=c, α=α)
    JLD2.jldopen("reference/mvnormal_inputs.jld2", true, true, true, IOStream) do file
        write(file, "para_subset", para_subset)
        write(file, "d_subset", d_subset)
        write(file, "α", α)
        write(file, "c", c)
    end
    JLD2.jldopen("reference/mvnormal_output.jld2", true, true, true, IOStream) do file
        write(file, "θ_new", test_θ_new)
        write(file, "new_mix_density", test_new_mix_density)
        write(file, "old_mix_density", test_old_mix_density)
    end
end
=#

file = JLD2.jldopen("reference/mvnormal_inputs.jld2")
para_subset = read(file, "para_subset")
d_subset    = read(file, "d_subset")
α           = read(file, "α")
c           = read(file, "c")
close(file)

test_θ_new, test_new_mix_density, test_old_mix_density = DSGE.mvnormal_mixture_draw(para_subset,
                                                                                d_subset; c=c, α=α)

#=JLD2.jldopen("reference/mvnormal_output.jld2", true, true, true, IOStream) do file
    write(file, "θ_new", test_θ_new)
    write(file, "new_mix_density", test_new_mix_density)
    write(file, "old_mix_density", test_old_mix_density)
end=#

file = JLD2.jldopen("reference/mvnormal_output.jld2")
    saved_θ_new = read(file, "θ_new")
    saved_new_mix_density = read(file, "new_mix_density")
    saved_old_mix_density = read(file, "old_mix_density")
close(file)

####################################################################
@testset "MvNormal Mixture Draw" begin
    @test test_θ_new == saved_θ_new
    @test test_new_mix_density == saved_new_mix_density
    @test test_old_mix_density == saved_old_mix_density
end

####################################################################
# Testing ESS Computation
####################################################################
#=
file = jldopen("reference/smc_sw_cloud_fix=true_blocks=3.jld2")
cloud = read(file, "cloud")
current_weights = read(file, "w")[:,3]
close(file)

n_part    = length(cloud.particles)
loglh     = [cloud.particles[i].loglh for i=1:n_part]
old_loglh = [cloud.particles[i].old_loglh for i=1:n_part]
ϕ_n       = 9.25022e-6
ϕ_n1      = 2.15769e-6

ess = DSGE.compute_ESS(loglh, current_weights, ϕ_n, ϕ_n1)

jldopen("reference/ess_inputs.jld2", true, true, true, IOStream) do file
    write(file, "loglh", loglh)
    write(file, "current_weights", current_weights)
    write(file, "ϕ_n", ϕ_n)
    write(file, "ϕ_n1", ϕ_n1)
    write(file, "old_loglh", old_loglh)
end

jldopen("reference/ess_output.jld2", true, true, true, IOStream) do file
    write(file, "ess", ess)
end
=#
file = JLD2.jldopen("reference/ess_inputs.jld2")
loglh = read(file, "loglh")
current_weights = read(file, "current_weights")
ϕ_n = read(file, "ϕ_n")
ϕ_n1 = read(file, "ϕ_n1")
close(file)

file = JLD2.jldopen("reference/ess_output.jld2")
saved_ESS = read(file, "ess")
close(file)

test_ESS = SMC.compute_ESS(loglh, current_weights, ϕ_n, ϕ_n1)
####################################################################
@testset "Compute ESS" begin
    @test test_ESS ≈ saved_ESS
end


text:
module BioServices

include("common.jl")
include("eutils.jl")
include("umls.jl")
include("gggenome.jl")
include("ebiproteins.jl")

end

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
cd(@__DIR__)

using Pkg

using Test, TestSetExtensions, SafeTestsets, Logging
using Genie

Logging.global_logger(NullLogger())

@testset ExtendedTestSet "Genie tests" begin
  @includetests ARGS
end


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

using Markdown
using InteractiveUtils

# ╔═╡ f79f691c-8e49-11eb-0aff-714a3f5248d5
using PlutoUI, Distributions, Plots

# ╔═╡ 70426596-8e49-11eb-2455-a9818cad032b
md"""
# Expectimax
"""

# ╔═╡ 543cd594-8e4a-11eb-1101-2d188eda8b88
md"""
## Application 1 (5.16)
"""

# ╔═╡ 64fbc4f6-8e4a-11eb-36d8-29ff6f774100
PlutoUI.LocalResource("./img/chancetree.png", :width=>400)

# ╔═╡ e1dffc86-8e49-11eb-10b7-e795ee87a053
md"""
This question considers pruning in games with chance nodes. Consider the complete game tree for a trivial game shown above.

Assume that the leaf nodes are to be evaluated in left-to-right order, and that before a leaf node is evaluated, we know nothing about its value, i.e. the range of possible values is $-\infty$ to $-\infty$.

1. Copy the figure, mark the value of all the internal nodes, and indicate the best move at the root with an arrow.
2. Given the values of the first six leaves, do we need to evaluate the seventh and eighth leaves?
3. Given the values of the first seven leaves, do we need to evaluate the eighth leaf? Explain your answers.
4. Suppose the leaf node values are known to lie between –2 and 2 inclusive. After the first two leaves are evaluated, what is the value range for the left-hand chance node?
5. Circle all the leaves that need not be evaluated under the assumption in (4). Explain your answers."""

# ╔═╡ ed342f68-8e4a-11eb-0173-a3cd935b42eb
md"""
## Application 2 (5.21)

In the following, a “max” tree consists only of max nodes, whereas an “expectimax” tree consists of a max node at the root with alternating layers of chance and max nodes. At chance nodes, all outcome probabilities are nonzero. The goal is to find the value of the root with a bounded-depth search. An illustration is shown below.
"""

# ╔═╡ 96eb968e-8e4b-11eb-0d82-f7cf952a4b55
PlutoUI.LocalResource("./img/expectimaxtree.png",:width=>200)

# ╔═╡ 11d98910-8e4b-11eb-2097-69a6a67b97b0
begin
	PlutoUI.LocalResource("./img/maxtree.png",:width=>200)
end

# ╔═╡ a719f96a-8e4b-11eb-1077-8be6a992691a
md"""
1. Assuming that leaf values are finite but unbounded, is pruning (as in alpha–beta) ever possible in a max tree? Give an example, or explain why not.
2. Is pruning ever possible in an expectimax tree under the same conditions? Give an example, or explain why not.
3. If leaf values are constrained to be non-negative, is pruning ever possible in a max tree? Give an example, or explain why not.
4. If leaf values are constrained to be non-negative, is pruning ever possible in an expectimax tree? Give an example, or explain why not.
5. If leaf values are constrained to be in the range $[0, 1]$, is pruning ever possible in a max tree? Give an example, or explain why not.
6. If leaf values are constrained to be in the range $[0, 1]$, is pruning ever possible in an expectimax tree? Give an example (qualitatively different from your example in (5), if any) or explain why not.
7. Consider the outcomes of a chance node in an expectimax tree. Which of the following evaluation orders is most likely to yield pruning opportunities: (i) Lowest probability first; (ii) Highest probability first; (iii) Doesn’t make any difference?
"""

# ╔═╡ 282f47d0-8e4c-11eb-3e1e-d15a008161b3
md"""
## Application 3 - game
* each  player gets to play $n$ times. At each round a player can decide if he wants to add one or two to his score.
* after each player, a coin is flipped. If it's heads, one of your own points goes to the other player, if not, nothing changes. Assume that the coin is not biased and you have a 50/50 chance between heads and tails.
* the value in a leaf node is given by the difference between max and min.

This game has been implemented following the layout you already know from previous sessions. Try to implement the expecitiminiax algorithm and use it on this game. For a game with limited length you can draw the game tree and reason on it.

"""

# ╔═╡ 387aaea8-8e4d-11eb-0abb-93c31193af2f
begin
	abstract type AbstractGame end
	abstract type AbstractGameState end

	"""
		ProbGameState

	Representation of a probabilistic game state

	# Arguments
	`turn::String`: current player
	`score::Tuple{Int64,Int64}`: tuple holding the scores for Max and Min
	`prevplay::String`: previous 'human' player (i.e. not a chance node)
	`nturns::Int64`: number of turns that have been completed
	"""
	struct ProbGameState <: AbstractGameState
		turn::String
		score::Tuple{Int64,Int64}
		prevplay::String
		nturns::Int64
	end

	function Base.show(io::IO, s::ProbGameState)
    	print(io, "ProbGameState\n-------------\nCurrent player: $(s.turn)\nscore: $(s.score)\nPrevious 'human' player: $(s.prevplay)\nNumber of turns played: $(round(Int,s.nturns/2))")
	end

	function ProbGameState()
		return ProbGameState("max",(0,0),"",0)
	end

	"""
		ProbGame

	Representation of a probabilistic game.

	# Notes
	Game starts with max

	# Arguments
	`initial::ProbGameState`: initial state
	`n::Int64`: number of turns for each player
	"""
	struct ProbGame <: AbstractGame
		initial::ProbGameState
		n::Int64
	end

	function ProbGame(initial::ProbGameState=ProbGameState(); n=3)
		return ProbGame(initial, n)
	end

	function Base.show(io::IO, g::ProbGame)
		print(io, "ProbGame with a depth of $(g.n)")
	end

	"""
		actions(game, state)

	Retrieve possible actions from the current state of the game.
	"""
	function actions(game::ProbGame, state::ProbGameState)
		if state.turn == "max"
			return [(1,0), (2,0)]
		elseif state.turn == "min"
			return [(0,1), (0,2)]
		elseif state.turn == "chance"
			if state.prevplay == "max"
				return [(-1,1), (0,0)]
			else
				return [(1,-1), (0,0)]
			end
		end
	end

	"""
		result(game, state, move)

	Given a game, its state and the move that is played, return the new state.
	"""
	function result(game::ProbGame, state::ProbGameState, move::Tuple{Int64,Int64})
		prevplay = state.prevplay
		nturns = state.nturns
		# determine next player
		if state.turn ≠ "chance"
			nturns += 1
			nextplayer = "chance"
			if state.turn == "max"
				prevplay = "max"
			else
				prevplay = "min"
			end
		else
			if state.turn == "chance" && state.prevplay == "max"
				nextplayer = "min"
			else
				nextplayer = "max"
			end
		end
		# determine new score
		newscore = state.score .+ move

		return ProbGameState(nextplayer, newscore, prevplay, nturns)
	end

	"""
		terminal_test(game, state)

	Verify a if game has ended
	"""
	function terminal_test(game::ProbGame, state::ProbGameState)
		return game.n == state.nturns/2
	end

	"""
		utility(game, state)

	Returns the utility for a player of a game in a state
	"""
	function utility(game::ProbGame, state::ProbGameState)
		r = state.score[1] - state.score[2]
		if state.prevplay == "max"
			return -r
		elseif state.prevplay == "min"
			return r
		end
	end

	"""
		to_move(game, state)

	Determine whos turn it is.
	"""
	function to_move(game::ProbGame, state::ProbGameState)
		return state.turn
	end

	"""
		prob(move::Tuple{Int64,Int64})

	Probablity of a given move (all equal to 0.5, cf. problem description)
	"""
	function prob(move::Tuple{Int64,Int64})
		return 0.5
	end
end

# ╔═╡ d3b39b3c-8e4d-11eb-123f-0307f1775d20
ProbGameState()

# ╔═╡ 4e82fd7a-8e4f-11eb-1e96-4daacfbfa65d
ProbGame(n=1)

# ╔═╡ 53e38136-8e4f-11eb-046f-95806459c80e
## illustration of a small game
let
	@info "Play a small game with random actions"
	# game setup
	game = ProbGame(n=1)
	state = game.initial
	# keep playing until the game is finished
	while !terminal_test(game,state)
		@show state
		action = rand(actions(game, state))
		@show action
    	state = result(game, state, action)
	end

	@show state
	@show utility(game, state)
end

# ╔═╡ 33b23a08-8e52-11eb-2c06-dfbd8939d670
md"""
Below you can find an implementation of the EXPECTIMINIMAX algorithm
"""

# ╔═╡ 9e1a905c-8e52-11eb-0aa0-89aaf3701a29
begin
	"""
		expectiminmaxchoice(game, state)

	Determine the best action (move) for a stochastic game in a state.
	"""
	function expectiminmaxchoice(game::T, state::M) where {T <: AbstractGame, M <: AbstractGameState}
		# list of tuples (move, expexted value)
		movescores = [(move, expectiminmaxvalue(game, result(game, state, move))) for move in actions(game, state)]

		if state.turn == "max"
			return sort(movescores, by=x->x[2], rev=true)[1][1]
		else
			return sort(movescores, by=x->x[2], rev=false)[1][1]
		end
	end

	"""
		expectiminmaxvalue(game, state)

	Determine the expectiminmax value for a game in a state.
	"""
	function expectiminmaxvalue(game::T, state::M) where {T <: AbstractGame, M <: AbstractGameState}
		player = state.turn

		# If we are in a terminal state, we return the utility of the player
		if (terminal_test(game, state))
			r = utility(game, state)
			return r
		end

		# otherwise, we keep exploring the three
		if player == "max"
			return maximum([expectiminmaxvalue(game, result(game, state, move)) for move in actions(game, state)])
		elseif player == "min"
			return minimum([expectiminmaxvalue(game, result(game, state, move)) for move in actions(game, state)])
		elseif player == "chance"
			r = sum([prob(move)*expectiminmaxvalue(game, result(game, state, move)) for move in actions(game,state)])
			return r
		end
	end
end

# ╔═╡ 633d00e0-8e53-11eb-3f5c-55bb14c56012
let
	@info "Determine best move for game of depth 1"
	# game setup
	game = ProbGame(n=1)
	@show game
	# best move for player (and expected utility)
	bestmove = expectiminmaxchoice(game,game.initial)
	@show bestmove
end

# ╔═╡ b5ce4ea0-8e54-11eb-36f1-bd256cc6cef7
let
	@info "illustration of a game of depth 1 starting with min"
	# game setup
	game = ProbGame(ProbGameState("min", (0,0), "", 0),n=1)
	@show game
	# best move for player (and expected utility)
	bestmove = expectiminmaxchoice(game,game.initial)
	@show bestmove
end

# ╔═╡ e5168334-8e55-11eb-32aa-056184ae94ed
md"""
# Making simple decision
*Note*: In section 16.3.4, a lot is said about human behavior. One of the authors quoted is Daniel Kahneman. He wrote an entire book about how easy it is to trick the mind (e.g. by framing, anchoring) and how hard is for us to make unbiased decisions. If you are interested, the book is called ["Thinking, Fast and Slow"](https://en.wikipedia.org/wiki/Thinking,_Fast_and_Slow) and is definitely worth a read.

## Application 1 (16.3)
In 1713, Nicolas Bernoulli stated a puzzle, now called the St. Petersburg paradox, which works as follows. You have the opportunity to play a game in which a fair coin is tossed repeatedly until it comes up heads. If the first heads appears on the n$^{th}$ toss, you win $2^n$ dollars.
1. Show that the expected monetary value of this game is infinite.
2. How much would you, personally, pay to play the game?
3. Nicolas’s cousin Daniel Bernoulli resolved the apparent paradox in 1738 by suggesting that the utility of money, $U(S_n)$, is measured on a logarithmic scale: $U(S_n) = a \log_2(m)+b$, where $S_n$ is the state of having $m$ dollars. What is the expected utility of the game under this assumption? Suppose you use the entirety of your capital to play this game.
4. How would you describe a rational player in this case? What is the maximum amount that would be rational to pay to play the game, assuming that (1) one's initial wealth is $k$ (2) one uses a quantity $c$ to play the game (3) a single game costs $1$ dollar?
"""

# ╔═╡ 6b0aead6-8e57-11eb-0ed3-f728f1bf72dd

# ╔═╡ 194dde60-8e57-11eb-0719-f18672c2fc48
md"""
## Application 2 (16.10) - Expected utility and post-decision disappointment
Let $X_1,\dots , X_k$ be continuous variables, independently distributed according to the same probability density function $f_X(x)$. Show that the density function for $\max \{ X_1,\dots,X_k \} $ is given by $kf_X(x)(F_X(x))^{k−1}$, where $F_X(x)$ is the cumulative distribution for $f_X$.

Some hints:
* What do you remember from independent events in terms of probability?
* What is the relation between $F_X(x)$ and $f_X(x)$?

Once you find the theoretical distribution, run a simulation to illustrate the bias that is introduced by using estimates of the expected utilty (cf. probability & statistics + modelling & simulation course).
1. try to experimentally reproduce figure 16.3.
2. use another distribution, look at the effect.
"""

# ╔═╡ 3e8fa6b2-8e58-11eb-1039-a9f966b12266
begin
	"""
		makedata(d::Distribution, k::Int64, n=1000)

	Generate a sample of length n from a distribution d for the maximum of k estimates

	"""
	function makedata(d::Distribution, k::Int64, n=1000)
		return maximum.([rand(d,k) for _ in 1:n])
	end

	"""
		makeplot(dist::Distribution, kvals::Array{Int64,1}=[3, 10, 30, 100], n=1000)

	Generates an plot of the distribution of the maximum of k utility estimates (for k in kval).

	Both the experimental and the theoretical distribution are shown.

	"""
	function makeplot(dist::Distribution, kvals::Array{Int64,1}=[3, 10, 30, 100], n=1000)
		f(x) = pdf(dist,x)
		F(x) = cdf(dist,x)
		Plots.histogram()
		for k in kvals
			s = makedata(dist,k,1000)
			x = range(minimum(s),maximum(s),length=500)
			Plots.histogram!(s, label="k = $(k)", alpha=0.6, normalize=:pdf)
			Plots.plot!(x,k*f.(x) .* F.(x) .^(k-1), label="", linecolor=:black, lw=4, ls=:dash)
		end
		Plots.xlabel!("x")
		Plots.ylabel!("experimental PDF")
		Plots.title!("Distribution of maxima for\n $(dist)")
	end
end

# ╔═╡ 64859af2-8e58-11eb-24c3-93dcbcf1901d
makeplot(Normal())

# ╔═╡ 6af7da96-8e58-11eb-376d-7b587a5c31f0
makeplot(Exponential(10))

# ╔═╡ Cell order:
# ╠═f79f691c-8e49-11eb-0aff-714a3f5248d5
# ╟─70426596-8e49-11eb-2455-a9818cad032b
# ╟─543cd594-8e4a-11eb-1101-2d188eda8b88
# ╟─64fbc4f6-8e4a-11eb-36d8-29ff6f774100
# ╟─e1dffc86-8e49-11eb-10b7-e795ee87a053
# ╟─ed342f68-8e4a-11eb-0173-a3cd935b42eb
# ╟─96eb968e-8e4b-11eb-0d82-f7cf952a4b55
# ╟─11d98910-8e4b-11eb-2097-69a6a67b97b0
# ╟─a719f96a-8e4b-11eb-1077-8be6a992691a
# ╟─282f47d0-8e4c-11eb-3e1e-d15a008161b3
# ╠═387aaea8-8e4d-11eb-0abb-93c31193af2f
# ╠═d3b39b3c-8e4d-11eb-123f-0307f1775d20
# ╠═4e82fd7a-8e4f-11eb-1e96-4daacfbfa65d
# ╠═53e38136-8e4f-11eb-046f-95806459c80e
# ╟─33b23a08-8e52-11eb-2c06-dfbd8939d670
# ╠═9e1a905c-8e52-11eb-0aa0-89aaf3701a29
# ╠═633d00e0-8e53-11eb-3f5c-55bb14c56012
# ╠═b5ce4ea0-8e54-11eb-36f1-bd256cc6cef7
# ╟─e5168334-8e55-11eb-32aa-056184ae94ed
# ╠═6b0aead6-8e57-11eb-0ed3-f728f1bf72dd
# ╟─194dde60-8e57-11eb-0719-f18672c2fc48
# ╟─3e8fa6b2-8e58-11eb-1039-a9f966b12266
# ╠═64859af2-8e58-11eb-24c3-93dcbcf1901d
# ╠═6af7da96-8e58-11eb-376d-7b587a5c31f0

max_stars_repo_path: PS/PS07 - expectimax and utilities.jl

max_issues_repo_path: PS/PS07 - expectimax and utilities.jl

max_forks_repo_path: PS/PS07 - expectimax and utilities.jl


text:
import REPL.REPLCompletions, FuzzyCompletions

struct CompletionSuggestion
  replacementPrefix::String
  # suggestion body
  text::String
  type::String
  icon::String
  rightLabel::String
  leftLabel::String
  description::String
  descriptionMoreURL::String
  # for `getSuggestionDetailsOnSelect` API
  detailtype::String
end

CompletionSuggestion(prefix, text, type, icon; rl = "", ll = "", desc = "", url = "", detail = "") =
  CompletionSuggestion(prefix, text, type, icon, rl, ll, desc, url, detail)

# as an heuristic, suppress completions if there are over 500 completions,
# ref: currently `completions("", 0, Main)` returns 1098 completions as of v1.5
const SUPPRESS_COMPLETION_THRESHOLD = 500

# autocomplete-plus only shows at most 200 completions
# ref: https://github.com/atom/autocomplete-plus/blob/master/lib/suggestion-list-element.js#L49
const MAX_COMPLETIONS = 200

const DESCRIPTION_LIMIT = 200

# threshold up to which METHOD_COMPLETION_CACHE can get fat, to make sure it won't eat up memory
# NOTE: with 2000 elements it used ~10MiB on my machine.
const MAX_METHOD_COMPLETION_CACHE = 3000

# baseline completions
# --------------------

handle("completions") do data
  @destruct [
    # general
    line,
    mod || "Main",
    path || nothing,
    # local context
    context || "",
    row || 1,
    startRow || 0,
    column || 1,
    # configurations
    is_fuzzy || true,
    force || false
  ] = data

  adapter = is_fuzzy ? fuzzycompletionadapter : replcompletionadapter
  return adapter(
    # general
    line, mod,
    # local context
    context, row - startRow, column,
    # configurations
    force
  )::Vector{CompletionSuggestion}
end

# NOTE:
# `replcompletionadapter` and `fuzzycompletionadapter` are really similar
# and maybe better to be refactored into a single function,
# but we define them as separate functions for type-stabilities for now

@inbounds function replcompletionadapter(
  # general
  line, m = "Main",
  # local context
  context = "", row = 1, column = 1,
  # configurations
  force = false
)
  mod = getmodule(m)

  cs, replace, shouldcomplete = try
    REPLCompletions.completions(line, lastindex(line), mod)
  catch err
    # might error when e.g. type inference fails
    REPLCompletions.Completion[], 1:0, false
  end
  prefix = line[replace]

  # suppress completions if there are still too many of them i.e. when invoked with `$|`, `(|`, etc.
  # XXX: the heuristics below mayn't be so robust to work for all the cases.
  !force && (isempty(prefix) && length(cs) > SUPPRESS_COMPLETION_THRESHOLD) && (cs = REPLCompletions.Completion[])

  # initialize suggestions with local completions so that they show up first
  comps = if force || !isempty(prefix)
    filter!(let p = prefix
      c -> startswith(c.text, p)
    end, localcompletions(context, row, column, prefix))
  else
    CompletionSuggestion[]
  end

  # FIFO cache refreshing
  if length(METHOD_COMPLETION_CACHE) ≥ MAX_METHOD_COMPLETION_CACHE
    for k in collect(keys(METHOD_COMPLETION_CACHE))[1:MAX_METHOD_COMPLETION_CACHE÷2]
      delete!(METHOD_COMPLETION_CACHE, k)
    end
  end

  cs = cs[1:min(end, MAX_COMPLETIONS - length(comps))]
  REPLCompletions.afterusing(line, Int(first(replace))) && filter!(!ispkgcomp, cs) # need `Int` for correct dispatch on x86
  append!(comps, completion.(Ref(mod), cs, prefix))

  return comps
end

@inbounds function fuzzycompletionadapter(
  # general
  line, m = "Main",
  # local context
  context = "", row = 1, column = 1,
  # configurations
  force = false
)
  mod = getmodule(m)

  cs, replace, shouldcomplete = try
    FuzzyCompletions.completions(line, lastindex(line), mod)
  catch err
    # might error when e.g. type inference fails
    FuzzyCompletions.Completion[], 1:0, false
  end
  prefix = line[replace]

  if !force
    filter!(c -> FuzzyCompletions.score(c) ≥ 0, cs) # filter negative scores

    # suppress completions if there are still too many of them i.e. when invoked with `$|`, `(|`, etc.
    # XXX: the heuristics below mayn't be so robust to work for all the cases.
    (isempty(prefix) && length(cs) > SUPPRESS_COMPLETION_THRESHOLD) && (cs = FuzzyCompletions.Completion[])
  end

  # initialize suggestions with local completions so that they show up first
  comps = localcompletions(context, row, column, prefix)
  if !force
    filter!(let p = prefix
      # NOTE: let's be a bit strict on local completions so that we avoid verbose completions, e.g. troublemaker cases
      comp -> FuzzyCompletions.fuzzyscore(p, comp.text) > 0
    end, comps)
  end

  # FIFO cache refreshing
  if length(METHOD_COMPLETION_CACHE) ≥ MAX_METHOD_COMPLETION_CACHE
    for k in collect(keys(METHOD_COMPLETION_CACHE))[1:MAX_METHOD_COMPLETION_CACHE÷2]
      delete!(METHOD_COMPLETION_CACHE, k)
    end
  end

  cs = cs[1:min(end, MAX_COMPLETIONS - length(comps))]
  REPLCompletions.afterusing(line, Int(first(replace))) && filter!(!ispkgcomp, cs) # need `Int` for correct dispatch on x86
  append!(comps, completion.(Ref(mod), cs, prefix))

  return comps
end

const ispkgcomp = Base.Fix2(isa, REPLCompletions.PackageCompletion)

# for functions below works both for REPLCompletions and FuzzyCompletions modules
for c in [:KeywordCompletion, :PathCompletion, :ModuleCompletion, :PackageCompletion,
          :PropertyCompletion, :FieldCompletion, :MethodCompletion, :BslashCompletion,
          :ShellCompletion, :DictCompletion]
  eval(:(const $c = Union{REPLCompletions.$c, FuzzyCompletions.$c}))
end
completion_text(c::REPLCompletions.Completion) = REPLCompletions.completion_text(c)
completion_text(c::FuzzyCompletions.Completion) = FuzzyCompletions.completion_text(c)

completion(mod, c, prefix) =
  CompletionSuggestion(
    prefix, completiontext(c), completiontype(c), completionicon(c);
    rl = completionmodule(mod, c), ll = completionreturntype(c), url = completionurl(c), detail = completiondetailtype(c)
  )

const MethodCompletionInfo = NamedTuple{(:f,:m,:tt),Tuple{Any,Method,Type}}
const MethodCompletionDetail = NamedTuple{(:rt,:desc),Tuple{String,String}}
"""
    METHOD_COMPLETION_CACHE::OrderedDict{String,Union{MethodCompletionInfo,MethodCompletionDetail}}
    MethodCompletionInfo = NamedTuple{(:f,:m,:tt),Tuple{Any,Method,Type}}
    MethodCompletionDetail = NamedTuple{(:rt,:desc),Tuple{String,String}}

[`OrderedCollections.OrderedDict`](@ref) that caches details of [`MethodCompletion`](@ref)s:
- result of a return type inference
- documentation description

This cache object has the following structure:
Keys are hash `String` of a `MethodCompletion` object,
  which can be used for cache detection or refered lazily by
  [autocompluete-plus's `getSuggestionDetailsOnSelect` API]
  (https://github.com/atom/autocomplete-plus/wiki/Provider-API#defining-a-provider)
Values can be either of:
- `MethodCompletionInfo`: keeps (temporary) information of this `MethodCompletion`
  + `f`: a function object
  + `m`: a `Method` object
  + `tt`: an input `Tuple` type
- `MethodCompletionDetail`: stores the lazily computed detail of this `MethodCompletion`
  + `rt`: `String` representing a return type of this method call
  + `desc`: `String` representing documentation of this method

If the second element is still a `MethodCompletionInfo` instance, that means its details
  aren't lazily computed yet.
If the second element is already a `MethodCompletionDetail` then we can reuse this cache.

!!! note
    If a method gets updated (i.e. redefined etc), key hash for the `MethodCompletion`
      will be different from the previous one, so it's okay if we just rely on the hash key.

!!! warning
    Within the current implementation, we can't reflect changes that happens in backedges.
"""
const METHOD_COMPLETION_CACHE = OrderedDict{String,Union{MethodCompletionInfo,MethodCompletionDetail}}()

function completion(mod, c::MethodCompletion, prefix)
  k = repr(hash(c))

  m = c.method
  v = get(METHOD_COMPLETION_CACHE, k, nothing)
  if v !== nothing
    if v isa MethodCompletionDetail
      # cache found
      rt, desc = v
      dt = ""
    else
      # MethodCompletion information has already been stored, but lazy return type inference and
      # description completement has not been done yet
      rt, desc = "", ""
      dt = k
    end
  else
    # store MethodCompletion information for lazy return type inference and description completement
    info = MethodCompletionInfo((
      f = c.func,
      m = m,
      tt = c.input_types
    ))
    METHOD_COMPLETION_CACHE[k] = info
    rt, desc = "", ""
    dt = k
  end

  return CompletionSuggestion(
    prefix, completiontext(c), completiontype(c), completionicon(c);
    rl = completionmodule(mod, c), ll = rt, desc = desc, url = completionurl(c), detail = dt
  )
end

completiontext(c) = completion_text(c)
completiontext(c::MethodCompletion) = begin
  ct = completion_text(c)
  m = match(r"^(.*) in .*$", ct)
  m isa Nothing ? ct : m[1]
end
completiontext(c::DictCompletion) = rstrip(completion_text(c), [']', '"'])
completiontext(c::PathCompletion) = rstrip(completion_text(c), '"')

completionreturntype(c) = ""
completionreturntype(c::PropertyCompletion) = begin
  isdefined(c.value, c.property) || return ""
  shortstr(typeof(getproperty(c.value, c.property)))
end
completionreturntype(c::FieldCompletion) =
  shortstr(fieldtype(c.typ, c.field))
completionreturntype(c::DictCompletion) =
  shortstr(valtype(c.dict))
completionreturntype(::PathCompletion) = "Path"

completionurl(c) = ""
completionurl(c::ModuleCompletion) = begin
  mod, name = c.parent, c.mod
  val = getfield′(mod, name)
  if val isa Module # module info
    urimoduleinfo(parentmodule(val) == val || val ∈ (Base, Core) ? name : "$mod.$name")
  else
    uridocs(mod, name)
  end
end
completionurl(c::MethodCompletion) = uridocs(c.method.module, c.method.name)
completionurl(c::PackageCompletion) = urimoduleinfo(c.package)
completionurl(c::KeywordCompletion) = uridocs("Main", c.keyword)

completionmodule(mod, c) = shortstr(mod)
completionmodule(mod, c::ModuleCompletion) = shortstr(c.parent)
completionmodule(mod, c::MethodCompletion) = shortstr(c.method.module)
completionmodule(mod, c::FieldCompletion) = shortstr(c.typ) # predicted type
completionmodule(mod, ::KeywordCompletion) = ""
completionmodule(mod, ::PathCompletion) = ""

completiontype(c) = "variable"
completiontype(c::ModuleCompletion) = begin
  ct = completion_text(c)
  ismacro(ct) && return "snippet"
  mod, name = c.parent, Symbol(ct)
  val = getfield′(mod, name)
  wstype(mod, name, val)
end
completiontype(::MethodCompletion) = "method"
completiontype(::PackageCompletion) = "import"
completiontype(::PropertyCompletion) = "property"
completiontype(::FieldCompletion) = "property"
completiontype(::DictCompletion) = "property"
completiontype(::KeywordCompletion) = "keyword"
completiontype(::PathCompletion) = "path"

completionicon(c) = ""
completionicon(c::ModuleCompletion) = begin
  ismacro(c.mod) && return "icon-mention"
  mod, name = c.parent, Symbol(c.mod)
  val = getfield′(mod, name)
  wsicon(mod, name, val)
end
completionicon(::DictCompletion) = "icon-key"
completionicon(::PathCompletion) = "icon-file"

completiondetailtype(c) = ""
completiondetailtype(::ModuleCompletion) = "module"
completiondetailtype(::KeywordCompletion) = "keyword"

function localcompletions(context, row, col, prefix)
  ls = locals(context, row, col)
  lines = split(context, '\n')
  return localcompletion.(ls, prefix, Ref(lines))
end
function localcompletion(l, prefix, lines)
  desc = if l.verbatim == l.name
    # show a line as is if ActualLocalBinding.verbatim is not so informative
    lines[l.line]
  else
    l.verbatim
  end |> s -> strlimit(s, DESCRIPTION_LIMIT)
  type = (type = static_type(l)) == "variable" ? "attribute" : type
  icon = (icon = static_icon(l)) == "v" ? "icon-chevron-right" : icon
  return CompletionSuggestion(
    prefix, l.name, type, icon;
    rl = l.root, desc = desc, detail = ""
  )
end

# completion details
# ------------------

function completiondetail!(comp)
  dt = comp["detailtype"]::String
  isempty(dt) && return comp
  word = comp["text"]::String

  if dt == "module"
    mod = comp["rightLabel"]::String
    completiondetail_module!(comp, mod, word)
  elseif dt == "keyword"
    completiondetail_keyword!(comp, word)
  else # detail for method completion
    completiondetail_method!(comp, dt)
  end

  comp["detailtype"] = "" # may not be needed, but empty this to make sure any further detail completion doesn't happen
  return comp
end

handle(completiondetail!, "completiondetail")

function completiondetail_module!(comp, mod, word)
  mod = getmodule(mod)
  cangetdocs(mod, word) || return
  comp["description"] = completiondescription(getdocs(mod, word))
end

completiondetail_keyword!(comp, word) =
  comp["description"] = completiondescription(getdocs(Main, word))

using JuliaInterpreter: sparam_syms
using Base.Docs

function completiondetail_method!(comp, k)
  v = get(METHOD_COMPLETION_CACHE, k, nothing)
  v === nothing && return # shouldn't happen but just to make sure and help type inference
  if v isa MethodCompletionDetail
    # NOTE
    # This path sometimes happens maybe because of some lags between Juno and getSuggestionDetailsOnSelect API,
    # especially just after we booted Julia, and there are not so much inference caches in the Julia internal.
    # But since details are computed, we can just use here as well
    comp["leftLabel"] = v.rt
    comp["description"] = v.desc
    return
  end
  f, m, tt = v.f, v.m, v.tt

  # return type inference
  rt = rt_inf(f, m, Base.tuple_type_tail(tt))
  comp["leftLabel"] = rt

  # description for this method
  mod = m.module
  fsym = m.name
  desc = if cangetdocs(mod, fsym)
    try
      docs = Docs.doc(Docs.Binding(mod, fsym), Base.tuple_type_tail(m.sig))
      completiondescription(docs)
    catch
      ""
    end
  else
    ""
  end
  comp["description"] = desc

  # update method completion cache with the results

  METHOD_COMPLETION_CACHE[k] = MethodCompletionDetail((rt = rt, desc = desc))
end

function rt_inf(@nospecialize(f), m, @nospecialize(tt::Type))
  try
    world = typemax(UInt) # world age

    # first infer return type using input types
    # NOTE:
    # since input types are all concrete, the inference result from them is the best what we can get
    # so here we eagerly respect that if inference succeeded
    if !isempty(tt.parameters)
      inf = Core.Compiler.return_type(f, tt, world)
      inf ∉ (nothing, Any, Union{}) && return shortstr(inf)
    end

    # sometimes method signature can tell the return type by itself
    sparams = Core.svec(sparam_syms(m)...)
    inf = @static if isdefined(Core.Compiler, :NativeInterpreter)
      Core.Compiler.typeinf_type(Core.Compiler.NativeInterpreter(), m, m.sig, sparams)
    else
      Core.Compiler.typeinf_type(m, m.sig, sparams, Core.Compiler.Params(world))
    end
    inf ∉ (nothing, Any, Union{}) && return shortstr(inf)
  catch err
    # @error err
  end
  return ""
end

using Markdown

completiondescription(docs) = ""
completiondescription(docs::Markdown.MD) = begin
  md = CodeTools.flatten(docs).content
  for part in md
    if part isa Markdown.Paragraph
      desc = Markdown.plain(part)
      occursin("No documentation found.", desc) && return ""
      return strlimit(desc, DESCRIPTION_LIMIT)
    end
  end
  return ""
end


text:
import LinearAlgebra: norm

@doc doc"""
    Rotations{N} <: Manifold

Special orthogonal manifold $\mathrm{SO}(n)$ represented by $n\times n$
real-valued orthogonal matrices with determinant $+1$.

# Constructor

    Rotations(n)

generates the $\mathrm{SO}(n) \subset \mathbb R^{n\times n}$
"""
struct Rotations{N} <: Manifold end
Rotations(n::Int) = Rotations{n}()

function representation_size(::Rotations{N}) where N
    return (N, N)
end

@doc doc"""
    manifold_dimension(M::Rotations)

returns the dimension of the manifold $\mathrm{SO}(n)$, i.e. $\frac{n(n-1)}{2}$.
"""
manifold_dimension(M::Rotations{N}) where {N} = div(N*(N-1), 2)

@doc doc"""
    inner(M, x, w, v)

compute the inner product of the two tangent vectors `w, v` from the tangent
plane at `x` on the special orthogonal space `M=`$\mathrm{SO}(n)$ using the
restriction of the metric from the embedding, i.e.

$(v, w)_x = \operatorname{tr}(v^T w)$.

Tangent vectors are represented by matrices.
"""
inner(M::Rotations, x, w, v) = dot(w, v)

norm(M::Rotations, x, v) = norm(v)

project_tangent!(M::Rotations, w, x, v) = w .= (v .- transpose(v))./2

function flat!(M::Rotations, v::FVector{CotangentSpaceType}, x, w::FVector{TangentSpaceType})
    copyto!(v.data, w.data)
    return v
end

function sharp!(M::Rotations, v::FVector{TangentSpaceType}, x, w::FVector{CotangentSpaceType})
    copyto!(v.data, w.data)
    return v
end

function exp!(M::Rotations, y, x, v)
    y .= x * exp(v)
    return y
end

function exp!(M::Rotations{2}, y, x, v)
    θ = vee(M, x, v)[1]
    @assert size(y) == (2, 2)
    @assert size(x) == (2, 2)
    @inbounds begin
        sinθ, cosθ = sincos(θ)
        y[1] = x[1] * cosθ + x[3] * sinθ
        y[2] = x[2] * cosθ + x[4] * sinθ
        y[3] = x[3] * cosθ - x[1] * sinθ
        y[4] = x[4] * cosθ - x[2] * sinθ
    end
    return y
end

function exp!(M::Rotations{3}, y, x, v)
    θ = norm(M, x, v) / sqrt(2)
    if θ ≈ 0
        a = 1 - θ^2 / 6
        b = θ / 2
    else
        a = sin(θ) / θ
        b = (1 - cos(θ)) / θ^2
    end
    y .= x .+ x * (a .* v .+ b .* (v^2))
    return y
end

@doc doc"""
    angles_4d_skew_sym_matrix(A)

The Lie algebra of $\mathrm{SO}(4)$ consists of 4x4 skew-symmetric matrices.
The unique imaginary components of their eigenvalues are the angles of the two
plane rotations. This function computes these more efficiently than
`eigvals`.

By convention, the returned values are sorted in decreasing order
(corresponding to the same ordering of _angles_ as
[`cos_angles_4d_rotation_matrix`](@ref)).
"""
function angles_4d_skew_sym_matrix(A)
    @assert size(A) == (4, 4)
    @inbounds begin
        halfb = (A[1,2]^2 + A[1,3]^2 + A[2,3]^2 + A[1,4]^2 + A[2,4]^2 + A[3,4]^2) / 2
        c = (A[1,2] * A[3,4] - A[1,3] * A[2,4] + A[1,4] * A[2,3])^2
    end
    sqrtdisc = sqrt(halfb^2 - c)
    return sqrt(halfb + sqrtdisc), sqrt(halfb - sqrtdisc)
end

@doc doc"""
    exp!(M::Rotations{4}, y, x, v)

Exponential map of tangent vector `v` at point `x` from $\mathrm{SO}(4)$
manifold `M`. Result is saved to `y`.

The algorithm used is a more numerically stable form of those proposed in
[[Gallier, 2003](#Gallier2003)] and [[Andrica, 2013](#Andrica2013)].
"""
function exp!(M::Rotations{4}, y, x, v)
    T = eltype(v)
    α, β = angles_4d_skew_sym_matrix(v)
    sinα, cosα = sincos(α)
    sinβ, cosβ = sincos(β)
    α² = α^2
    β² = β^2
    Δ = β² - α²
    if !isapprox(Δ, 0; atol=1e-6)  # Case α > β ≥ 0
        sincα = sinα / α
        sincβ = β == 0 ? one(T) : sinβ / β
        a₀ = (β² * cosα - α² * cosβ) / Δ
        a₁ = (β² * sincα - α² * sincβ) / Δ
        a₂ = (cosα - cosβ) / Δ
        a₃ = (sincα - sincβ) / Δ
    elseif α == 0 # Case α = β = 0
        a₀ = one(T)
        a₁ = one(T)
        a₂ = T(1/2)
        a₃ = T(1/6)
    else  # Case α ⪆ β ≥ 0, α ≠ 0
        sincα = sinα / α
        r = β / α
        c = 1 / (1 + r)
        d = α * (α - β) / 2
        if α < 1e-2
            e = @evalpoly(α², T(1/3), T(-1/30), T(1/840), T(-1/45360))
        else
            e = (sincα - cosα) / α²
        end
        a₀ = (α * sinα + (1 + r - d) * cosα) * c
        a₁ = ((3 - d) * sincα - (2 - r) * cosα) * c
        a₂ = (sincα - (1 - r) / 2 * cosα) * c
        a₃ = (e + (1 - r) * (e - sincα / 2)) * c
    end

    v² = v * v
    y .= a₀ .* x .+ x * (a₁ .* v .+ a₂ .* v² .+ a₃ .* (v² * v))
    return y
end

@doc doc"""
    log!(M::Rotations, v, x, y)

compute the logarithmic map on the [`Rotations`](@ref) manifold
`M`$=\mathrm{SO}(n)$, which is given by
```math
\log_{x} y =
  \frac{1}{2} \bigl(\operatorname{Log}(x^{\mathrm{T}}y)
  - (\operatorname{Log} x^{\mathrm{T}}y)^{\mathrm{T}}),
```
where $\operatorname{Log}$ denotes the matrix logarithm, and save the result
to `v`.

For antipodal rotations the function returns one of the tangent vectors that
point at `y`.
"""
function log!(M::Rotations, v, x, y)
    U = transpose(x) * y
    v .= real(log_safe(U))
    project_tangent!(M, v, x, v)
    return v
end

function log!(M::Rotations{2}, v, x, y)
    U = transpose(x) * y
    @assert size(U) == (2, 2)
    @inbounds θ = atan(U[2], U[1])
    hat!(M, v, x, θ)
    return v
end

function log!(M::Rotations{3}, v, x, y)
    U = transpose(x) * y
    cosθ = (tr(U) - 1) / 2
    if cosθ ≈ -1
        eig = eigen_safe(U)
        ival = findfirst(λ -> isapprox(λ, 1), eig.values)
        vi = SVector{3}(1:3)
        ax = eig.vectors[vi, ival]
        hat!(M, v, x, π * ax)
    else
        v .= (U .- transpose(U)) ./ (2 * usinc_from_cos(cosθ))
    end
    return v
end

@doc doc"""
    cos_angles_4d_rotation_matrix(R)

4D rotations can be described by two orthogonal planes that are unchanged by
the action of the rotation (vectors within a plane rotate only within the
plane). The cosines of the two angles of rotation about these planes may be
obtained from the distinct real parts of the eigenvalues of the rotation
matrix. This function computes these more efficiently by solving the system

```math
\begin{aligned}
\cos\alpha + \cos\beta &= \frac{1}{2} \operatorname{tr}(R)\\
\cos\alpha + \cos\beta &= \frac{1}{8} \operatorname{tr}(R)^2
                           - \frac{1}{16} \operatorname{tr}((R - R^T)^2) - 1.
\end{aligned}
```

By convention, the returned values are sorted in increasing order. See
[`angles_4d_skew_sym_matrix`](@ref). For derivation of the above, see
[[Gallier, 2013]](#Gallier2003).
"""
function cos_angles_4d_rotation_matrix(R)
    trR = tr(R)
    a = trR / 4
    b = sqrt(clamp(tr((R .- transpose(R))^2) / 16 - a^2 + 1, 0, Inf))
    return (a + b, a - b)
end

function log!(M::Rotations{4}, v, x, y)
    U = transpose(x) * y
    cosα, cosβ = cos_angles_4d_rotation_matrix(U)
    α = acos(clamp(cosα, -1, 1))
    β = acos(clamp(cosβ, -1, 1))
    if α ≈ π && β ≈ 0
        A² = Symmetric((U - I) ./ 2)
        P = eigvecs(A²)
        E = similar(U)
        fill!(E, 0)
        α = acos(clamp(cosα, -1, 1))
        @inbounds begin
            E[2, 1] = -α
            E[1, 2] = α
        end
        v .= P * E * transpose(P)
    else
        v .= real(log_safe(U))
    end
    project_tangent!(M, v, x, v)
    return v
end

injectivity_radius(M::Rotations, x) = injectivity_radius(M)

injectivity_radius(M::Rotations) = π*sqrt(2.0)

@doc doc"""
    PolarRetraction

Retraction on the rotations manifold using the polar method.
This SVD-based retraction is a second-order approximation of the
exponential map. Let

$USV = x + xv$

be the singular value decomposition, then the formula reads

$\operatorname{retr}_x v = UV^\mathrm{T}$

Retraction is performed by the function [`retract!(::Rotations, y, x, v, ::PolarRetraction)`](@ref)
"""
struct PolarRetraction <: AbstractRetractionMethod end

@doc doc"""
    retract_polar!(M::Rotations, y, x, v, method::PolarRetraction)

Compute the SVD-based retraction [`PolarRetraction`](@ref), a second-order
approximation of the exponential map.
"""
function retract!(M::Rotations, y, x, v, method::PolarRetraction)
    A = x + x*v
    S = svd(A)
    y .= S.U * transpose(S.V)
    return y
end

"""
    QRRetraction

Retraction on the rotations manifold using the QR method, a first order
approximation of the exponential map.
"""
struct QRRetraction <: AbstractRetractionMethod end

@doc doc"""
    retract!(M, y, x, v, method::QRRetraction)

This QR-based retraction is a first-order approximation of the exponential map.
"""
function retract!(M::Rotations, y::AbstractArray{T}, x, v, method::QRRetraction) where T
    A = x + x*v
    qr_decomp = qr(A)
    d = diag(qr_decomp.R)
    D = Diagonal( sign.(d .+ convert(T, 0.5)) )
    y .= qr_decomp.Q * D
    return y
end

retract!(M::Rotations, y, x, v) = retract!(M, y, x, v, QRRetraction())

"""
    PolarInverseRetraction

Inverse retraction on the rotations manifold using the polar method.
"""
struct PolarInverseRetraction <: AbstractInverseRetractionMethod end

injectivity_radius(::Rotations, x, ::PolarRetraction) = π*sqrt(2.0)/2

@doc doc"""
    inverse_retract!(M, v, x, y, ::PolarInverseRetraction)

Compute a vector from the tagent space $T_x\mathrm{SO}(n)$
of the point `x` on the [`Rotations`](@ref) manifold `M`
with which the point `y` can be reached by the
[`PolarRetraction`](@ref) from the point `x` after time 1.

The formula reads
$v = -\frac{1}{2}(x^{\mathrm{T}}ys - (x^{\mathrm{T}}ys)^{\mathrm{T}})$

where $s$ is the solution to the Sylvester equation

$x^{\mathrm{T}}ys + s(x^{\mathrm{T}}y)^{\mathrm{T}} + 2\mathrm{I}_n = 0.$
"""
function inverse_retract!(M::Rotations, v, x, y, method::PolarInverseRetraction)
    A = transpose(x) * y
    H = 2 * one(x)
    try
        B = sylvester(collect(A), collect(transpose(A)), collect(H))
        C = A * B
        v .= (transpose(C) .- C)./2
    catch e
        if isa(e, LinearAlgebra.LAPACKException)
            throw(OutOfInjectivityRadiusError())
        else
            rethrow()
        end
    end
    return v
end

"""
    QRInverseRetraction

Inverse retraction on the rotations manifold using the QR method.
"""
struct QRInverseRetraction <: AbstractInverseRetractionMethod end

@doc doc"""
    inverse_retract!(M::Rotations, x, y, ::QRInverseRetraction)

Compute a vector from the tagent space $T_x\mathrm{SO}(n)$
of the point `x` on the [`Rotations`](@ref) manifold `M`
with which the point `y` can be reached by the
[`QRRetraction`](@ref) from the point `x` after time 1.
"""
function inverse_retract!(M::Rotations{N}, v, x, y, ::QRInverseRetraction) where N
    A = transpose(x) * y
    R = zero(v)
    for i = 1:N
        b = zeros(i)
        b[end] = 1
        b[1:(end-1)] = - transpose(R[1:(i-1), 1:(i-1)]) * A[i, 1:(i-1)]
        R[1:i, i] = A[1:i, 1:i] \ b
    end
    C =  A * R
    v .= (C .- transpose(C))./2
    return v
end

zero_tangent_vector(M::Rotations, x) = zero(x)
zero_tangent_vector!(M::Rotations, v, x) = (v .= zero(x))

function hat!(M::Rotations{2}, Ω, x, θ::Real)
    @assert length(Ω) == 4
    @inbounds begin
        Ω[1] = 0;  Ω[3] = -θ
        Ω[2] = θ;  Ω[4] = 0
    end
    return Ω
end

hat!(M::Rotations{2}, Ω, x, ω) = hat!(M, Ω, x, ω[1])

@doc doc"""
    hat!(M::Rotations, Ω, x, ω)

Convert the unique tangent vector components $\omega$ at point $x$ on rotations
group $\mathrm{SO}(n)$ to the matrix representation $\Omega$ of the tangent
vector. See [`vee!`](@ref) for the conventions used.
"""
function hat!(M::Rotations{N}, Ω, x, ω) where {N}
    @assert size(Ω) == (N, N)
    @assert length(ω) == manifold_dimension(M)
    @inbounds begin
        Ω[1,1] = 0;      Ω[1,2] = -ω[3];  Ω[1,3] = ω[2]
        Ω[2,1] = ω[3];   Ω[2,2] = 0;      Ω[2,3] = -ω[1]
        Ω[3,1] = -ω[2];  Ω[3,2] = ω[1];   Ω[3,3] = 0
        k = 4
        for i=4:N
            for j=1:i-1
                Ω[i,j] = ω[k]
                Ω[j,i] = -ω[k]
                k += 1
            end
            Ω[i,i] = 0
        end
    end
    return Ω
end

vee!(M::Rotations{2}, ω, x, Ω) = (ω[1] = Ω[2])

@doc doc"""
    vee!(M::Rotations, ω, x, Ω)

Extract the unique tangent vector components $\omega$ at point $x$ on rotations
group $\mathrm{SO}(n)$ from the matrix representation $\Omega$ of the tangent
vector.

The basis on the Lie algebra $\mathfrak{so}(n)$ is chosen such that for
$\mathrm{SO}(2)$, $\omega=\theta=\Omega_{21}$ is the angle of rotation, and
for $\mathrm{SO}(3)$,
$\omega = (\Omega_{32}, \Omega_{13}, \Omega_{21}) = \theta u$ is the
angular velocity and axis-angle representation, where $u$ is the unit vector
along the axis of rotation.

For $\mathrm{SO}(n)$ where $n \ge 4$, the additional elements of $\omega$ are
$\omega_{i (i - 3)/2 + j + 1} = \Omega_{ij}$, for $i \in [4, n], j \in [1,i)$.
"""
function vee!(M::Rotations{N}, ω, x, Ω) where {N}
    @assert size(Ω) == (N, N)
    @assert length(ω) == manifold_dimension(M)
    @inbounds begin
        ω[1] = Ω[3,2]
        ω[2] = Ω[1,3]
        ω[3] = Ω[2,1]

        k = 4
        for i=4:N
            for j=1:i-1
                ω[k] = Ω[i,j]
                k += 1
            end
        end
    end
    return ω
end

"""
    check_manifold_point(M,x; kwargs...)

checks, whether `x` is a valid point on the [`Rotations`](@ref) `M`,
i.e. is an array of size [`manifold_dimension`](@ref)`(M)` and represents a
valid rotation.
The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_manifold_point(M::Rotations{N},x; kwargs...) where {N}
    if size(x) != (N, N)
        return DomainError(size(x),
            "The point $(x) does not lie on $M, since its size is not $((N, N)).")
    end
    if !isapprox(det(x), 1; kwargs...)
        return DomainError(det(x), "The determinant of $x has to be +1 but it is $(det(x))")
    end
    if !isapprox(transpose(x)*x, one(x); kwargs...)
        return DomainError(norm(x), "$x has to be orthogonal but it's not")
    end
    return nothing
end

"""
    check_tangent_vector(M,x,v; kwargs... )

checks whether `v` is a tangent vector to `x` on the [`Rotations`](@ref)
space `M`, i.e. after [`check_manifold_point`](@ref)`(M,x)`, `v` has to be of same
dimension as `x` and orthogonal to `x`.
The tolerance for the last test can be set using the ´kwargs...`.
"""
function check_tangent_vector(M::Rotations{N},x,v; kwargs...) where N
    perr = check_manifold_point(M,x)
    perr === nothing || return perr
    if size(v) != (N, N)
        return DomainError(size(v),
            "The array $(v) is not a tangent to a point on $M since its size does not match $((N, N)).")
    end
    if !isapprox(transpose(v)+v, zero(v); kwargs...)
        return DomainError(size(v),
            "The array $(v) is not a tangent to a point on $M since it is not skew-symmetric.")
    end
    return nothing
end

"""
    normal_tvector_distribution(M::Rotations, x, σ)

Normal distribution in ambient space with standard deviation `σ`
projected to tangent space at `x`.
"""
function normal_tvector_distribution(M::Rotations, x, σ)
    d = Distributions.MvNormal(reshape(zero(x), :), σ)
    return ProjectedFVectorDistribution(TangentBundleFibers(M), x, d, project_vector!, x)
end

"""
    NormalRotationDistribution(M::Rotations, d::Distribution, x::TResult)

Distribution that returns a random point on the manifold [`Rotations`](@ref)
`M`. Random point is generated using base distribution `d` and the type
of the result is adjusted to `TResult`.

See [`normal_rotation_distribution`](@ref) for details.
"""
struct NormalRotationDistribution{TResult, TM<:Rotations, TD<:Distribution} <: MPointDistribution{TM}
    manifold::TM
    distr::TD
end

function NormalRotationDistribution(M::Rotations, d::Distribution, x::TResult) where TResult
    return NormalRotationDistribution{TResult, typeof(M), typeof(d)}(M, d)
end

function _fix_random_rotation(A::AbstractMatrix)
    s=diag(sign.(qr(A).R))
    D=Diagonal(s)
    C = qr(A).Q*D
    if det(C) < 0
        C[:,[1,2]] = C[:,[2,1]]
    end
    return C
end

function rand(rng::AbstractRNG, d::NormalRotationDistribution{TResult,Rotations{N}}) where {TResult,N}
    if N==1
        return convert(TResult, ones(1,1))
    else
        A = reshape(rand(rng, d.distr), (N, N))
        return convert(TResult, _fix_random_rotation(A))
    end
end

function _rand!(rng::AbstractRNG, d::NormalRotationDistribution{TResult,Rotations{N}}, x::AbstractArray{<:Real}) where {TResult,N}
    x .= rand(rng, d)
    return x
end

@doc doc"""
    normal_rotation_distribution(M::Rotations, x, σ::Real)

return a random point on the manifold [`Rotations`](@ref) `M`
by generating a (Gaussian) random orthogonal matrix with determinant $+1$. Let

$QR = A$

be the QR decomposition of a random matrix $A$, then the formula reads

$x = QD$

where $D$ is a diagonal matrix with the signs of the diagonal entries of $R$,
i.e.

$D_{ij}=\begin{cases} \operatorname{sgn}(R_{ij}) & \text{if} \; i=j \\ 0 & \, \text{otherwise} \end{cases}.$

It can happen that the matrix gets -1 as a determinant. In this case, the first
and second columns are swapped.

The argument `x` is used to determine the type of returned points.
"""
function normal_rotation_distribution(M::Rotations{N}, x, σ::Real) where N
    d = Distributions.MvNormal(zeros(N*N), σ)
    return NormalRotationDistribution(M, d, x)
end


text:
using MemoryInspector
using Test

@testset "MemoryInspector.jl" begin
    include("MemoryInspector.jl")
end


text:
module Timestamps

import Dates
import Dates: Date,
              Time,
              TimeType,
              DateTime,
              Hour,
              Minute,
              Second,
              Nanosecond,
              UTInstant,
              totaldays,
              argerror,
              daysinmonth,
              days,
              value,
              year, month, day, hour, minute, second, millisecond, microsecond, nanosecond

export Timestamp, TimestampTZ, TimeZoneOffset

"Nanosecond precision Date and Time."
abstract type AbstractTimestamp <: TimeType end

struct Timestamp <: AbstractTimestamp
    date::Date
    time::Time
end

# https://github.com/oracle/odpi/issues/38
# When fetched, tzMinuteOffset is negative if tzHourOffset is negative
# and is positive when tzHourOffset is positive.
struct TimeZoneOffset
    hour::Int8
    minute::Int8

    function TimeZoneOffset(hour::Int8, minute::Int8)
        @assert minute == 0 || signbit(hour) == signbit(minute) "`hour` and `minute` Time Zone offset must have the same sign."
        return new(hour, minute)
    end
end

TimeZoneOffset(hour::Integer, minute::Integer) = TimeZoneOffset(Int8(hour), Int8(minute))

struct TimestampTZ{L} <: AbstractTimestamp
    ts::Timestamp
    tz_offset::TimeZoneOffset

    function TimestampTZ{L}(ts::Timestamp, tz_offset::TimeZoneOffset) where {L}
        @assert isa(L, Bool)
        return new{L}(ts, tz_offset)
    end
end

TimestampTZ(local_time_zone::Bool, ts::Timestamp, tz_offset::TimeZoneOffset) = TimestampTZ{local_time_zone}(ts, tz_offset)

ts_date(t::Timestamp) = t.date
ts_time(t::Timestamp) = t.time

ts_date(t::TimestampTZ) = ts_date(t.ts)
ts_time(t::TimestampTZ) = ts_time(t.ts)

function is_ltz(t::TimestampTZ{L}) :: Bool where {L}
    return L
end

"""
    Timestamp(y, [m, d, h, mi, s, ns]) -> Timestamp
Construct a `Timestamp` type by parts. Arguments must be convertible to [`Int64`](@ref).
"""
function Timestamp(y::Integer, m::Integer=1, d::Integer=1,
                   h::Integer=0, mi::Integer=0, s::Integer=0, ns::Integer=0)
    date = Date(y,m,d)
    time = Time(Nanosecond( ns + 1E9 * ( s + 60*mi + 3600*h ) ))
    return Timestamp(date, time)
end

function TimestampTZ(local_time_zone::Bool, y::Integer, m::Integer=1, d::Integer=1,
                     h::Integer=0, mi::Integer=0, s::Integer=0, ns::Integer=0,
                     tz_hour_offset::Integer=0, tz_minute_offset::Integer=0)
    ts = Timestamp(y, m, d, h, mi, s, ns)
    return TimestampTZ(local_time_zone, ts, TimeZoneOffset(tz_hour_offset, tz_minute_offset))
end

Base.eps(t::AbstractTimestamp) = Nanosecond(1)
Base.:(==)(t1::Timestamp, t2::Timestamp) = t1.date == t2.date && t1.time == t2.time
Base.:(==)(o1::TimeZoneOffset, o2::TimeZoneOffset) = o1.hour == o2.hour && o1.minute == o2.minute
Base.:(==)(t1::TimestampTZ, t2::TimestampTZ) = t1.ts == t2.ts && t1.tz_offset == t2.tz_offset
Base.hash(ts::Timestamp) = 1 + hash(ts.date) + hash(ts.time)
Base.hash(ts::TimestampTZ) = 2 + hash(ts.ts) + hash(ts.tz_offset)
Base.hash(o::TimeZoneOffset) = 4 + hash(o.hour) + hash(o.minute)

function Base.:(==)(ts::Timestamp, dt::DateTime)
    if microsecond(ts) == 0 && nanosecond(ts) == 0
        return (
               days(ts) == days(dt)
            && hour(ts) == hour(dt)
            && minute(ts) == minute(dt)
            && second(ts) == second(dt)
            && millisecond(ts) == millisecond(dt)
        )
    else
        # Timestamp has nanosecond precision. DateTime has millisecond precision.
        return false
    end
end

Base.:(==)(dt::DateTime, ts::Timestamp) = ts == dt

Dates.days(t::AbstractTimestamp) = Dates.days(ts_date(t))
nanoseconds(t::AbstractTimestamp) = Dates.value(ts_time(t))

for fun in (:hour, :minute, :second, :millisecond, :microsecond, :nanosecond)
    @eval begin
        ($fun)(t::AbstractTimestamp) = ($fun)(ts_time(t))
    end
end

function Timestamp(dt::DateTime)
    return Timestamp(year(dt), month(dt), day(dt), hour(dt), minute(dt), second(dt), Int64(millisecond(dt)*1E6))
end

function Base.convert(::Type{DateTime}, ts::Timestamp)
    return DateTime(ts.date, ts.time)
end

function Base.convert(::Type{Date}, ts::Timestamp)
    return convert(Date, convert(DateTime, ts))
end

end # module Timestamps


text:
"""
    oddsratio(a, b, c, d; alpha=.05, exact=false)
    oddsratio(A, B; alpha=.05, exact=false)	# todo

Calculate odds ratio of 2x2 table (Case-Control Studies).
If A, B are vectors, a two-way table is generated using `crossTable` and OR is calculated.

Exact method is not implemented.

|             | case | non-cases |         |
|:------------|:-----|:----------|:--------|
| exposed     | a    | b         | a+b     |
| non-exposed | c    | d         | c+d     |
|             | a+c  | b+d       | a+b+c+d |

## Example outputs

```julia
B = 2×2 Named Array{Int64,2}
Exposed ╲ Case │  +   -
───────────────┼───────
+              │ 25  12
-              │  9  69
1×3 Named Array{Float64,2}
      ╲  │      OR    lower    upper
─────────┼──────────────────────────
OR 95%CI │ 15.9722  6.00717  42.4679
```
"""
function oddsratio(a, b, c, d; alpha=0.05, exact=false)
		if any(i==0 for i in [a b c d]) @error("Zero cell(s), cannot compute OR's variance") end
		B = NamedArray([a b; c d], ( ["+", "-"], ["+", "-"] ), ("Exposed", "Case"))
		rowname = string("OR ", string(100*(1-alpha)), "%CI")
		A = NamedArray([.0 .0 .0], ( [rowname], ["OR", "lower", "upper"] ), ("", ""))
		A[1] = (a * d) / (b * c)
		v = sum([1/i for i in [a, b, c, d]])
		z = quantile(Normal(), 1-alpha/2)
		A[2], A[3] = A[1] * [exp(-z*sqrt(v)) exp(z*sqrt(v))]
		@show B
		A
end

"""
    riskratio(a, b, c, d; alpha=.05)
    riskratio(A, B; alpha=.05, exact=false)	# todo

Calculate risk ratio (RR) of 2x2 table; this can also be used for calculating prevalence ratio (PR).

If A, B are vectors, a two-way table is generated using `crossTable` and RR is calculated.

|             | case | non-cases |         |
|:------------|:-----|:----------|:--------|
| exposed     | a    | b         | a+b     |
| non-exposed | c    | d         | c+d     |
|             | a+c  | b+d       | a+b+c+d |

## Example outputs

```julia
julia> riskratio(122,355,22,965, alpha=0.01)
B = 2×2 Named Array{Int64,2}
Exposed ╲ Case │   +    -
───────────────┼─────────
+              │ 122  355
-              │  22  965
1×3 Named Array{Float64,2}
        ╲  │      RR    lower    upper
───────────┼──────────────────────────
RR 99.0%CI │ 11.4746  6.43047  20.4752
```
"""
function riskratio(a, b, c, d; alpha=0.05)
		if any(i==0 for i in [a c]) error("Zero cell(s), cannot compute RR's variance") end
		B = NamedArray([a b; c d], ( ["+", "-"], ["+", "-"] ), ("Exposed", "Case"))
		rowname = string("RR ", string(100*(1-alpha)), "%CI")
		A = NamedArray([.0 .0 .0], ( [rowname], ["RR", "lower", "upper"] ), ("", ""))
		p11, p12 = (a/(a+b), c/(c+d))
		A[1] = p11 / p12
		v = (1-p11)/a + (1-p12)/c
		z = quantile(Normal(), 1-alpha/2)
		A[2], A[3] = A[1] * [exp(-z*sqrt(v)) exp(z*sqrt(v))]
		@show B
		A
end


text:
abstract type AbstractServiceFormulation end

abstract type AbstractReservesFormulation<:AbstractServiceFormulation end

abstract type AbstractRegulationReserveFormulation<:AbstractReservesFormulation end

struct RampLimitedReserve<:AbstractReservesFormulation end

struct LoadProportionalReserve<:AbstractReservesFormulation end


text:
if VERSION < v"1.0.0"
    error("VS Code julia language server only works with julia 1.0.0+")
end

@info "Starting the Julia Language Server"

using InteractiveUtils

include("error_handler.jl")

try
    if length(Base.ARGS) != 4
        error("Invalid number of arguments passed to julia language server.")
    end

    conn = stdout
    (outRead, outWrite) = redirect_stdout()

    if Base.ARGS[2] == "--debug=no"
        const global ls_debug_mode = false
    elseif Base.ARGS[2] == "--debug=yes"
        const global ls_debug_mode = true
    end

    using LanguageServer, Sockets, SymbolServer

    server = LanguageServerInstance(stdin, conn, ls_debug_mode, Base.ARGS[1], Base.ARGS[4], global_err_handler)
    run(server)
catch e
    global_err_handler(e, catch_backtrace())
end

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode