# Maxwell fish eye for polarized light

Mher Davtyan  Institute of Radiophysics and Electronics, Ashtarak-2, 0203,
Armenia<EMAIL_ADDRESS>Zhyrair Gevorkian  Institute of Radiophysics
and Electronics, Ashtarak-2, 0203, Armenia  Yerevan Physics Institute, 2
Alikhanian Brothers St., 0036, Yerevan, Armenia<EMAIL_ADDRESS>Armen
Nersessian  Yerevan Physics Institute, 2 Alikhanian Brothers St., 0036,
Yerevan, Armenia  Yerevan State University, 1 Alex Manoogian St., Yerevan,
0025, Armenia  Bogoliubov Laboratory of Theoretical Physics, Joint Institute
for Nuclear Research, Dubna, Russia<EMAIL_ADDRESS>
###### Abstract

We consider the propagation of polarized light in the medium with Maxwell fish
eye refraction index profile. We show that polarization violates the
additional symmetries of medium, so that ray trajectories no longer remain
closed. Then we suggest a modified, polarization-dependent Maxwell fish eye
refraction index which restores all symmetries of initial profile and yields
closed trajectories of polarized light. Explicit expressions for the
polarization dependent integrals of motion and the solutions of corresponding
ray trajectories are presented.

##  I Introduction

It is well-known that minimal action principle came in physics from the
geometric optics. Initially, it was invented for the description of the
propagation of light and is presently known as the Fermat principle

${\cal S}_{Fermat}=\frac{1}{\lambdabar_{0}}\int d{l},\qquad
dl:=n(\mathbf{r})|d\mathbf{r}/d\tau|d\tau$ (1)

where $n(\mathbf{r})$ is the refraction index, and $\lambdabar_{0}$ is the
wavelength in vacuum. This action could be interpreted as the action of the
system on the three-dimensional curved space equipped with the “optical
metrics” or the Fermat metrics (see. arnold ) of Euclidean signature

$dl^{2}=n^{2}(\mathbf{r})d\mathbf{r}d\mathbf{r}\;.$ (2)

Thus, the symmetries of the system which describe the propagation of light in
a particular medium are coming from the symmetries of the optical metrics of
that particular medium. In systems with maximal number of functionally
independent integrals of motion ($2N-1$ integrals for $N$-dimensional system),
all the trajectories of the system become closed. Closeness of the
trajectories makes such profiles highly relevant in the study of cloaking and
perfect imaging phenomena. The most well-known profile of this sort is the so-
called “Maxwell fish eye” profile which is defined by the metrics of (three-
dimensional) sphere or pseudosphere (under pseudosphere we mean the upper (or
lower) sheet of the two-sheet hyperboloid)

$n_{Mfe}(\mathbf{r})=\frac{n_{0}}{|1+\kappa\mathbf{r}^{2}|},\qquad\kappa=\pm\frac{1}{4r^{2}_{0}}$
(3)

where plus/minus sign in the expression for $\kappa$ corresponds to the
sphere/pseudosphere with the radius $r_{0}$, and $n_{0}>0$. Apart from
applications in cloaking and perfect imaging phenomena Pendry06 ; Leonhardt06
; leonhardconformal , Maxwell fish eye is a common profile in quantum optics
with single atoms and photons perczel2018 , optical resonators turks14 ,
discrete spectrum radiation discretespectrumradiation  etc. Moreover, there
are many experimental implementations of the Maxwell fish eye lenses
zhao2021bifunctional ; lu2020millimeter ; koala2020 .

In most of the listed studies, scalar wave approximation was used and
polarization of electromagnetic waves was not taken into account. In these
systems the ray trajectories belong to the plane which is orthogonal to the
angular momentum. Introduction of spin (polarization) turns this plane to the
constant angle proportional to spin, moreover, it breaks the non-rotational
symmetries of the optical systems with Maxwell fish eye profile, so that
photon trajectories no longer remain closed gevdav20 . Hence, taking into
account of spin violates the key property of Maxwell fish-eye profile which
make it relevant in cloaking and perfect imaging phenomena.

In the present paper we continue our study of a polarized light passing
through the Maxwell fish eye profile within the geometrical optics
approximation. The key point of our study is that we introduce a polarization-
dependent deformation of Maxwell fish-eye profile

$n^{s}_{Mfe}(\mathbf{r})=\frac{n_{Mfe}(\mathbf{r})}{2}\left(1+\sqrt{1-\frac{4\kappa
s^{2}\lambdabar^{2}_{0}}{n_{0}}\frac{1}{n_{Mfe}(\mathbf{r})}}\right),$ (4)

where $n_{Mfe}(\mathbf{r})$ is original Maxwell fish eye profile given by (3),
and $s$ is a polarization of light. For the linearly/circularly polarized we
respectively have $s=0$/$s=1$. Proposed deformation restores all the
symmetries of the optical Hamiltonian with Maxwell fish eye profile which were
broken after the inclusion of polarization. The proposed profile ensures
closeness of the trajectories for polarized photons and can be used for
cloaking and perfect imaging of polarized photons.

The paper is organized as follows. In Section 2, we describe the Hamiltonian
formulation of the geometric optical system given by the action (1). We also
present some other textbook fact on duality between Coulomb and free-particle
system on (pseudo)sphere which allows to relate the Maxwell fish eye and
Coulomb profiles and will be used in our further consideration.

In Section 3, we present the Hamiltonian formalism for the polarized light
propagating in optical medium, and propose the general scheme of the
deformation of isotropic refraction index profile which allows us to restore
the initial symmetries after the inclusion of polarization.

In Section 4, we use the proposed scheme for the construction of ”polarized
Maxwell fish eye” profile (4) which inherits all the symmetries of the
original Maxwell fish eye profile(3) when light polarization is taken into
account. We present the explicit expressions for the symmetry generators of
respective Hamiltonian system and find the expressions of the Casimirs of
their symmetry algebra.

In Section 5 the explicit expressions for the trajectories of polarized light
are presented. It is shown that these trajectories are no longer orthogonal to
the angular momentum but turn to a fixed angle relative to it. Despite
deviations from circles these trajectories remain closed.

Through the text we will use the notation $r:=|\mathbf{r}|$,
$\mathbf{r}:=(x_{1},x_{2},x_{3})$, $\mathbf{p}:=(p_{1},p_{2},p_{3})$,
$p:=|\mathbf{p}|$, and so on.

##  II Scalar waves

Due to reparametrization-invariance of the action (1), the Hamiltonian
constructed by the standard Legendre transformation is identically zero.
However, the constraint between momenta and coordinates appears there

$\Phi:=\frac{{\bf p}^{2}}{n^{2}({\bf r})}-\lambdabar^{-2}_{0}=0.$ (5)

Hence, in accordance with the Dirac’s constraint theory dirac  the respective
Hamiltonian system is defined by the canonical Poisson brackets

$\\{x_{i},p_{j}\\}=\delta_{ij},\quad\\{p_{i},p_{j}\\}=\\{x_{i},x_{j}\\}=0,$
(6)

and by the Hamiltonian

$\mathcal{H}_{0}=\alpha({\bf p},{\bf r})\Phi=\alpha({\bf p},{\bf
r})\left(\frac{{p}^{2}}{n^{2}({\bf r})}-\lambdabar^{-2}_{0}\right)\approx 0.$
(7)

Here $\alpha$ is the Lagrangian multiplier which could be an arbitrary
function of coordinates and momenta, and $i,j=1,2,3$. The notation “weak
zero”, $\mathcal{H}_{0}\approx 0$, means that when writing down the
Hamiltonian equations of motion, we should take into account the constraint
(5) only after the differentiation,

$\frac{df({\bf r},{\bf
p})}{d\tau}=\\{f,\mathcal{H}_{0}\\}=\\{f,\alpha\\}\Phi+\alpha\\{f,\Phi\\}\approx\alpha\\{f,\Phi\\}.$
(8)

The arbitrariness in the choice of the function $\alpha$ reflects the
reparametrization-invariance of (1). For the description of the equations of
motion in terms of arc-length of the original Euclidian space one should
choose (see, e.g. bliokh )

$\alpha=\frac{n^{2}({\bf r})}{{p}+\lambdabar^{-1}_{0}n({\bf
r})},\qquad\Rightarrow\quad\mathcal{H}_{\rm Opt}=p-\lambdabar^{-1}_{0}n({\bf
r}).$ (9)

With this choice, the equations of motion take the conventional form ko

