What is the functional derivative? I do not understand, if the functional derivative is


*

*a function

*a generalized function (distribution)

*a functional itself

*something different (see Euler-Lagrange)


To clarify my question, I have seen multiple instances of functional derivative definitions
Functionals
When the Functional gets Taylor expanded (here using a "good" $\eta(x)$) we get
$$F[y(x)+\epsilon \eta(x)] = F[y(x)] + \frac{dF[y(x) + \epsilon \eta(x)]}{d\epsilon}\Big|_{\epsilon=0}\cdot \epsilon + ...$$
as I understood, the term on the RHS is the functional derivative. But since the LHS is a functional and the RHS is a functional + a real number ($\epsilon$) times the functional derivative, I conclude that the functional derivative must also be a functional.
Functions/Distributions
The english wikipedia page [2] states, that the functional derivative is defined as
$$\int{\frac{\delta F}{\delta \rho} (x)\phi(x)dx}=\frac{dF[\rho(x) + \epsilon \phi(x)]}{d\epsilon}\Big|_{\epsilon=0}$$
notice that the RHS is equivalent to the functional derivative defined above. However, it is $$\frac{\delta F}{\delta \rho} (x)$$ that is defined to be the functional derivative, and not the RHS (as I concluded above). Therefore I may as well assume that the functional derivative is a function/distribution.
Something else
The solution to the Euler-Lagrange Equation (one dimensional for simplicity) given an Energy Functional $J[y] = \int_{a}^{b}{L(x,y,y')}$ is
$$\frac{\delta J}{\delta y} = \frac{dL}{dy} - \frac{d}{dx}(\frac{dL}{dy'}) = 0$$
here, $\frac{\delta J}{\delta y}$ is supposedly the fractional derivative of the integral, which has to be stationary. RHS tells me that the functiona derivative is a differential equation - which has a function as a solution - but I am now completely unsure what the functional derivative in itself actualy is. 
I have seen multiple viewpoints, each and every one cluttering my intuition even more. For instance the wikipedia article claims that $\frac{\delta F}{\delta \rho} (x)$ has to be seen as a "gradient" (which is a vector in multivariate calculus), while $\int{\frac{\delta F}{\delta \rho} (x)\phi(x)dx}$ has to be thought of like a directional derivative (which is the inner product of the gradient and the direction vector). But since there are no bounds on the integral the "directional derivative" is also a function, or am I mistaken?
[1] http://lab.sentef.org/wp-content/uploads/2016/11/Tutorial_02.pdf page 4
[2] https://en.wikipedia.org/wiki/Functional_derivative
 A: The expression
$\delta F[\rho,\phi] := \frac{dF[\rho(x) + \epsilon \phi(x)]}{d\epsilon}\Big|_{\epsilon=0},$
when defined, is a functional of $\rho$ and $\phi.$ The dependency on $\rho$ is usually non-linear, while the dependency on $\phi$ is usually linear.
If the expression is restricted to $\phi \in C_c^\infty(\mathbb R^n)$ and the dependency on $\phi$ is linear, then the mapping $\phi \mapsto \delta F[\rho,\phi]$ is usually a distribution. Often this distribution can be identified with a function.
Thus, $\delta F[\rho,\phi]$ is a functional, usually a distribution, and often a function.
Often we have $F[\rho] = \int L(x, \rho(x), \rho'(x)) \, dx$ for some Lagrangian $L.$ Then, if $\phi$ vanishes on the boundary of the domain, 
$$
\delta F[\rho,\phi] = \int \left( \frac{\partial L}{\partial \rho} \phi(x) + \frac{\partial L}{\partial \rho'} \phi'(x) \right) dx
= \int \left( \frac{\partial L}{\partial \rho} - \frac{d}{dx}\frac{\partial L}{\partial \rho'} \right) \phi(x) \, dx.
$$
In this case, $\delta F[\rho,\phi]$ is given by an integral of a function (the parenthesis) times $\phi.$ Thus this falls into the case "Often this distribution can be identified with a function".
A: Functionals of Smooth Fields
Suppose you have a smooth multivariable scalar field ${f}\in\mathcal{C}^{\infty}(\mathbb{R}^{n})$. Then a functional ${W}\in\mathrm{Funct}(\mathcal{C}^{\infty}(\mathbb{R}^{n}))$ is a mapping from $\mathcal{C}^{\infty}(\mathbb{R}^{n})$ to $\mathbb{R}$ that takes the form
\begin{align*}
{W}[f]=\int\mathcal{L}(x)[f(x)]{\,}\mathrm{d}^{n}{x}{\,}{,}
\end{align*}
where $\mathcal{L}(x)[f(x)]$ is the Lagrangian, that depends both on the coordinate variables ${x}$, as well as functionally depends on the field in coordinates ${f(x)}$ and its consecutive derivatives $\big\lbrace\frac{\partial{f(x)}}{\partial{x}^{\alpha}},\frac{\partial^{2}{f(x)}}{\partial{x}^{\alpha_{1}}\partial{x}^{\alpha_{2}}},\frac{\partial^{3}{f(x)}}{\partial{x}^{\alpha_{1}}\partial{x}^{\alpha_{2}}\partial{x}^{\alpha_{3}}},\cdots\big\rbrace$. Now, rather then looking at the variational derivative as a limit, try thinking of $\frac{\delta}{\delta{f(x)}}\in\mathrm{Der}(\mathrm{Funct}(\mathcal{C}^{\infty}(\mathbb{R}^{n})))$ as the derivative operator on $\mathrm{Funct}(\mathcal{C}^{\infty}(\mathbb{R}^{n}))$,
similar to how $\frac{\partial}{\partial{x}^{\alpha}}\in\mathrm{Der}(\mathcal{C}^{\infty}(\mathbb{R}^{n}))$ is the derivative operator on $\mathcal{C}^{\infty}(\mathbb{R}^{n})$, that satisfies the algebra
\begin{align*}
(\mathrm{I}){\,\,}&\frac{\delta({a}\mathcal{L}_{\alpha}(x)[f(x)]+{b}\mathcal{L}_{\beta}(x)[f(x)])}{\delta{f(y)}}={a}\frac{\delta\mathcal{L}_{\alpha}(x)[f(x)]}{\delta{f(y)}}+{b}\frac{\delta\mathcal{L}_{\beta}(x)[f(x)]}{\delta{f(y)}}{\quad}{a}{,}{b}\in\mathbb{R}{\,}{,}\\
(\mathrm{II}){\,\,}&\frac{\delta(\mathcal{L}_{\alpha}(x)[f(x)]{\,}\mathcal{L}_{\beta}(x)[f(x)])}{\delta{f(y)}}=\mathcal{L}_{\alpha}(x)[f(x)]\frac{\delta\mathcal{L}_{\beta}(x)[f(x)]}{\delta{f(y)}}+{\,}\mathcal{L}_{\beta}(x)[f(x)]\frac{\delta\mathcal{L}_{\alpha}(x)[f(x)]}{\delta{f(y)}}{\,}{,}\\
(\mathrm{III}){\,\,}&\frac{\delta\mathcal{L}(x)[g(f(x))]}{\delta{f(y)}}=\frac{\delta\mathcal{L}(x)[g(f(x))]}{\delta{g(f(y))}}\frac{\mathrm{d}{g(z)}}{\mathrm{d}{z}}\bigg\vert_{{z}={f(y)}}{\,}{.}
\end{align*}
The functional derivative can now be written as
\begin{align*}
\frac{\delta{W}[f]}{\delta{f(y)}}&=\int\frac{\delta\mathcal{L}(x)[f(x)]}{\delta{f(y)}}{\,}\mathrm{d}^{n}{x}=\int\sum_{{m}={0}}^{\infty}\frac{\partial\mathcal{L}(x)[f(x)]}{\partial\partial^{m}{f(x)}/\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\frac{\delta}{\delta{f(y)}}\bigg[\frac{\partial^{m}{f(x)}}{\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\bigg]{\,}\mathrm{d}^{n}{x}\\
&=\int\sum_{{m}={0}}^{\infty}\bigg(\frac{\partial^{m}}{\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\bigg[\frac{\partial\mathcal{L}(x)[f(x)]}{\partial\partial^{m}{f(x)}/\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\frac{\delta{f(x)}}{\delta{f(y)}}\bigg]\\
&+(-{1})^{m}\frac{\partial^{m}}{\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\bigg[\frac{\partial\mathcal{L}(x)[f(x)]}{\partial\partial^{m}{f(x)}/\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\bigg]\frac{\delta{f(x)}}{\delta{f(y)}}\bigg){\,}\mathrm{d}^{n}{x}{\,}{,}
\end{align*}
where $\frac{\delta{f(x)}}{\delta{f(y)}}=\Delta({x}-{y})$ is the Dirac delta function that is simply defined as
\begin{align*}
\int{f(x)}\Delta({x}-{y}){\,}\mathrm{d}^{n}{x}={f(y)}.
\end{align*}
Furthermore, since $\frac{\partial^{m}}{\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\big[\frac{\partial\mathcal{L}(x)[f(x)]}{\partial\partial^{m}{f(x)}/\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\Delta({x}-{y})\big]$ is a total derivative, it vanishes when being integrated over a set without boundary. And with all that combined, the functional derivative finally becomes
\begin{align*}
\frac{\delta{W}[f]}{\delta{f(y)}}&=\int\sum_{{m}={0}}^{\infty}(-{1})^{m}\frac{\partial^{m}}{\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\bigg[\frac{\partial\mathcal{L}(x)[f(x)]}{\partial\partial^{m}{f(x)}/\partial{x}^{\alpha_{1}}\cdots\partial{x}^{\alpha_{m}}}\bigg]\Delta({x}-{y}){\,}\mathrm{d}^{n}{x}\\
&=\sum_{{m}={0}}^{\infty}(-{1})^{m}\frac{\partial^{m}}{\partial{y}^{\alpha_{1}}\cdots\partial{y}^{\alpha_{m}}}\bigg[\frac{\partial\mathcal{L}(y)[f(y)]}{\partial\partial^{m}{f(y)}/\partial{y}^{\alpha_{1}}\cdots\partial{y}^{\alpha_{m}}}\bigg]{\,}{.}
\end{align*}
Euler Lagrange Equation
The Euler Lagrange equation of an action is given by simply setting the variational derivative of the functional equal to zero, being
\begin{align*}
{0}\equiv\frac{\delta{W}[f]}{\delta{f(y)}}=\sum_{{m}={0}}^{\infty}(-{1})^{m}\frac{\partial^{m}}{\partial{y}^{\alpha_{1}}\cdots\partial{y}^{\alpha_{m}}}\bigg[\frac{\partial\mathcal{L}(y)[f(y)]}{\partial\partial^{m}{f(y)}/\partial{y}^{\alpha_{1}}\cdots\partial{y}^{\alpha_{m}}}\bigg]{\,}{.}
\end{align*}
Application in classical mechanics
In classical mechanics, the action functional of a curve has the form
\begin{align*}
{W}[x]=\int\mathcal{L}(t)[x(t)]{\,}\mathrm{d}{t}=\int\bigg(\frac{m}{2}\bigg\Vert\frac{\mathrm{d}{x(t)}}{\mathrm{d}{t}}\bigg\Vert^{2}-{m}\Phi(x(t))\bigg){\,}\mathrm{d}{t}{\,}{.}
\end{align*}
With the standard Euler Lagrange equation
\begin{align*}
{0}\equiv\frac{\delta{W}[x]}{\delta{{x}^{\alpha}(t)}}=\frac{\partial\mathcal{L}(t)[x(t)]}{\partial{{x}^{\alpha}(t)}}-\frac{\mathrm{d}}{\mathrm{d}{t}}\bigg[\frac{\partial\mathcal{L}(t)[x(t)]}{\partial\mathrm{d}{{x}^{\alpha}(t)}/\mathrm{d}{t}}\bigg]=-{m}\frac{\partial\Phi(x(t))}{\partial{{x}^{\alpha}(t)}}-{m}\frac{\mathrm{d}^{2}{{x}^{\alpha}(t)}}{\mathrm{d}{t}^{2}}{\,}{.}
\end{align*}