$\frac{d{\bf{p}}}{dl}=\lambdabar^{-1}_{0}{\bf\nabla}n({\bf
r}),\qquad\frac{d{\bf{r}}}{dl}=\frac{\bf p}{{p}},$ (10)

where $dl:=\alpha({\bf r},{\bf p})d\tau$ is the element of arc-length. These
equations describe the motion of a wave package with center coordinate ${\bf
r}$ and momentum ${\bf p}$ in the medium with refraction index $n({\bf r})$.

Assume we have the Hamiltonian system given by the Poisson bracket (6) and by
the Hamiltonian

$H=\frac{{p}^{2}}{g({\bf r})}+V({\bf r}).$ (11)

In accordance with the Mopertuit principle, after fixing the energy surface
$H=E$, we can relate its trajectories with the optical Hamiltonian (7) with
the refraction index

$n({\bf r})=\lambdabar_{0}\sqrt{2g({\bf r})(E-V({\bf r}))}.$ (12)

Clearly, the optical Hamiltonian (7) (as well as the Hamiltonian (9)) with the
refraction index (12) inherits all symmetries and constants of motion of the
Hamiltonian (11).

Canonical transformations preserve the symmetries of the Hamiltonians and
their level surfaces. Hence, we are able to construct the physically non-
equivalent optical Hamiltonians (and refraction indices) with the same
symmetry algebra. The simplest illustration is the well-known relation between
the Coulomb Hamiltonian which defines the so-called Coulomb refraction index
profile and the free-particle Hamiltonian on the three-dimensional sphere,
which defines the ”Maxwell fish eye” refraction index (see e.g. perelomov ).
Firstly, we fix the energy surface of the Coulomb Hamiltonian and get the
respective refraction index

$H_{Coul}-E:=\frac{{p}^{2}}{2}-\frac{{\gamma}}{r}-E=0,\quad\Rightarrow\quad
n_{Coul}=\lambdabar_{0}\sqrt{2(E+\gamma/r)},\qquad{\rm where}\quad\gamma>0.$
(13)

The constants of motion of the Coulomb problem (and of the respective optical
Hamiltonian) are given by the rotational momentum and by the Runge-Lenz vector

$\mathbf{L}=\mathbf{r}\times\mathbf{p},\qquad{\bf A}={\bf
L}\times\mathbf{p}+\gamma\frac{{\bf r}}{{r}}$ (14)

which form the algebra

$\\{A_{i},A_{j}\\}=-2\varepsilon_{ijk}H_{Coul}L_{k},\quad\\{A_{i},L_{j}\\}=\varepsilon_{ijk}A_{k},\quad\\{L_{i},L_{j}\\}=\varepsilon_{ijk}L_{k}.$
(15)

Now, let us perform a simple canonical transformation,

$({\bf p},{\bf r})\to(-{\bf r},{\bf p}).$ (16)

As a result, the first equation in (13) reads

${r}^{2}-\frac{2{\gamma}}{p}-2E=0\quad\Rightarrow\quad
p-\frac{2\gamma}{{r}^{2}-2E}=0$ (17)

Interpreting the second equation as an optical Hamiltonian, we get the
refraction index profile known as the “Maxwell fish eye” (3) with the
parameters $\kappa,n_{0}$ defined as follows

$\kappa:=-\frac{1}{2E},\qquad\frac{n_{0}}{\lambdabar_{0}}:=2\epsilon\kappa\gamma,$
(18)

where $\epsilon=-{\rm sgn}(r^{2}+1/\kappa)$.

The integrals of motion (14) result in the symmetry generators of the optical
Hamiltonian with the Maxwell fish eye refraction index

$\mathbf{L}\to\mathbf{L},\quad\mathbf{A}\to\frac{\mathbf{T}}{2\kappa},\quad\mathbf{T}=\left(1-\kappa
r^{2}\right)\mathbf{p}+2\kappa(\mathbf{rp})\mathbf{r}=\left(2-\frac{n_{0}}{n_{Mfe}({\bf
r})}\right)\mathbf{p}+2\kappa(\mathbf{rp})\mathbf{r}.$ (19)

These integrals form the $so(4)$ algebra for $\kappa>0$, and $so(1.3)$ algebra
for $\kappa<0$:

$\\{L_{i},L_{j}\\}=\varepsilon_{ijk}L_{k},\quad\\{T_{i},L_{j}\\}=\varepsilon_{ijk}T_{k},\quad\\{T_{i},T_{j}\\}=4\kappa\varepsilon_{ijk}L_{k}\;.$
(20)

In the next Sections we will use the above described duality for the
construction of the Maxwell fish eye profile for polarized light.

##  III Inclusion of polarization

Let us briefly discuss the inclusion of polarization.

To this end we should add to the scalar Lagrangian $L_{0}={\bf
p\dot{r}}-p+\lambda_{0}^{-1}n$ the additional term
$L_{1}=-s\mathbf{A}(\mathbf{p})\dot{\mathbf{p}}$, where $s$ is spin of photon,
and $\mathbf{A}$ is the the vector-potential of “Berry monopole” (i.e. the
potential of the magnetic (Dirac) monopole located at the origin of momentum
space) bliokh

$\mathbf{F}:=\frac{\partial}{\partial\mathbf{p}}\times\mathbf{A}(\mathbf{p})=\frac{\mathbf{p}}{{p}^{3}}$
(21)

From the Hamiltonian viewpoint this means to preserve the form of the
Hamiltonian (7) and replace the canonical Poisson brackets (6) by the twisted
ones

$\\{x_{i},p_{j}\\}=\delta_{ij},\qquad\\{x_{i},x_{j}\\}=s\varepsilon_{ijk}F_{k}(\mathbf{p}),\qquad\\{x_{i},x_{j}\\}=0,$
(22)

where $i,j,k=1,2,3$, and $F_{k}$ are the components of Berry monopole (21). On
this phase space the rotation generators take the form

$\mathbf{J}=\mathbf{r}\times\mathbf{p}+s\frac{{\mathbf{p}}}{{p}}$ (23)

while the equations of motion read

$\frac{d\mathbf{p}}{dl}=\lambdabar_{0}^{-1}\mathbf{\nabla}n(\mathbf{r}),\qquad\frac{d\mathbf{r}}{dl}=\frac{\mathbf{p}}{p}-\frac{s}{\lambdabar_{0}}\mathbf{F}\times\mathbf{\nabla}n({\bf
r}),$ (24)

However, the above procedure, i.e. twisting the Poisson bracket with
preservation of the Hamiltonian, violates the non-kinematical (hidden)
symmetry of the system. To get the profiles admitting the symmetries in the
presence of polarization, we use the following observation lnp  (see mardoyan
for its quantum counterpart). Assume we have the three-dimensional
rotationally-invariant system

$\mathcal{H}_{0}=\frac{{p}^{2}}{2g(r)}+V(r),\qquad\\{p_{i},x_{j}\\}=\delta_{ij},\quad\\{p_{i},p_{j}\\}=\\{x_{i},x_{j}\\}=0.$
(25)

For the inclusion of interaction with magnetic monopole, we should transit
from the canonical Poisson brackets to the twisted ones:

$\\{p_{i},x_{j}\\}=\delta_{ij},\qquad\\{p_{i},p_{j}\\}=s\varepsilon_{ijk}\frac{x_{k}}{r^{3}},\qquad\\{x_{i},x_{j}\\}=0.$
(26)

The rotation generators then read

$\mathbf{J}=\mathbf{r}\times\mathbf{p}+s\frac{{\mathbf{r}}}{r}\;:\quad\\{J_{i},J_{j}\\}=\varepsilon_{ijk}J_{k}.$
(27)

By modifying the initial Hamiltonian to

$\mathcal{H}_{s}=\frac{{p}^{2}}{2g(r)}+\frac{s^{2}}{2g(r){r}^{2}}+V(r),$ (28)

we find that trajectories of the system preserve their form, but the plane
which they belong to, fails to be orthogonal to the the axis ${\mathbf{J}}$.
Instead, it turns to the constant angle

$\cos\theta_{0}=\frac{s}{|\mathbf{J}|}.$ (29)

For the systems with hidden symmetries one can find the appropriate
modifications of the hidden symmetry generators respecting the inclusion of
the monopole field.

For applying this observation on the systems with polarized light, we should
choose the appropriate integrable system with magnetic monopole, and then
perform the canonical transformation (16) which yields the Poisson brackets
for polarized light (26). Afterwards we need to solve the following equation

$r^{2}+\frac{s^{2}}{p^{2}}-2g(p)(E-V(p))=0,\quad\Rightarrow\quad
p=\frac{n^{s}_{inv}(r)}{\lambdabar_{0}}.$ (30)

For example, to get the “polarized Coulomb profile” we have to start from the
free-particle Hamiltonian on three-dimensional sphere/hiperboloid interacting
with Dirac monopole:

$H_{s}=\frac{(1+\kappa
r^{2})^{2}}{2}\left({p}^{2}+\frac{s^{2}}{r^{2}}\right).$ (31)

Then, after fixing the energy surface $H_{s}=E$ and performing canonical
transformation (16) we arrive to the third-order (with respect to ${p}^{2}$)
algebraic equation:

$(1+\kappa
p^{2})^{2}\left(r^{2}+\frac{s^{2}}{p^{2}}\right)=2E(>0),\qquad{\Leftrightarrow}\qquad
y^{3}u-y^{2}(u-\kappa s^{2})-Ey+E=0,$ (32)

with $y:=1+\kappa p^{2}$, $u:=r^{2}$.

This equation could have either one real and two complex solutions or three
real solutions, which will describe the “polarized Coulomb profiles”. When we
start from the opposite direction, that is when we start from the Coulomb
problem with Dirac monopole we will arrive to the ”polarized Maxwell fish
eye”, i.e. the deformation of the ”Maxwell fish eye” which preserves, in the
presence of polarized light, all symmetries of initial scalar system. The
latter is considered in detail in the next Section.

##  IV Polarized Maxwell fish eye

Let us consider the Coulomb system with Dirac monopole which is known as
“MICZ-Kepler system” MICZ . It is defined by the twisted Poisson brackets (26)
and by the Hamiltonian

$H_{MICZ}=\frac{{p}^{2}}{2}+\frac{s^{2}}{2r^{2}}-\frac{\gamma}{r}.$ (33)

Besides the conserved angular momentum (27), this system has the conserved
Runge-Lenz vector

$\mathbf{A}_{s}=\mathbf{J}\times\mathbf{p}+\gamma\frac{\mathbf{r}}{r},$ (34)

which form the symmetry algebra of Coulomb problem (20) (with the replacement
$(\mathbf{L},\mathbf{A})\to(\mathbf{J},\mathbf{A}_{s})$). After performing
canonical transformation (16), we get

$H_{MICZ}=E\quad\Leftrightarrow\quad
r^{2}+\frac{s^{2}}{p^{2}}-\frac{2\gamma}{p}-2E=0.$ (35)

Solvinq this quadratic (on $p$) equation, we get the refraction index given by
the expression (4), where the notation (18) is used.

The rotation generator (27) transforms to (23), and the Runge-Lenz vector (34)
transforms to ${\mathbf{T}_{s}}/{\kappa}$, where

$\mathbf{T}_{s}=\Big{(}2-\frac{n_{0}}{n^{s}_{Mfe}(\mathbf{r})}\Big{)}\mathbf{p}+2\kappa(\mathbf{rp})\mathbf{r}+\frac{2\kappa
s}{n^{s}_{Mfe}(\mathbf{r})}\mathbf{J}.$ (36)

Along with (23), these generators form the symmetry algebra of the original
Maxwell fish eye profile (20) (where $(\mathbf{L},\mathbf{T})$ are replaced by
$(\mathbf{J},\mathbf{T}_{s})$. The Casimirs of symmetry algebra are given by
the expressions

$\mathbf{T}^{2}_{s}+4\kappa(\mathbf{J}^{2}-s^{2})=\frac{n^{2}_{0}}{\lambdabar^{2}_{0}},\qquad\mathbf{T}_{s}\cdot\mathbf{J}=\frac{sn_{0}}{\lambdabar_{0}}.$
(37)

Hence, for $\kappa>0$ the vectors $\sqrt{4\kappa}\mathbf{J}$ and
$\mathbf{T}_{s}$ form the parallelogram with the fixed lengths of diagonals

$|\mathbf{T}_{s}\pm\sqrt{4\kappa}\mathbf{J}|=|\frac{n_{0}}{\lambdabar_{0}}\pm\sqrt{4\kappa}s|.$
(38)

This immediately leads to the conclusion that for $\kappa>0$ the generators
$\mathbf{T}_{s}$ and $\mathbf{J}$ reach the lower/upper bounds being parallel
to each other

$\left(|\mathbf{J}|_{\rm min}=s,\;|\mathbf{T}_{s}|_{\rm
max}=\frac{n_{0}}{\lambdabar_{0}}\right),\qquad\left(|\mathbf{J}|_{\rm
max}=\frac{n_{0}}{\lambdabar_{0}\sqrt{4\kappa}}\;,|\mathbf{T}_{s}|_{\rm
min}=\sqrt{4\kappa}s\right).$ (39)

Notice also, that for $\kappa>0$ we get a restriction of rays in the finite
domain

$\kappa>0\;:\quad
r\leq\sqrt{\frac{n^{2}_{0}}{4s^{2}\lambdabar^{2}_{0}\kappa^{2}}-\frac{1}{\kappa}}.$
(40)

One can also note that spin appears in the expression for the refraction index
(4) along with the factor
$\kappa{\lambda_{0}}^{2}=({\lambda_{0}}/{2r_{0}})^{2}$. In order to stay
within the bounds of geometrical optics approximation, this factor must be
reasonably small. Therefore, the influence of the spin will be far more
notable within certain range of distance from the core of the fish eye. The
latter happens when the condition ${4\kappa
s^{2}\lambdabar^{2}_{0}}/{n_{0}}\approx n_{Mfe}(\mathbf{r})$ holds. At these
distances the refraction index in presence of spin can be two times smaller as
compared to the refraction index with zero spin.

Figure 1: Maxwell fish eye refraction index profile for $s=0$ and $s=1$ when
$n_{0}=1.5$, $\lambda_{0}=1$, $r_{0}=5$.

##  V Trajectories

Let us study the ray trajectories of the polarized light propagating in the
medium with above constructed profile (4). One can see that

$\mathbf{r}\cdot\mathbf{J}=s\frac{\mathbf{rp}}{p},\quad\mathbf{r}\cdot\mathbf{T}_{s}=\frac{n_{0}}{\lambdabar_{0}}\frac{\mathbf{rp}}{p}\qquad\Rightarrow\qquad\mathbf{r}\cdot\left(\mathbf{J}-\frac{s\lambdabar_{0}}{n_{0}}\mathbf{T}_{s}\right)=0.$
(41)

Hence, ray trajectories are orthogonal to the axis

$\mathbf{E}_{3}=\mathbf{J}-\frac{s\lambdabar_{0}}{n_{0}}\mathbf{T}_{s},$ (42)

and, therefore the trajectories belong to the plane spanned by the following
vectors

$\mathbf{E}_{1}=\mathbf{T}_{s}\times\mathbf{J},\qquad\mathbf{E}_{2}=\mathbf{E}_{3}\times\mathbf{E}_{1}=\Big{(}\mathbf{J^{2}}-s^{2}\Big{)}\left(\mathbf{T}_{s}-\frac{4s\lambdabar_{0}\kappa}{n_{0}}\mathbf{J}\right)\quad:\quad\mathbf{E}_{3}\cdot\mathbf{E}_{2}=\mathbf{E}_{3}\cdot\mathbf{E}_{1}=0.$
(43)

Then, from the expression
$\mathbf{J}\cdot\big{(}\mathbf{r}\times\mathbf{T_{s}}\big{)}$ we immediately
obtain the solution for the ray trajectories:

$\mathbf{r}\cdot\big{(}\mathbf{T_{s}}\times\mathbf{J}\big{)}=\Big{(}J^{2}-s^{2}\Big{)}\Big{(}2-\frac{n_{0}}{n^{s}_{mfe}}\Big{)}.$
(44)

This prompts us to introduce the following orthogonal frame

$\mathbf{e}_{i}=\frac{\mathbf{E}_{i}}{|\mathbf{E}_{i}|}\;:\quad\mathbf{e}_{i}\cdot\mathbf{e}_{j}=\delta_{ij},$
(45)

where

$|\mathbf{E}_{1}|^{2}=\left(\mathbf{J}^{2}-s^{2}\right)\left(\frac{n^{2}_{0}}{\lambdabar^{2}_{0}}-{4\kappa\mathbf{J}^{2}}\right),\qquad|\mathbf{E}_{3}|^{2}=\left(\mathbf{J}^{2}-s^{2}\right)\left(1-\frac{4s^{2}\lambdabar^{2}_{0}\kappa}{n^{2}_{0}}\right),\qquad|\mathbf{E}_{2}|^{2}=|\mathbf{E}_{1}|^{2}|\mathbf{E}_{3}|^{2}$
(46)

Decomposing $\mathbf{r}$ over this frame, we introduce the polar coordinates

$\mathbf{r}=x_{1}\mathbf{e}_{1}+x_{2}\mathbf{e}_{2},\qquad
x_{1}=r\cos\varphi,\quad x_{2}=r\sin\varphi.$ (47)

Then, having in mind Eqs. (43) and (46), we can immediately rewrite the
equation (44) in polar coordinates

$1-|\kappa||\mathbf{a}_{s}|r\cos\varphi=\frac{1+\kappa
r^{2}}{1+\sqrt{1-\frac{4\kappa s^{2}\lambdabar^{2}_{0}}{n^{2}_{0}}(1+\kappa
r^{2})}},$ (48)

where

$|\mathbf{a}_{s}|^{2}:=R^{2}_{s}-\frac{1}{\kappa},\qquad
R^{2}_{s}:=\frac{n^{2}_{0}-4\kappa
s^{2}\lambdabar^{2}_{0}}{4\lambdabar^{2}_{0}\kappa^{2}(J^{2}-s^{2})}.$ (49)

So, trajectories of polarized light are not circles anymore, in contrast to
the case of scalar waves. However they can be attributed by the parameters
${\bf a_{s}}$ and $R_{s}$ (49) which result center coordinate and radius of
circle, respectively, when $s\to 0$. Indeen, for $s=0$ the solution (44)
results in the the equation for circle with the center located at
$\mathbf{e}_{1}$ axis

$s=0\;:\quad\mathbf{r}\cdot\big{(}\mathbf{T}\times\mathbf{L}\big{)}=L^{2}(1-\kappa
r^{2})\quad\Rightarrow\quad\left(\mathbf{r}-\mathbf{a}_{0}\right)^{2}=R^{2}_{0},$
(50)

where

$\mathbf{a}_{0}:=\frac{\mathbf{T}\times\mathbf{L}}{2\kappa
L^{2}},\quad|\mathbf{a}_{0}|^{2}=R^{2}_{0}-\frac{1}{\kappa},\quad
R_{0}:=\frac{n_{0}}{2|\kappa|\lambdabar_{0}L}.$ (51)

Figure 2: Deformations of the ray trajectories for different values of $r_{0}$
when $n_{0}=1.5$, $\lambdabar_{0}=1$. The black curves correspond to the basic
trajectories where $T=T_{min}=s/r_{0}$. The red curves correspond to
trajectories with intermediary value of $T$. Dashed curves are the
trajectories corresponding to the same value of $T$ but for circularly
polarized light ($s=1$). The first figure ($r_{0}=20$) corresponds to the case
when the deformations of the profile only result to the shift of the centers
of the trajectories not affecting their shapes. Conversely, in the second
figure ($r_{0}=2$), the deformations of the Maxwell fish eye result to highly
deformed trajectories.

For a small $s\lambdabar_{0}\ \ll n_{0}/\kappa$ the deformation of circle is
negligible, which is not the case for $s\lambdabar_{0}\ \sim
n_{0}/\sqrt{\kappa}$.

As we can see from Fig.2, for $r_{0}=20$ the only notable manifestation of the
polarization is the shift of the center of the trajectory. However, it is
worth noting that since the dashed lines are not circles anymore. Talking
about the center of the trajectory for $s=1$, we refer to the point which
becomes the center of the circular trajectory when we pass from $s=1$ to
$s=0$. The second picture illustrates the circular trajectories and their
deformations for $s\lambdabar_{0}\ \sim n_{0}/\sqrt{\kappa}$. In this case the
original profile (3) and the deformed one (4) differ drastically, so that
circular trajectories are notably deformed.

###### Acknowledgements.

Authors acknowledge partial financial support from Armenian Committee of
Science, project 20RF-023 (A.N., Zh.G.), and from the Russian Foundation of
Basic Research grant 20-52-12003 (A. N.). This work was done within ICTP
projects NT-04 and AF-04.

## References

  * (1) V. I. Arnold, _Mathematical methods in classical mechanics_ , Graduate Texts in Mathematics, Vol. 60 (Second edition), 1989, Springer-Verlag, New York
  * (2) J. Pendry, D. Schurig, and D. Smith, _Controlling electromagnetic fields_ , Science 312, 1780 (2006).
  * (3) U. Leonhardt, _Optical conformal mapping_ , Science 312, 1777 (2006);
  * (4) U. Leonhardt,_Notes on conformal invisibility devices_ , New J. Physics,8, 118 (2006).
  * (5) Z. Gevorkian and M. Davtyan, _Discrete spectrum radiation from a charged particle moving in a medium with Maxwell’s fish-eye refraction-index profile,_ Phys. Rev. A 102, 063504 (2020)
  * (6) J. Perczel, P. K$\acute{o}$m$\acute{a}$r, and M. D. Lukin,_Quantum optics in Maxwell’s fish eye lens with single atoms and photons,_ Phys. Rev. A 98, 033803 (2018).
  * (7) K. Dadashi, H. Kurt, K. Ustun and R. Esen, _Graded index optical microresonators: analytical and numerical analyses,_ J. Opt. Soc. Am. B 31, 2239 (2014).
  * (8) Zh. Jin, Y. D. Wang, L. Z. Yin, F. Y. Han, T. J. Huang and P. K. Liu, _Bifurcational Luneburg-fish eye lens based on the manipulation of spoof surface plasmons,_ Opt. Lett. 46, 1389 (2021).
  * (9) H. Lu, G. Wu, Y. Liu and X. Lv, _A Millimeter-Wave Fully Metallic Six-Channel Crossover Based on Maxwell Fish-Eye Lens,_ IEEE Microwave and Wireless Components Letters, 30, 1041 (2020).
  * (10) R.A.S.D.Koala, D.Headland,Y.Yamagami, F.Masayuki and T.Nagatsuma, _Broadband Terahertz Dielectric Rod Antenna Array with Integrated Half-Maxwell Fisheye Lens,_ 2020 International Topical Meeting on Microwave Photonics (MWP) (pp. 54-57), (2020) DOI: 10.23919/MWP48676.2020.9314424
  * (11) Z. Gevorkian, M. Davtyan and A. Nersessian, _Extended symmetries in geometrical optics,_ Phys. Rev. A 101, 023840 (2020)
  * (12) P.A.M. Dirac, The principles of quantum mechanics, (Oxford, Clarendon, 4th ed. 1958) A. Deriglazov, Classical Mechanics: Hamiltonian and lagrangian Formalism, Springer, 2010
  * (13) K. Y. Bliokh, Avi Niv, V. Kleiner, E. Hasman, _Geometrodynamics of spinning light_ , Nature Photonics, 2, 748-753 (2008). K. Y. Bliokh, _Geometrodynamics of polarized light:Berry phase and spin Hall effect in a gradient-index medium_ , J. Opt. A: Pure Appl. Opt. 11, 094009 (2009).
  * (14) Y. A. Kravtsov and Y. I. Orlov, _Geometrical optics of inhomogeneous medium_ , Springer, 1990
  * (15) A. M. Perelomov, _Integrable Systems of Classical Mechanics and Lie Algebras_ , Springer, 1990
  * (16) A. Nersessian, _Elements of (super-)Hamiltonian formalism,_ Lect. Notes Phys. 698, 139 (2006) [hep-th/0506170].
  * (17) L. Mardoyan, A. Nersessian and A. Yeranyan, _Relationship between quantum mechanics with and without monopoles,_ Phys. Lett. A 366, 30 (2007)
  * (18) D. Zwanziger, _Quantum field theory of particles with both electric and magnetic charges,_ Phys. Rev. 176, 1489 (1968) H. V. Mcintosh and A. Cisneros, _Degeneracy in the presence of a magnetic monopole,_ J. Math. Phys. 11, 896 (1970) .


# [Reproducibility Report] Path Planning using Neural A* Search

Shreya Bhatt
<EMAIL_ADDRESS>Aayush Jain
<EMAIL_ADDRESS>Parv Maheshwari
<EMAIL_ADDRESS>Animesh Jha
<EMAIL_ADDRESS>Debashish Chakravarty
<EMAIL_ADDRESS>
## Reproducibility Summary

The following paper is a reproducibility report for "Path Planning using
Neural A* Search" [1] published in ICML 2021 as part of the ML Reproducibility
Challenge 2021. The source code 111https://github.com/shreya-
bhatt27/NeuralAstar-ported for our reimplementation and additional experiments
performed is available for running.

### Scope of Reproducibility

The original paper proposes the Neural A* planner, and claims it achieves an
optimal balance between the reduction of node expansions and path accuracy. We
verify this claim by reimplementing the model in a different framework and
reproduce the data published in the original paper. We have also provided a
code-flow diagram to aid comprehension of the code structure. As extensions to
the original paper, we explore the effects of (1) generalizing the model by
training it on a shuffled dataset, (2) introducing dropout (3) implementing
empirically chosen hyperparameters as trainable parameters in the model, (4)
altering the network model to Generative Adversarial Networks (GANs) to
introduce stochasticity, (5) modifying the encoder from Unet to Unet++ (6)
incorporating cost maps obtained from the Neural A* module in other variations
of A* search.

### Methodology

We reimplemented the publicly available source code provided by the authors in
Pytorch Lightning to encourage reproducibilty of the code and flexibility over
different hardware setups. We reproduced the results published by the authors
and also conducted additional experiments. The training code was run on Kaggle
with GPU (Tesla P100-PCIE-16GB) and CPU (13GB RAM + 2-core of Intel Xeon).

### Results

The claims of the original paper were successfully reproduced and validated
within 3.2% of the reported values. Results for additional experiments
mentioned above have also been included in the report.

### What was easy

The code provided in the original repository was well structured and
documented making it easy to understand and reimplement. The authors also
provide the source code for dataset generation which made the task of
reproducing the results fairly simple .

### What was difficult

Experimentation on some datasets took a considerable amount of time, limiting
our experiments to the MP Dataset. Results of runtime calculation could not be
reproduced as they are affected by various factors, including dissimilarity in
datasets, hardware environment and A* search implementation.

### Communication with original authors

The authors were contacted via email regarding the computational requirements
of training and errors faced, to which prompt and helpful replies were
received.

##  1 Introduction

As described in the original paper, A* search is a path planning algorithm
that is frequently used in vehicle navigation, and robot arm manipulation to
solve the lowest cost path problems. The paper then goes on to introduce
Neural A*, a unique data-driven path planning system. The Neural A* planner
uses a differentiable A* module with a fully convolutional encoder. The
encoder produces a costmap, which is used by the search algorithm. The network
learns to recognize visual cues from input maps which are effective in
producing ground truth paths.

In this work, we reimplement the Neural A* planner (originally in Pytorch) in
Pytorch Lightning [2] to make it easier to reproduce the research in the
future. This allowed us to train the code flexibly over different platforms
and automate the optimization process. We utilized the functionalities of
PyTorch Lightning for efficient checkpointing and logging on Wandb. We check
the reproducibility of the claims published in the original paper and perform
multiple ablations. We then go on to explore the performance of different
variations of A* coupled with the Neural A* costmap.

One of the limitations of the original work was the lack of diversity in
generated paths, even though there are often multiple optimal paths. Hence in
an attempt to extend the planner to account for stochastic sampling of paths
we experimented with GAN [3]. We also added an implementation with Unet++
encoder [4].

##  2 Scope of Reproducibility

The original paper attempts to compare the performance of the proposed Neural
A* planner with other data-driven and imitation based planners. In this work,
we check the reproducibility of the following claims.

  * •

Comparison with other planners: Neural A* successfully finds near-optimal
paths while significantly reducing the number of nodes explored compared to
other planners. The differentiable A* reformulation of the Neural A* planner
is key in achieving optimal performance.

  * •

Start and goal maps: The provision of start and goal maps to the Neural A*
encoder enables the extraction of visual cues properly conditioned to start
and goal locations, leading to significant improvement in path optimality and
reduction of nodes explored.

  * •

VGG-16 backbone: Using VGG-16 backbone [5] leads to better performance than
ResNet-18 [6].

  * •

Path length optimality: Neural A* produces nearly optimal path lengths over
several datasets.

  * •

Planning on raw images: Neural A* facilitates path planning on raw images
without semantic pixel-wise labelling of the environment. The model
outperforms several other planners in predicting pedestrian trajectories from
raw surveillance data.

We validate the above claims by reproducing all results from the main article
and appendix (excluding SDD dataset) except for runtime calculations (further
explained under Section 6.2). In addition, we had performed several
experiments to study the efficiency of the planner when coupled with different
network models, mixed and shuffled datasets and variations of the simple A*
search algorithm. We also performed hyperparameter tuning and implemented
empirically chosen constants to instead be trainable parameters learnt through
backpropagation. Section 5 contains more details.

##  3 Methodology

The Neural A* and Neural BF planners were reimplemented by us, taking
inspiration from the official codebase 222https://github.com/omron-
sinicx/neural-astar

###  3.1 Model descriptions

The simple A* search operates by maintaining an open list of nodes. The search
then alternates between searching the open list for the node most likely to
lead to a low cost path and expanding the list by appending the neighbours of
the selected node until the goal is reached. The search uses a combination of
the total best path cost and a heuristic function for node selection. The A*
search algorithm has been explained in great detail in the original paper. An
encoder in the Neural A* model converts a problem instance to a guidance map
which assigns a guidance cost to each node. The Differentiable A* module then
performs an A* search using the guidance cost along with the heuristic as part
of the total cost. As inspired by previous chapters of MLRC [7] the codeflow
for the model has been provided (Fig. 1).

Figure 1: Codeflow of Neural A*

Description of the planner

The Neural A* planner consists of 2 modules:

  1. 1.

Encoder : UNet [8], a fully convolutional encoder network, is employed. Visual
cues such as the contours of dead ends and by-passes are learned by the
encoder. The encoder input is a concatenation of the environmental map with
the start and goal maps (binary maps indicating the locations of start and
goal nodes). The encoder then generates a guidance map, which is a scalar
valued map representation of the problem instance. During A* search, the
guidance map is used as part of the heuristic.

  2. 2.

Differentiable A* module : The differentiable A* module is a differentiable
reformulation of the classic A* search method. This was accomplished by using
a discretized activation technique inspired by [9] with simple matrix
operations. The module performs an A* search in the forward pass and then
backpropagates the losses to other trainable modules after each epoch. This
module’s goal is to teach the encoder how to make guidance maps that reduce
the difference between search histories and ground-truth paths.

###  3.2 Loss function

Let us consider the closed list $C$, which is a binary matrix containing all
the searched nodes $V$. The optimal path P is a binary matrix containing the
nodes which belong to the ground-truth path. The mean L1 loss between $C$ and
the optimal path map $P$ is calculated as follows [1]:

$\displaystyle L=\frac{||C-P||_{1}}{|V|}$

The loss defined as above penalizes both paths that stray from the optimal
path (calculated by running a simple A* search, referred to as Vanilla A*) and
the exploration of nodes not part of the optimal path. This loss motivates
Neural A* to choose paths that are similar to the ground-truth path while also
reducing node explorations.

###  3.3 Datasets

The MP, TiledMP and CSM datasets are used and preprocessed identical to the
original implementation. We were unable to setup the pipeline for
reimplementing results for SDD dataset due to resource constraints which is
further discussed in section 6.2. The source code333https://github.com/omron-
sinicx/planning-datasets for generating the datasets was made publicly
available by the authors. In addition, we created a mixed dataset, mixing maps
from various subsets of the MP dataset to check the generalizability of the
model. This is described further in section 5.1.

###  3.4 Hyperparameters

The hyperparameters and their respective ranges are same as mentioned in the
original paper. The hyperparameters were set after conducting random search
using Wandb [10]. The best hyperparameters are $g\\_ratio=0.2$ and rest of the
hyperparameters same as the paper. As can be observed from the values, the
changes in hyperparameters are in accordance with weighted A*, decreasing the
weight of the heuristic (= $1-g\\_ratio)$ , the number of explorations and
path optimality increase. On the other hand increasing the weight of the
heuristic leads to a decrease in both number of explorations and path
optimality. Note, number of explorations and EXP are not the same. Exp
(Section 3.5) is inversely proportional to number of explorations. Hence, it
can be seen that Exp and Opt (Section 3.5) have an inverse relationship, and
that it is important to find hyperparameters that maintain the equilibrium
between them.Encoder input can be of the type m+ (includes start maps, goal
maps and map designs), or m (includes just the map designs). Note that the
original paper does not report results for Resnet backbone with Encoder input
of type $m$.

Table 1: VGG Backbone with the best hyperparameters Input Type |  Version |  Opt |  Exp |  Hmean |  Suc |
---|---|---|---|---|---|---
m | Best | 64.95 (60.22, 69.87) | 42.80(39.52, 46.07) | 45.44(42.03, 48.90) | 100.0 |
| Original | 67.0 (65.1, 68.8) | 36.8(35.6, 38.1) | 41.5(40.2, 42.7) | 100.0 |
m+ | Best | 78.97(75.23, 82.97) | 44.49 (41.28, 47.68) | 52.42 (49.14, 55.77) | 100.0 |
| Original | 87.7 (86.6, 88.9) | 40.1(38.9, 41.3) | 52.0(50.7, 53.1) | 100.0 |
Table 2: Resnet Backbone with best hyperparameters Input Type |  Version |  Opt |  Exp |  Hmean |  Suc |
---|---|---|---|---|---|---
m | Best | 64.10 (59.34, 69.07) | 42.40 (39.14, 45.76) | 44.60 (41.15, 48.12) | 100.0 (100.0, 100.0) |
| Original | — | — | — | — |
m+ | Best | 75.57 (71.88, 79.83) | 43.71 (40.57, 46.85) | 50.80 (47.77, 53.93) | 100.0 (100.0, 100.0) |
| Original | 79.8(78.1, 81.5) | 41.4(40.2, 42.7) | 49.2(47.9, 50.5) | 100.0 (100.0, 100.0) |

###  3.5 Experimental setup and code

The models were trained using the RMSProp optimizer with a batch size of 100
with a learning rate 0.001. The models were trained for 100 epochs on the MP
Dataset and 400 epochs for the CSM and TiledMP datasets [11, 12, 13]. The
following metrics are as defined in the original paper and evaluated for each
trained model on their test set performance.

  * •

Path optimality ratio (Opt): measures the percentage of shortest path
predictions for each map.

  * •

Reduction ratio of node explorations (Exp): measures the number of search
steps reduced by a model as compared to vanilla A* search. Let $E^{*}$ be the
number of nodes explored by vanilla A* and $E$ be the number of nodes explored
by a particular model. Then Exp is defined as
max$(100(\frac{E^{*}-E}{E^{*}},0))$.

  * •

Harmonic mean of Opt and Exp (Hmean): Indicates the trade-off between Opt and
Exp achieved.

A higher Hmean value is desired and indicates that a model achieves satisfying
optimality and search efficiency. The central planners implemented in this
work are Neural A* and Neural BF (i.e. Neural Best First which is simply the
Neural A* planner with hyperparameter g_ratio set to 0). Various other
planners were used as baselines for judging the performance of Neural A*
including SAIL [14], BBA* [15, 1], Weighted A* [16], Best First Search, VIN
[17], GPPN [18] and MMP [19]. The values of these metrics for baseline
planners were obtained from running the official source code and were
identical to those published in the original paper.

###  3.6 Computational requirements

The training code was run on Kaggle with GPU (Tesla P100-PCIE-16GB) and CPU
(13gb RAM + 2-core of Intel Xeon). The average training times were $1.25$, 72
and 58 hours for the MP, Tiled MP and CSM datasets respectively.

##  4 Results

Detailed description of the experiments and results supporting the claims made
in Section 2 is given below

###  4.1 Results reproducing original paper

We compare the Bootstrap means and 95% confidence intervals of path optimality
ratio, reduction ratio of node explorations and the harmonic mean of these
values from the ported code (R) and original paper (O).

Table 3: Ablation Results Dataset |  Planner |  Version |  Opt |  Exp |  Hmean |  Suc |
---|---|---|---|---|---|---|---
MP | Neural A* | R | 86.15 (85.00, 87.37) | 39.96 (38.74, 41.18) | 51.50 (50.23, 52.80) | 100.0 |
|  | O | 87.7 (86.6, 88.9) | 40.1 (38.9, 41.3) | 52.0 (50.7, 53.3) | 100.0 |
| Neural BF | R | 74.29 (72.70, 75.90) | 44.97 (43.66, 46.27) | 50.93 (49.64, 52.23) | 100.0 |
|  | O | 75.5 (73.8, 77.1) | 45.9 (44.6, 47.2) | 52.0 (50.7, 53.4) | 100.0 |
Tiled MP | Neural A* | R | 74.84 (71.20, 78.69) | 48.86 (45.25, 52.44) | 55.66 (52.56, 58.73) | 100.0 |
|  | O | 63.0 (60.7, 65.2) | 55.8 (54.1, 57.5) | 54.2 (52.6, 55.8) | 100.0 |
| Neural BF | R | 44.17 (40.95, 48.34) | 61.3 (57.74, 64.92) | 44.61 (39.37, 49.00) | 100.0 |
|  | O | 43.7 (41.4, 46.1) | 61.5 (59.7, 63.3) | 44.4 (42.5, 46.2) | 100.0 |
CSM | Neural A* | R | 77.82 (74.07, 81.75) | 34.99 (31.06, 38.86) | 42.64 (38.79, 46.51) | 100.0 |
|  | O | 73.5 (71.5, 75.5) | 37.6 (35.5, 39.7) | 43.6 (41.7, 45.5) | 100.0 |
| Neural BF | R | 82.70(81.29, 84.15) | 39.16 (37.92, 40.38) | 49.27 (47.98, 50.55) | 100.0 |
|  | O | 79.8 (78.1, 81.5) | 41.4(40.2, 42.7) | 49.2(47.9, 50.5) | 100.0 |

####  4.1.1 Path length optimality

Path length optimality is another metric used to judge the performance of the
planner. Path length optimality P is defined as the percent ratio of optimal
path length and predicted path length.

Table 4: Path length optimality comparison Dataset |  Planner |  Reimplementation |  Original Paper |
---|---|---|---|---
MPD | Neural A* | 99.23 (98.99, 99.51) | 99.1 (99.0, 99.2) |
| Neural BF | 98.04 (97.65, 98.45) | 97.5 (97.3, 97.8) |
TiledMP | Neural A* | 99.12 (98.89, 99.39) | 98.4 (98.3, 98.6) |
| Neural BF | 96.85 (96.25, 97.47) | 95.0 (94.7, 95.4) |
CSM | Neural A* | 99.31 (99.09, 99.55) | 98.9 (98.8, 99.0) |
| Neural BF | 97.62 (97.14, 98.13) | 97.4 (97.1, 97.6) |

####  4.1.2 Ablation Results from the Original Paper

Comparison of results of reproducing ablations from original paper.

  1. 1.

The encoder was not provided start and goal node locations

  2. 2.

Resnet-18 backbone was used instead of VGG-16

Table 5: Ablation results comparison Ablation |  Version |  Opt |  Exp |  Hmean |  Suc |
---|---|---|---|---|---|---
Ablation 1 | R | 70.07 (68.14, 72.01) | 39.56(38.27, 40.83) | 45.89(44.54, 47.27) | 100.0 |
| O | 67.0 (65.1, 68.8) | 36.8(35.6, 38.1) | 41.5(40.2, 42.7) | 100.0 |
Ablation 2 | R | 82.70(81.29, 84.15) | 39.16 (37.92, 40.38) | 49.27 (47.98, 50.55) | 100.0 |
| O | 79.8 (78.1, 81.5) | 41.4(40.2, 42.7) | 49.2(47.9, 50.5) | 100.0 |

##  5 Additional Experiments

For additional experiments we suggest the use of "number of explorations" as
an additional metric. Since Exp is a relative metric that compares expansions
by a particular model with Vanilla A*, the degraded performance of Vanilla A*
can sometimes cause misleading results. On the other hand, "number of
explorations" is an intuitive and independent metric making it easier to
compare models. The author also suggested the use of number of explorations as
a proxy metric that is likely to be used for neural planners. We use NA* Exps
as number of explorations by Neural A* and VA* Exps as number of explorations
by Vanilla A*.

###  5.1 Shuffling Dataset

In the original implementation, the authors train the Neural A* model on
subsets of the MP dataset as individual environments. This leads to the model
working efficiently on few maps, but generalisation performance of the model
suffers. In this experiment we choose some maps from each of the MP subsets
and create a mixed set. We then train the model on this dataset and compare
it’s results with those on individual environments. From Table 6 we observe
that the model achieved better performance on a generalized test set as
compared to individually trained models.

Table 6: Mixed set comparison Sub-dataset | NA* Exps | VA* Exps | Opt | Exp | Hmean |
---|---|---|---|---|---|---
alternating gaps | 58.51 | 102.89 | 62.24 (56.06, 68.85) | 36.21 (32.48, 39.83) | 38.14 (33.96, 42.34) |
bugtrap forest | 76.23 | 98.73 | 73.14 (68.51, 78.13) | 24.93 (21.52, 28.20) | 32.38 (28.82, 35.99) |
forest | 66.56 | 91.15 | 73.80 (69.62, 78.22) | 27.18 (24.13, 30.09) | 35.69 (32.22, 39.11) |
gaps and forest | 61.77 | 89.80 | 71.48 (66.95, 76.16) | 30.35 (26.84, 33.92) | 39.20 (35.42, 42.92) |
mazes | 79.79 | 104.95 | 71.83 (67.19, 76.68) | 27.10 (23.28, 30.72) | 34.84 (30.92, 38.82) |
multiple bugtraps | 75.04 | 99.11 | 77.39 (73.20, 81.89) | 25.70 (22.10, 29.17) | 33.44 (29.61, 37.30) |
shifting gaps | 61.04 | 99.55 | 70.56 (65.17, 76.27) | 32.56 (28.52, 36.78) | 38.70 (33.96, 43.39) |
single bugtrap | 101.8 | 94.02 | 71.91 (66.95, 77.04) | 18.11 (14.60, 21.27) | 23.13 (19.44, 26.82) |
Average for MP dataset | 72.59 | 97.52 | 71.44 (69.65, 73.30) | 27.29 (26.00, 28.60) | 33.98 (32.53, 35.40) |
Mixed set | 57.74 | 94.06 | 84.32 (80.50, 88.39) | 35.51 (32.05, 38.85) | 45.97 (42.25, 49.71) |

###  5.2 Dropout

Experimented by adding dropout to UNet taking inspiration from [20] where it
is claimed that the dropout introduces stochasticity at the planning step. We
believe that further performance improvements can be made by implementing a
replanning step as explained in section 7. Here, g_ratio = 0.5 in R
(implementation with dropout) as comparing with O (original paper that has
g_ratio = 0.5)

Table 7: Dropout Results for g_ratio = 0.5  Input Type |  Version |  NA* Exps |  VA* Exps |  Opt |  Exp |  Hmean |  Suc
---|---|---|---|---|---|---|---
m | R | 66.63 | 103.01 | 72.58 (68.25, 77.17) | 31.43(28.14, 34.65) | 38.80(35.34, 42.23) | 100.0
| O | 68.18 | 103.01 | 67.00 (65.10, 68.80) | 36.80(35.60, 38.10) | 41.50(40.20, 42.70) | 100.0
m+ | R | 59.08 | 103.01 | 84.71(81.61, 88.05) | 37.89(34.83, 40.92) | 49.03(45.89, 52.21) | 100.0
| O | 56.27 | 103.01 | 87.70 (86.60, 88.90) | 40.10(38.90, 41.30) | 52.00(50.70, 53.10) | 100.0

###  5.3 Suggested Experiments by the Authors

The authors have mentioned in the paper that if there are several paths to the
goal location, the planner fails to identify actual pedestrian trajectories.
They also suggested that adopting a generative framework [21] that permits
numerous paths to be stochastically sampled could be a suitable extension to
address this issue [22]. While communicating with the authors, we were
suggested to experiment by changing the encoder to a generative model based
architecture like GAN. The authors also suggested that learning the value of
temperature via backpropagation is an interesting direction to explore.

####  5.3.1 Temperature Function

The temperature function $\tau$ used in the differentiable A* module is
empirically chosen by the authors. In this experiment we implement $\tau$ as a
trainable parameter of the Neural A* model whose value gets updated by
backpropagation. It was observed that on initialization of $\tau$ to the
square root of map size, the value remains almost constant over the training
process. In addition, we also ran sweeps over a range of values of $\tau$. The
sweeps have only been conducted on one sub dataset of the MP dataset
(bugtrap_forest_32) due to computational restraints.

Table 8: Hyperparameter Sweep on Temperature Function with g_ratio = 0.2 Temperature |  NA* Exps |  VA* Exps |  Opt |  Exp |  Hmean |  Suc
---|---|---|---|---|---|---
$\sqrt{32}$(Original) | 46.02 | 69.63 | 78.32 ( 75.02, 81.84) | 38.95(36.37, 41.32) | 50.15(47.71, 52.63) | 100.0
$2\sqrt{32}$ | 46.48 | 84.22 | 78.20 ( 74.30, 82.31) | 36.50(33.50, 39.39) | 46.64 (43.73, 49.52) | 100.0
$3\sqrt{32}$ | 52.96 | 75.84 | 73.13(68.95, 77.60) | 41.98(39.53, 44.43) | 50.64 (47.92, 53.34) | 100.0
$4\sqrt{32}$ | 45.00 | 81.34 | 75.41(71.61, 79.44) | 41.34(38.71, 43.86) | 50.94 (48.31, 53.61) | 100.0
${\sqrt{32}}/{2}$ | 57.51 | 88.85 | 79.50 (75.74, 83.53) | 31.87 (29.48,34.18) | 43.20(40.60, 45.80) | 100.0
${\sqrt{32}}/{3}$ | 50.05 | 76.61 | 85.33 (82.21, 88.86) | 32.46 (30.22, 34.67) | 45.10(42.46, 47.76) | 100.0

####  5.3.2 GANs

General Adversarial Networks are a class of generative models that involves
two sub-models, a Generator which generates new examples from the domain and a
Discriminator which classifies the generated output as real or generated.
Generative models have been used for path planning before [23]. For our
experiments we use the vanilla GANs from Pytorch Lighting-Bolts [24] as the
encoder inplace of UNet and train it on the . From the results mentioned in
Table 9 and 10 we can infer that our GANs experiment combined with best
hyperparameters (g_ratio = 0.5) outperforms the original paper for maps with
input type as $m$. Further, for the same hyperparameters, it can be seen that
GANs with input type as $m$ performs slightly better than that with $m+$,
contrary to the original paper where input type $m+$ outperforms $m$ by a
large margin. We believe that in cases with multiple possible routes the
inherent stochasticity of the generator would help in producing guidance maps
that tempt the planner to try different possibilities and make it more robust.

###  5.4 UNet++

As advised by the author, we experimented with UNet++[4] architecture while
undertaking an architectural search and assessing its scope for further
performance improvements. UNet++ is similar to UNet, with an ensemble of
variable depth UNets which partially share an encoder and co-learn
simultaneously using deep supervision. However, we find that UNet performs
better than UNet++, showing that powerful segmentation models might not always
be better at learning guidance maps for Neural A*.

Table 9: Encoder Architecture Experimentation on MP Dataset with ’m’ type encoder input Encoder Type |  NA* Exps |  VA* Exps |  Opt |  Exp |  Hmean |  Suc
---|---|---|---|---|---|---
UNet | 68.18 | 103.01 | 67.00 (65.10, 68.80) | 36.80(35.60, 38.10) | 41.50(40.20, 42.70) | 100.0
UNet++ | 64.99 | 102.75 | 76.41(72.27, 80.79) | 31.79(28.52, 35.00) | 39.68(36.38, 43.02) | 100.0
GAN | 68.40 | 103.25 | 66.44(61.53, 71.54) | 38.27 (34.81, 41.68) | 42.12(38.6, 45.65) | 100.0
Table 10: Encoder Architecture Experimentation on MP Dataset with ’m+’ type encoder input Encoder Type |  NA* Exps |  VA* Exps |  Opt |  Exp |  Hmean |  Suc
---|---|---|---|---|---|---
UNet | 56.27 | 103.01 | 87.70 (86.60, 88.90) | 40.10(38.90, 41.30) | 52.00(50.70, 53.10) | 100.0
UNet++ | 57.13 | 102.75 | 87.46(84.67, 90.48) | 36.41(33.30, 39.49) | 48.17(44.97, 51.46) | 100.0
GAN | 67.74 | 103.25 | 65.54(60.60, 70.68) | 38.44(34.97, 41.86) | 41.74(38.20, 45.31) | 100.0

###  5.5 Variations of A*

In this experiment we test the performance of various versions of A* search
coupled with and without the costmaps generated by the Neural A* module. We
can infer from the obtained results that the costmaps generated by the Neural
A* planner can be used with other variations of the A* search. However, a
particular model works best with variations of A* using the weight of
heuristic similar to the model. As a result, we decided against incorporating
the costmaps generated using dynamic weighted A* because the model’s weight
varies while the search is ongoing.

Table 11: Variations of A* search on mazes (MP Dataset)  A* Variation |  Opt |  Exp |  Hmean |
---|---|---|---|---
Vanilla A*[25] | 57.99 (54.06, 62.01) | 21.48 (18.86, 23.95) | 24.71 (21.02, 28.29) |
Weighted A* | 44.45 (39.56, 49.36) | 17.80 (13.38, 22.04) | 18.15 (14.92, 21.18) |
Beam search[26] | 58.74 (55.04, 62.43) | 20.38 (17.85, 22.78) | 26.07 (22.69, 29.40) |

##  6 Discussion

Our results reproduce the findings of [1] and strengthen the claims made in
section 2. From section 4.1 we conclude that the Neural A* planner
significantly reduces the number of nodes explored while still maintaining
sufficiently large path optimality. The results obtained from the
reimplementation were within 3.2% of the official implementation. From section
4.1.2 we also validate that provision of start and goal position to the
encoder plays a crucial role in guidance map generation. We also conclude that
ResNet-18 architecture leads to a drop in performance. For additional
experiments (Section 5), we have worked on the suggestions of the authors. The
generalizability of the model was displayed by the use of generic datasets and
coupling the planner with different A* implementations. We add UNet++ and GANs
implementations and experiment with dropout.

###  6.1 What was easy

Porting the code to Lightning was easy as the original source code is very
structured and well commented. Ablations reported in the paper are well
organized, and results are clearly reported making them easier to replicate.
Writing the code for additional experiments was also not difficult as the
original code is written well, allowing for fast execution.

###  6.2 What was difficult

Runtime calculations posed many difficulties due to differences in training
datasets, hardware and A* implementation. The paper had several parameters
that could be tuned hence it was difficult to ascertain their relative
importance. Experimentation was not possible on SDD Dataset which composed of
raw images due to training time constraints. Similarly, the TiledMP and CSM
datasets also took considerable time to train, and hence we only reproduced
the experiments from the original paper on these sets. All other experiments
were conducted on the MP dataset.

###  6.3 Communication with original authors

The conversation was based on the difficulties mentioned above along with
advice on how to approach the future improvements mentioned in the paper. The
authors’ replies were very prompt and thorough, enabling us to better
reproduce their work. We also implemented several of their suggestions, as
explained in Section 5.3.

##  7 Future Work

For future reproduction of this paper we suggest the use of "number of
explorations" as an additional metric. Exp and Hmean values alone can often be
misleading and difficult to analyze. On the other hand, the number of
explorations is an intuitive metric and makes it easier to compare models. The
author also suggested the use of number of explorations as a proxy metric that
is likely to be used for neural planners. We list down few experiments which
were planned but couldn’t be completed. We believe that future works can be
built on these two extensions:

  * •

Stochastic Path Sampling in Neural A* : For this, we believe that implementing
a replanning step coupled with dropout as suggested in [20] will be an
interesting addition to the planner. Similar to GANs, Variational Auto
Encoders can also be used. This would force the planner to propose multiple
paths for the same start and goal positions, thus overcoming its current
limitation.

  * •

Extension of Neural A* to high dimensions: For high dimensional spaces we
think that using an architecture like 3D UNet [27] or adding a feature
extractor similar to [28] can be a plausible solution.

## References

  * [1] R. Yonetani, T. Taniai, M. Barekatain, M. Nishimura, and A. Kanezaki, “Path planning using neural a* search,” in _International Conference on Machine Learning_. PMLR, 2021, pp. 12 029–12 039.
  * [2] W. Falcon et al., “Pytorch lightning,” _GitHub. Note: https://github.com/PyTorchLightning/pytorch-lightning_ , vol. 3, 2019.
  * [3] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio, “Generative adversarial nets,” _Advances in neural information processing systems_ , vol. 27, 2014.
  * [4] Z. Zhou, M. M. R. Siddiquee, N. Tajbakhsh, and J. Liang, “Unet++: A nested u-net architecture for medical image segmentation,” in _Deep learning in medical image analysis and multimodal learning for clinical decision support_. Springer, 2018, pp. 3–11.
  * [5] K. Simonyan and A. Zisserman, “Very deep convolutional networks for large-scale image recognition,” _arXiv preprint arXiv:1409.1556_ , 2014.
  * [6] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2016, pp. 770–778.
  * [7] A. K. Issar, K. Mali, A. Mehta, K. Uppal, S. Mishra, and D. Chakravarty, “Reproducibility of "fda: Fourier domain adaptation forsemantic segmentation,” 2021.
  * [8] O. Ronneberger, P. Fischer, and T. Brox, “U-net: Convolutional networks for biomedical image segmentation,” in _International Conference on Medical image computing and computer-assisted intervention_. Springer, 2015, pp. 234–241.
  * [9] I. Hubara, M. Courbariaux, D. Soudry, R. El-Yaniv, and Y. Bengio, “Binarized neural networks,” _Advances in neural information processing systems_ , vol. 29, 2016.
  * [10] L. Biewald, “Experiment tracking with weights and biases,” 2020, software available from wandb.com. [Online]. Available: https://www.wandb.com/
  * [11] M. Bhardwaj, S. Choudhury, and S. Scherer, “Learning heuristic search via imitation,” in _Conference on Robot Learning_. PMLR, 2017, pp. 271–280.
  * [12] N. R. Sturtevant, “Benchmarks for grid-based pathfinding,” _IEEE Transactions on Computational Intelligence and AI in Games_ , vol. 4, no. 2, pp. 144–148, 2012.
  * [13] A. Robicquet, A. Sadeghian, A. Alahi, and S. Savarese, “Learning social etiquette: Human trajectory understanding in crowded scenes,” in _European conference on computer vision_. Springer, 2016, pp. 549–565.
  * [14] S. Choudhury, M. Bhardwaj, S. Arora, A. Kapoor, G. Ranade, S. Scherer, and D. Dey, “Data-driven planning via imitation learning,” _The International Journal of Robotics Research_ , vol. 37, no. 13-14, pp. 1632–1672, 2018.
  * [15] M. Vlastelica, A. Paulus, V. Musil, G. Martius, and M. Rolínek, “Differentiation of blackbox combinatorial solvers,” _arXiv preprint arXiv:1912.02175_ , 2019.
  * [16] I. Pohl, “Heuristic search viewed as path finding in a graph,” _Artif. Intell._ , vol. 1, pp. 193–204, 1970.
  * [17] A. Tamar, Y. Wu, G. Thomas, S. Levine, and P. Abbeel, “Value iteration networks,” _arXiv preprint arXiv:1602.02867_ , 2016.
  * [18] L. Lee, E. Parisotto, D. S. Chaplot, E. Xing, and R. Salakhutdinov, “Gated path planning networks,” in _International Conference on Machine Learning_. PMLR, 2018, pp. 2947–2955.
  * [19] N. D. Ratliff, J. A. Bagnell, and M. A. Zinkevich, “Maximum margin planning,” in _Proceedings of the 23rd international conference on Machine learning_ , 2006, pp. 729–736.
  * [20] A. H. Qureshi, A. Simeonov, M. J. Bency, and M. C. Yip, “Motion planning networks,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 2118–2124.
  * [21] A. Gupta, J. Johnson, L. Fei-Fei, S. Savarese, and A. Alahi, “Social gan: Socially acceptable trajectories with generative adversarial networks,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 2255–2264.
  * [22] T. Salzmann, B. Ivanovic, P. Chakravarty, and M. Pavone, “Trajectron++: Dynamically-feasible trajectory forecasting with heterogeneous data,” in _European Conference on Computer Vision_. Springer, 2020, pp. 683–700.
  * [23] Z. Li, J. Wang, and M. Q.-H. Meng, “Efficient heuristic generation for robot path planning with recurrent generative model,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_ , 2021, pp. 7386–7392.
  * [24] W. Falcon and K. Cho, “A framework for contrastive self-supervised learning and designing a new approach,” _arXiv preprint arXiv:2009.00104_ , 2020.
  * [25] P. E. Hart, N. J. Nilsson, and B. Raphael, “A formal basis for the heuristic determination of minimum cost paths,” _IEEE Transactions on Systems Science and Cybernetics_ , vol. 4, no. 2, pp. 100–107, 1968.
  * [26] E. Rich and K. Knight, _Artificial Intelligence_ , 2nd ed. McGraw-Hill Higher Education, 1990.
  * [27] Ö. Çiçek, A. Abdulkadir, S. S. Lienkamp, T. Brox, and O. Ronneberger, “3d u-net: Learning dense volumetric segmentation from sparse annotation,” _ArXiv_ , vol. abs/1606.06650, 2016.
  * [28] P. Biasutti, V. Lepetit, J.-F. Aujol, M. Brédif, and A. Bugeau, “Lu-net: An efficient network for 3d lidar point cloud semantic segmentation based on end-to-end-learned 3d features and u-net,” 2019.