[1]Technical University of Berlin
[2]Science of Intelligence Cluster of Excellence, TUB
[3]Computer Eng. Dept., Bilkent University. This work was supported by TUBITAK under 2232 program with project number 121C148 (“LiRA“).

In this paper, we propose using deep neural architectures (i.e., vision transformers and ResNet) as heuristics for sequential decision-making in robotic manipulation problems.
This formulation enables predicting the subset of objects that are relevant for completing a task.
Such problems are often addressed by task and motion planning (TAMP) formulations combining symbolic reasoning and continuous motion planning.
In essence, the action-object relationships are resolved for discrete, symbolic decisions that are used to solve manipulation motions (e.g., via nonlinear trajectory optimization).
However, solving long-horizon tasks requires consideration of all possible action-object combinations which limits the scalability of TAMP approaches.
To overcome this combinatorial complexity, we introduce a visual perception module integrated with a TAMP-solver.

Given a task and an initial image of the scene, the learned model outputs the relevancy of objects to accomplish the task.
By incorporating the predictions of the model into a TAMP formulation as a heuristic, the size of the search space is significantly reduced.
Results show that our framework finds feasible solutions more efficiently when compared to a state-of-the-art TAMP solver.

arrow = [draw, -latex']
line = [draw, -]
input = [rectangle, draw, fill=white!20,
text width=5em, text centered, rounded corners, minimum height=1cm, minimum width=2cm]
cnn = [rectangle, draw, fill=white!20,
text width=1em, text centered, rounded corners, minimum height=2.5cm, minimum width=0.5cm]
resnetl = [rectangle, draw, fill=white!20,
text width=1em, text centered, rounded corners, minimum height=2cm, minimum width=0.5cm]
resnets = [rectangle, draw, fill=white!20,
text width=1em, text centered, rounded corners, minimum height=1cm, minimum width=0.5cm]
t = [rectangle, draw, fill=white!20,
text width=1em, text centered, rounded corners, minimum height=.5cm, minimum width=.5cm]
output_sigma = [circle, draw, fill=white!20,
text width=1em, text centered, rounded corners, minimum size=.2cm]

§ INTRODUCTION

Task and Motion Planning (TAMP) usually combines a discrete search on a symbolic, logical level with a geometric planner.
Discrete search is employed to find the high-level symbolic action sequence, which, in turn, informs the geometric planner to solve for a feasible motion path that satisfies the goal state.
Due to the huge computational overhead of solving geometric problems and the combinatorial complexity of discrete decisions, solving a TAMP problem often requires a large amount of computation [Wells et al., 2019, Toussaint et al., 2018, Rodrıguez et al., 2019, Driess et al., 2019].

The process of task planning often depends on the number of objects that are available in the environment and the operations that the robot can perform.
The algorithm builds the corresponding search tree based on the combination of objects and robot action set and expands the tree step by step (with action-object(s) tuples) until it reaches the goal state.

Exponential growth rate of the search tree tremendously diminishes the efficiency of the TAMP methods.
In reality, there are often multiple objects in the same environment, and a large number of operational possibilities by the autonomous agent, which impede the usability of TAMP solvers in such complex scenarios.
Fortunately though, the task can usually be achieved by only considering a small subset of these objects.

Machine learning methods using visual data offer efficient ways to solve robotic manipulation problems [Kase et al., 2020, Driess et al., 2021].
However, one challenge in integrating machine learning into TAMP is how to set up the problem (for example) by encoding objects in the scene as fixed-length inputs and then feeding them to the learning algorithm [Driess et al., 2020].
In general, learning methods are often integrated into robotic problems for a single task only, where fixed-size feature representation is often sufficient.
On the other hand, TAMP formulations are designed to tackle a broad set of tasks in different environments with varying numbers of objects and goals.
This hinders the usage of such fixed size features for solving sequential robotic manipulation problems.

In this study, we argue that task-relevant objects can be inferred directly from the sensor space of the robot, instead of in a feature space.
The sensor space has a fixed dimensionality, making it suitable for standard learning algorithms. Furthermore, the sensory information provides details about what the robot can extract from the current world state.
However, it does not account for historical data or other prior information.
Therefore, this paper aims to answer the following question: Is it is possible to learn the relevance of an object in the scene given the goal state?
This is critical because, as mentioned previously, the complexity of the TAMP algorithm belongs to the NP-hard class, and therefore the efficiency of the algorithm can be greatly improved if its search space can be reduced by learning methods.

[Scene from camera]

Test scenario with six objects: (a) shows the overview of this test scenario, (b) is the view from the camera located on the robot end-effector, (c, d and e) are the canonical views of objects, the first two are related to the goal predicate, and the third is the object whose relevance is queried. Take the predicate (on-left-side) as an example, using these canonical views as inputs would query the model whether object in (e) should be considered when the goal is to put object in (c) on the left-side of object in (d).

We propose a deep neural network that takes a scene image, canonical views of objects included in the goal predicate, and the canonical view of another object in the scene to predict the relevance of that object for achieving the task (Fig. <ref>).
The relevance of objects can then be used to improve the time complexity of mixed-integer programs, which is one way to realize TAMP, where the discrete action (integer) variable represents the abstract decision and the resulting nonlinear trajectory optimization problem represents the geometric part.

In the experiments, we consider the problem of grasping the target object with one robot arm and placing it on the position that is indicated by the goal predicate, as shown in Fig. <ref>.
An object can be considered as irrelevant due the non-interference of the trajectory of reaching the goal state. In contrast, if an object blocks the trajectory (i.e., there might be a collision), then the symbolic planner has to offer a different solution that instructs the robot to move this object out of the current place, which in turn makes this object relevant.

One of the experiment is shown in Fig. <ref>.
As an example, if our goal is to place object c to the left of object d, then the symbolic planner, without prior knowledge, will traverse all possibilities and eventually find a path that satisfies the “object c is to the left of object d” state, or if there is no possible solution, then it returns non-reachable.
However, as humans, we can quickly evaluate that no other object needs to be moved except for these two objects.
This is what we want to achieve with our trained classifier, thus greatly limiting the expansion of the search space.

In general, completing the goal requires actions on very few objects, and the majority of the objects that appear in the scene are irrelevant to the goal state.
Hence, we train several neural network architectures, including vision transformers (ViT) and residual neural networks (ResNet), to predict the relevancy of objects.
The learned classifier then guides the discrete search.

The main contributions are as follows:

* We propose using vision-based deep neural networks to predict the relevance of the objects from visual input for sequential manipulation planning.

* We propose a flexible and scalable predicate extractor to represent the spatial relationship.

* We develop the use of the network to guide the discrete search as a heuristic.

We support those contributions by showing that the trained models reach high accuracy, and our approach significantly reduces the number of motion planning problems in general, and thus improving computational efficiency of a state-of-the-art TAMP solver.
By evaluating multiple deep neural network models on the same dataset, we did not observe significant performance differences, e.g., using a ViT does not appear to be a better choice [Yuan et al., 2022].

context_patches = [rectangle, draw, fill=blue!20, inner sep=0.0mm,
text width=10em, text height = 1.8cm, text centered, rounded corners, minimum height=2cm, minimum width=11cm]
context_image = [rectangle, draw, fill=orange!20, inner sep=0.0mm,
text width=20em, text centered, rounded corners=2pt, minimum height=1.2cm, minimum width=10cm]
context_patch = [rectangle, draw, fill=orange!20, inner sep=0.0mm,
text width=0em, text centered, rounded corners=2pt, minimum height=0.5cm, minimum width=0.6cm]
canonical_views = [rectangle, draw, fill=red!20, inner sep=0.0mm,
text width=10em, text height = 1.8cm, text centered, rounded corners, minimum height=2cm, minimum width=5.1cm]
goal_pred = [rectangle, draw, fill=gray!15, inner sep=0.0mm,
text width=5em, text height = 1cm, text centered, rounded corners, minimum height=1.2cm, minimum width=2.3cm]
query = [rectangle, draw, fill=gray!5, inner sep=0.0mm,
text width=4em, text height = 1cm, text centered, rounded corners, minimum height=1.2cm, minimum width=0.8cm]
resnet_large = [rectangle, draw, fill=orange!20, inner sep=0.0mm,
text width=20em, text centered, rounded corners, minimum height=1cm, minimum width=11cm]
resnet_small = [rectangle, draw, fill=orange!20, inner sep=0.0mm,
text width=10em, text centered, rounded corners, minimum height=1cm, minimum width=5.1cm]
projection_patches = [rectangle, draw, fill=orange!20, inner sep=0.0mm,
text width=20em, text centered, rounded corners, minimum height=1cm, minimum width=16.5cm]
projection_patch = [rectangle, draw, fill=orange!20, inner sep=0.0mm,
text width=0em, text centered, rounded corners=2pt, minimum height=0.5cm, minimum width=0.3cm]
projection_patch_description = [rectangle, draw, fill=red!20, inner sep=0.0mm,
text width=0.2cm, text centered, rounded corners=2pt, minimum height=0.5cm, minimum width=0.3cm]
transformer = [rectangle, draw, fill=gray!20, inner sep=0.0mm,
text width=20em, text centered, rounded corners, minimum height=1cm, minimum width=16.5cm]
embedding = [rectangle, draw, fill=purple!20, inner sep=0.0mm,
text width=0.2cm, text centered, rounded corners=2pt, minimum height=0.809cm, minimum width=0.3cm]
readout = [rectangle, draw, fill=green!10, dashed,
text width=20em, text depth = 2.6cm, anchor=north, text centered, rounded corners, minimum height=2.6cm, minimum width=16.5cm]
dataset = [rectangle, draw, fill=lime!10, text=black, dashed,
text width=20em, text depth = 2.6cm, anchor=north, text centered, rounded corners, minimum height=2cm, minimum width=16.5cm]
mlp = [rectangle, draw, fill=white!20,
text width=5em, text centered, rounded corners, minimum height=1cm, minimum width=3.4cm]
output = [rectangle, draw, fill=lime!20,
text width=5em, text centered, rounded corners, minimum height=1cm, minimum width=2cm]
arrow_txt = [text width=2cm,midway,right=0em,align=left, font=8pt12]

§ RELATED WORK

§.§ Task and Motion Planning

A TAMP problem is typically specified in two parts: the first part is a symbolic discrete decision, and the second part is a continuous motion planner.
The symbolic discrete decision proposes a discrete node to reach the goal state, while the continuous motion planner is intended to solve the feasibility problem between discrete states.
There are three mainstream solutions for a motion planner: sampling based [Alili et al., 2010, Srivastava et al., 2014, de Silva et al., 2013, Dantam et al., 2018, Siméon et al., 2004], discretization based on the configuration/action space [Erdem et al., 2011, Lagriffoul et al., 2012, Lagriffoul et al., 2014, Lozano-Pérez and Kaelbling, 2014], and the optimization based methods [Shoukry et al., 2016, Hadfield-Menell et al., 2016].
Logic Geometric Programming (LGP) [Toussaint, 2015] also falls into the optimization-based category, where the logic imposes a skeleton of active constraints on a nonlinear program [Toussaint et al., 2020, Ha et al., 2020].
Formalizing this hybrid nature of the contact and interaction modes in robotics  results in a mixed-integer program [Hogan et al., 2018, Deits and Tedrake, 2014], and LGP is formulated as such for TAMP problems [Toussaint et al., 2018].

A major problem of LGP, though, is the combinatorial complexity that is introduced by the discrete variables (i.e., due to the number of actions or objects in the scene).
In this study, we propose to address this complexity by reducing the search space using a vision-based deep learning model that predicts a small set of objects to be interacted with to realize a task.

§.§ Learning for Task and Motion Planning

Advances in deep networks have made it possible to build a planning system that is based on sensor inputs.
For example, attention mechanisms offer unique opportunities to learn representations of the environment [Vaswani et al., 2017, Dosovitskiy et al., 2020, Carion et al., 2020, Guo et al., 2022].
Such models can be employed to reason about the spatial relationship between objects, thus the symbolic planner can leverage this information for high level planning [Yuan et al., 2022, Nguyen et al., 2020].
Our model is based on vision transformer and residual neural network [Vaswani et al., 2017, He et al., 2016], and uses supervised learning to give the model the ability to infer task-relevant objects based on an image in the scene.
The aim is to use visual inference to guide the symbolic planner to consider the option that is more likely to succeed in a given situation.
The viability of image-based deep networks for direct service to planning and control systems has also been confirmed by many studies [Kase et al., 2020, Driess et al., 2020, Driess et al., 2021].

In spirit, the motivation with SORNET [Yuan et al., 2022] is most similar to ours.
The input to SORNET comprises image patches and canonical views of the objects to be queried.
The output is the spatial relationship between the two objects queried in the scene.

However, this information is often not directly used in task and motion planning, because knowing the spatial relationship between objects does not help much in motion planning.
SORNET uses the embedding output of the vision transformer (ViT) to train additional networks that models the predicates.
We, instead, use ViT to directly train a classifier that outputs the probability of relevancy of objects for a given goal.
This allows us to directly integrate the model output into our TAMP formulation.

In general, the method provided in this study inferred the task-related objects on the basis of images, which can limit the space of the discrete search part of task and motion planning, thus improving the overall search efficiency.

§ MIXED-INTEGER PROGRAM FOR ROBOT TRAJECTORY PLANNING

In this section, we explain our task and motion planning formulation based on  [Toussaint, 2015, Toussaint and Lopes, 2017].
Use cases for this formulation are two-fold: (i) we solve a large number of complex rearrangement problems for obtaining labels of task-relevancy of objects,
(ii) we integrate our learned classifiers for improving the computational efficiency of solving such sequential manipulation problems through the trajectory optimization framework that we utilize to create labels as training data and supervise deep learning as it solves the existing problem.

Consider $s$ a symbolic state variable drawn from the set $S$.
This set $S$ characterizes the scene, which is formulated by $\mathcal{X} \subset \mathbb{R}^{n(S)} \times SE(3)^{m(s, S)} \times \mathbb{R}^{6 \cdot n_{\textrm{cp}}(s, S)}$.
Here, the scene configuration integrates robot joint positions ($n(S)$-dimensional), poses of $m(s, S)$ rigid objects, and the dynamics of wrench interactions at the contact points ($n_{\textrm{cp}}(s, S)$).

To identify a globally consistent and feasible path $x: [0, KT] \rightarrow \mathcal{X}(s_{k(t)}, S)$ within this configuration, our objective is to minimize

\begin{align}
P(g, S) = \min_{\substack{ K \in \mathbb{N} \\ x: [0, KT] \rightarrow \mathcal{X} \\ a_{1:K}, s_{1:K}}} & \int_{0}^{KT} c( x(t), \dot{x}(t), \ddot{x}(t), s_{k(t)}, S)dt \label{eq:1a} \\
s.t. \qquad & \nonumber \\
\forall_{t \in [0, KT]} : \ & h_{\textrm{eq}}(x(t), \dot{x}(t), s_{k(t)}, S) = 0 \label{eq:1b} \\
\forall_{t \in [0, KT]} : \ & h_{\textrm{ineq}}(x(t), \dot{x}(t), s_{k(t)}, S) \leq 0 \label{eq:1c}\\
\forall_{k=1,...,K} : \ & h_{\textrm{sw}}(x(kT), \dot{x}(kT), a_{k}, S) = 0 \label{eq:1d}\\
\forall_{k=1,...,K} : \ & a_{k} \in \mathbb{A}(s_{k-1}, S) \label{eq:1e}\\
\forall_{k=1,...,K} : \ & s_{k} = succ(s_{k-1}, a_{k}) \label{eq:1f}\\
& x(0) = \tilde{x}_{0}(S) \label{eq:1g}\\
& s_{0} = \tilde{s}_{0}(S) \label{eq:1h}\\
& s_{K} \in \mathcal{S}_{\textrm{goal}}(g). \label{eq:1i}
\end{align}

There are $K \in \mathbb{N}$ phases, each having a uniform duration $T > 0$.
We assume global continuity for the path $(x~\in~C([0,~TK]))$, and smoothness for each phase $x~\in~C^{2}([(k-1)T,~kT])$.
Note that the value of $K$ is part of the decision problem itself.

These phases are often referred to as kinematic modes, as mentioned in [Toussaint et al., 2018, Mason, 1985].
Notably, as the number of degrees of freedom for the objects and the number of contact interactions vary in accordance with the symbolic state $s_{k(t)}$, where $k(t)=\left\lfloor{t/T}\right\rfloor$, the dimensionality of the path can also fluctuate between different stages.

Equations (<ref>), (<ref>), and (<ref>) collectively define the waypoint in phase $k$ of the entire path.
Here, $s_{k} \in \mathcal{S}(S)$ represents the state within the symbolic domain.
Switching constraints between various kinematic modes are characterized by  (<ref>).
The end of each path phase is marked by an action transition  (<ref>).
It is important to note that we make the assumption that $c, h_{eq}$ and $h_{ineq}$ are differentiable for a fixed $s$.

The discrete action $a_{k}$ (<ref>) paired with an object is called an action-object tuple.
During tree search to find the action sequence, LGP expands the current node based on a possible action-object tuple.
Once the scene becomes complex, the number of possible combinations of action-objects grows exponentially.
This makes traversing the whole search space infeasible.
Adding reasonable prior knowledge then becomes critical to prune irrelevant nodes.

Given an initial symbolic state $s_{0}=\tilde{s}_{0}(S)$ (<ref>), the action sequence $a_{1:K}$ establishes the sequence of symbolic states $s_{0:K}$. The goal state is symbolically defined by the set $\mathcal{S}_{\textrm{goal}}(g)$, where the symbolic goal is represented by $g \in \mathbb{G}(S)$.

A solution to the LGP is represented by an action sequence, denoted as $a_{1:K}$, which guides the symbolic states from their initial configuration towards the symbolic goal state. Additionally, it is essential to ensure that the constraints within and between each stage are satisfied.

We characterize the feasibility of an action sequence $a_{1:K}=(a_{1},...,a_{K})$ through the presence of a corresponding path,

\begin{equation} \label{eq:2}
F_{S}(a_{1:K}) =
\begin{cases}
1,  & \text{if } \exists x : [0, KT] \rightarrow \mathcal{X} : \eqref{eq:1b} - \eqref{eq:1h}\\
0,  & \text{otherwise}
\end{cases}.
\end{equation}

We also define a function that describes whether an object is involved in solving the task,

\begin{equation} \label{eq:3}
Z_{S}(o) =
\begin{cases}
1,  & \text{if } F_{S}(a_{1:K}) = 1 \rightarrow \exists a_{k \in K}: \theta(a_{k}) = o\\
0,  & \text{otherwise}
\end{cases}.
\end{equation}

$\theta(a)$ takes an action as parameter and retrieves the object involved in the action. All actions have only one operation object, thus, the outcome of $\theta$ is a single object.

The comprehensive LGP formulation, as given in (<ref>), aims to identify not only a viable solution but also a sequence of actions that results in the lowest trajectory costs, as expressed in (<ref>), when compared to all alternative sequences leading to the goal.
A viable solution is referred to as a solution to the TAMP/LGP problem, while a solution that minimizes the cost is termed an optimal solution of the LGP.
Given the multitude of possible discrete action sequences, rather than the optimal solutions, our primary focus is on acquiring an optimal solution.

Equation  (<ref>) represents the feasibility of the nonlinear program (<ref>) mathematically.
In practice, though,  (<ref>) is solved by discretizing $x$ and using a Gauss-Newton based augmented-Lagrangian method.

Specifically, the action sequence $F_{S}(a_{1:K})$ is constructed such that violations $(h_{eq} \neq 0, h_{ineq} > 0)$ are allowed below a predefined threshold for each waypoint.
If these violations exceed the threshold, the problem is deemed infeasible.

[node distance = 0cm, auto]
(scene_img)[inner sep=0pt] ;
(scene_img_caption) [above=0cm of scene_img] scene $S$;

[node distance = 0cm, auto]
(scene_img) [input] ;
(scene_img_caption) [above=-0.1cm of scene_img] input scene image;

(canonical_view) [input, below=0.5cm of scene_img, fill=red!20]

(canonical_view_caption) [above=-0.1cm of canonical_view] canonical view;

(scene_img_patches) [input, right=0.5cm of scene_img, fill=blue!20]
[width=.25,trim=.0 .66 .66 .0,clip]img/scene_img.png
[width=.25,trim=.33 .66 .33 .0,clip]img/scene_img.png
[width=.25,trim=.66 .66 .0 .0,clip]img/scene_img.png

[width=.25,trim=.0 .33 .66 .33,clip]img/scene_img.png
[width=.25,trim=.33 .33 .33 .33,clip]img/scene_img.png
[width=.25,trim=.66 .33 .0 .33,clip]img/scene_img.png

[width=.25,trim=.0 .0 .66 .66,clip]img/scene_img.png
[width=.25,trim=.33 .0 .33 .66,clip]img/scene_img.png
[width=.25,trim=.66 .0 .0 .66,clip]img/scene_img.png
(scene_img_patches_caption) [above=-0.1cm of scene_img_patches] scene image patches;

(cnn) [cnn, right=0.5cm of scene_img_patches, yshift=-0.55cm, fill=orange!20] 90CNN;
(label) [cnn, right=0.5cm of cnn] 90Attach position label;
(transformer) [cnn, right=0.5cm of label, fill=gray!20] 90Transformer;
(t1) [t, right=0.5cm of transformer, yshift=0.6cm, fill=red!50] 90;
(t2) [t, right=0.5cm of transformer, yshift=0cm, fill=green!50] 90;
(t3) [t, right=0.5cm of transformer, yshift=-0.6cm, fill=cyan!50] 90;
(embedding_caption) [above=0.01cm of t1] obj embedding;
(readout) [cnn, right=0.5cm of t2, yshift=0cm, fill=green!10] 90Predicate Extractor;
(output_sigma) [output_sigma, right=0.5cm of readout] $\sigma$;
(output_sigma_caption) [right=0.1cm of output_sigma] $t(I, p_{0}, p_{1}, q)$;

[arrow] (scene_img.east) – (scene_img_patches.west);
[arrow] (scene_img_patches.east) – (scene_img_patches-|cnn.west);
[arrow] (canonical_view.east) – (canonical_view-|cnn.west);
[arrow] (cnn.east) – (label.west);
[arrow] (label.east) – (transformer.west);
[arrow] ([yshift=0.6cm,xshift=0cm]transformer.east) – (t1.west);
[arrow] ([yshift=0cm,xshift=0cm]transformer.east) – (t2.west);
[arrow] ([yshift=-0.6cm,xshift=0cm]transformer.east) – (t3.west);

[arrow] (t1.east) – ([yshift=0.6cm,xshift=0cm]readout.west);
[arrow] (t2.east) – ([yshift=0cm,xshift=0cm]readout.west);
[arrow] (t3.east) – ([yshift=-0.6cm,xshift=0cm]readout.west);

[arrow] (readout.east) – (output_sigma.west);

Proposed neural net architecture based on transformer to predict the involvement of objects in the physical interaction for a given goal predicate.

[node distance = 0cm, auto]

(cps) [context_patches] Context Patches;
(cvs) [canonical_views, right=0.4cm of cps] Canonical Object Views;
(gp) [goal_pred, above=-1.28cm of cvs, xshift=-0.85cm] Predicate;
(q) [query, right=0.4cm of gp] Query;
(pps) [projection_patches, above=0.3cm of cps, xshift=2.75cm] Linear Projection of Flattened Patches;

(cp1) [context_patch, above=0cm of pps, xshift=-7.15cm, yshift=-2cm] ;
(cp2) [context_patch, right=0.5cm of cp1] ;
(cp3) [context_patch, right=0.5cm of cp2] ;
(cp4) [context_patch, right=0.5cm of cp3] ;
(cp5) [context_patch, right=0.5cm of cp4] ;
(cp6) [context_patch, right=0.5cm of cp5] ;
(cp7) [context_patch, right=0.5cm of cp6] ;
(cp8) [context_patch, right=0.5cm of cp7] ;
(cp9) [context_patch, right=0.5cm of cp8] ;

(qp1) [context_patch, above=0cm of pps, xshift=4.3cm, yshift=-2cm] ;
(qp2) [context_patch, right=0.5cm of qp1] ;
(qp3) [context_patch, right=1.1cm of qp2] ;

(pp1) [projection_patch, above=0.2cm of pps, xshift=-7cm] ;
(ppd1) [projection_patch_description, left=0cm of pp1] 1;

(pp2) [projection_patch, right=0.8cm of pp1] ;
(ppd2) [projection_patch_description, left=0cm of pp2] 2;

(pp3) [projection_patch, right=0.8cm of pp2] ;
(ppd3) [projection_patch_description, left=0cm of pp3] 3;

(pp4) [projection_patch, right=0.8cm of pp3] ;
(ppd4) [projection_patch_description, left=0cm of pp4] 4;

(pp5) [projection_patch, right=0.8cm of pp4] ;
(ppd5) [projection_patch_description, left=0cm of pp5] 5;

(pp6) [projection_patch, right=0.8cm of pp5] ;
(ppd6) [projection_patch_description, left=0cm of pp6] 6;

(pp7) [projection_patch, right=0.8cm of pp6] ;
(ppd7) [projection_patch_description, left=0cm of pp7] 7;

(pp8) [projection_patch, right=0.8cm of pp7] ;
(ppd8) [projection_patch_description, left=0cm of pp8] 8;

(pp9) [projection_patch, right=0.8cm of pp8] ;
(ppd9) [projection_patch_description, left=0cm of pp9] 9;

(q1) [projection_patch, above=0.2cm of pps, xshift=4.45cm] ;
(qm1) [projection_patch_description, left=0cm of q1] 0;

(q2) [projection_patch, right=0.8cm of q1] ;
(qm2) [projection_patch_description, left=0cm of q2] 0;

(q3) [projection_patch, right=1.4cm of q2] ;
(qm3) [projection_patch_description, left=0cm of q3] 0;

(transformer) [transformer, above=1cm of pps] Transformer;

(eb1) [embedding, above=0.2cm of transformer, xshift=4.45cm, fill=red!50] ;
(eb2) [embedding, right=0.8cm of eb1, fill=green!50] ;
(eb3) [embedding, right=1.4cm of eb2, fill=cyan!50] ;

(readout) [readout, above=1.2cm of transformer] PEs;
(mlp1) [mlp, above=2.5cm of transformer, xshift=-5.9cm, fill=red!25] PE1;
(mlp2) [mlp, right=0.5cm of mlp1, fill=teal!25] PE2;
(mlp3) [mlp, right=0.5cm of mlp2, fill=blue!25] PE3;
(mlp4) [mlp, right=0.5cm of mlp3, fill=brown!25] PE4;

(output1) [output, above=1.5cm of mlp1, fill=red!50] Left;
(output2) [output, above=1.5cm of mlp2, fill=teal!50] Right;
(output3) [output, above=1.5cm of mlp3, fill=blue!50] ...;
(output4) [output, above=1.5cm of mlp4, fill=brown!50] ...;

[arrow] (cp1.north) – (cp1|-pps.south);
[arrow] (cp2.north) – (cp2|-pps.south);
[arrow] (cp3.north) – (cp3|-pps.south);
[arrow] (cp4.north) – (cp4|-pps.south);
[arrow] (cp5.north) – (cp5|-pps.south);
[arrow] (cp6.north) – (cp6|-pps.south);
[arrow] (cp7.north) – (cp7|-pps.south);
[arrow] (cp8.north) – (cp8|-pps.south);
[arrow] (cp9.north) – (cp9|-pps.south);

[arrow] (qp1.north) – (qp1|-pps.south);
[arrow] (qp2.north) – (qp2|-pps.south);
[arrow] (qp3.north) – (qp3|-pps.south);

[line] (pps.north-|pp1) – (pp1.south);
[line] (pps.north-|pp2) – (pp2.south);
[line] (pps.north-|pp3) – (pp3.south);
[line] (pps.north-|pp4) – (pp4.south);
[line] (pps.north-|pp5) – (pp5.south);
[line] (pps.north-|pp6) – (pp6.south);
[line] (pps.north-|pp7) – (pp7.south);
[line] (pps.north-|pp8) – (pp8.south);
[line] (pps.north-|pp9) – (pp9.south);

[line] (pps.north-|q1) – (q1.south);
[line] (pps.north-|q2) – (q2.south);
[line] (pps.north-|q3) – (q3.south);

[arrow] (pp1.north) – (pp1|-transformer.south);
[arrow] (pp2.north) – (pp2|-transformer.south);
[arrow] (pp3.north) – (pp3|-transformer.south);
[arrow] (pp4.north) – (pp4|-transformer.south);
[arrow] (pp5.north) – (pp5|-transformer.south);
[arrow] (pp6.north) – (pp6|-transformer.south);
[arrow] (pp7.north) – (pp7|-transformer.south);
[arrow] (pp8.north) – (pp8|-transformer.south);
[arrow] (pp9.north) – (pp9|-transformer.south);

[arrow] (q1.north) – (q1|-transformer.south);
[arrow] (q2.north) – (q2|-transformer.south);
[arrow] (q3.north) – (q3|-transformer.south);

[line] (transformer.north-|eb1) – (eb1.south);
[line] (transformer.north-|eb2) – (eb2.south);
[line] (transformer.north-|eb3) – (eb3.south);

[arrow, line width=0.4mm, color=red!50] (eb1.north) – ++(0,0.4) -| ([yshift=0cm,xshift=-1cm]mlp1.south);
[arrow, line width=0.4mm, color=red!50] (eb1.north) – ++(0,0.4) -| ([yshift=0cm,xshift=-1cm]mlp2.south);
[arrow, line width=0.4mm, color=red!50] (eb1.north) – ++(0,0.4) -| ([yshift=0cm,xshift=-1cm]mlp3.south);
[arrow, line width=0.4mm, color=red!50] (eb1.north) – ++(0,0.4) -| ([yshift=0cm,xshift=-1cm]mlp4.south);

[arrow, line width=0.4mm, color=teal] (eb2.north) – ++(0,0.7) -| ([yshift=0cm,xshift=0cm]mlp1.south);
[arrow, line width=0.4mm, color=teal] (eb2.north) – ++(0,0.7) -| ([yshift=0cm,xshift=0cm]mlp2.south);
[arrow, line width=0.4mm, color=teal] (eb2.north) – ++(0,0.7) -| ([yshift=0cm,xshift=0cm]mlp3.south);
[arrow, line width=0.4mm, color=teal] (eb2.north) – ++(0,0.7) -| ([yshift=0cm,xshift=0cm]mlp4.south);

[arrow, line width=0.4mm, color=cyan] (eb3.north) – ++(0,1.0) -| ([yshift=0cm,xshift=1cm]mlp1.south);
[arrow, line width=0.4mm, color=cyan] (eb3.north) – ++(0,1.0) -| ([yshift=0cm,xshift=1cm]mlp2.south);
[arrow, line width=0.4mm, color=cyan] (eb3.north) – ++(0,1.0) -| ([yshift=0cm,xshift=1cm]mlp3.south);
[arrow, line width=0.4mm, color=cyan] (eb3.north) – ++(0,1.0) -| ([yshift=0cm,xshift=1cm]mlp4.south);

[arrow, line width=0.2mm] (mlp1.north) – ([yshift=0cm,xshift=0cm]output1.south);
[arrow, line width=0.2mm] (mlp2.north) – ([yshift=0cm,xshift=0cm]output2.south);
[arrow, line width=0.2mm] (mlp3.north) – ([yshift=0cm,xshift=0cm]output3.south);
[arrow, line width=0.2mm] (mlp4.north) – ([yshift=0cm,xshift=0cm]output4.south);

The model structure. The input to the model consists of two major parts: (i) Context Patches: the cropped regions of the scene image, and (ii) Canonical Object Views: the descriptive image of the object to be queried. In the middle is the default Transformer model. Finally, there are four independent PEs, which correspond to different Predicates.

§ VISUAL DEEP LEARNING MODELS AS PLANNING HEURISTICS

Here, we describe the deep learning models that we have developed.
These models accept images as input and output information that assist TAMP to improve planning time.

§.§ Embedding Model

§.§.§ The Vision Transformer (ViT)

The inputs to the model are cropped image patches.

These patches are of two types: the first type is the scene image and the second one is the canonical views of the objects that are queried (Fig. <ref>).
The patches are first converted into an embedding, e.g., by a convolutional network.
Then, the corresponding positional information is concatenated, which is later adopted in the transformer module.

In the default ViT architecture [Dosovitskiy et al., 2020], the embedding at position zero is critical.
It is not derived from the input patch, but instead can be learned and then used for the classification.
This part is where our model differs from the original ViT.
Our implementation extends the use of location labels to enable the model to perform the query task in a given situation, i.e., to determine whether the query object needs to be involved in the process of completing the task given the target state.

The location information and the embedding are then passed into the transformer.
The transformer here is identical to the original transformer structure used for natural language processing [Vaswani et al., 2017].
However, unlike the traditional transformer, the vision transformer discards all the final embeddings of the patches from the scene image, and uses the final embedding of the zero position (Fig. <ref>).
A multilayer perceptron is finally implemented to use this embedding to predict the relevancy of the queried object for several predicates - we call this component the Predicate Extractor (PE).

§.§.§ The Residual Neural Network (ResNet)

ResNets extract features by convolutional and residual blocks [He et al., 2016].

We built two versions (Fig. <ref>):

For the first one ((D)efault-ResNet), two different sizes of ResNets are used to analyze the scene image and the canonical views separately.

A fully linked network then produces three embeddings similar to the ViT-based model.
The PE eventually receives the generated embedding data to extract the relevance information for the query object given a task described by canonical object views.
In the second one, which we named three-dimensional ResNet, we first segment the scene image into patches as well, and stack them with the canonical view.
Then, a residual neural network based on a 3D convolutional algorithm is used to extract the features (Fig. <ref>).

The features output from the residual-based neural network model can also be considered as embedded information and utilized by the Predicate Extractor (PE).

[node distance = 0cm, auto]
(scene_img) [input] ;
(scene_img_caption) [above=-0.1cm of scene_img] input scene image;

(canonical_view) [input, below=0.5cm of scene_img, fill=red!20]

(canonical_view_caption) [above=-0.1cm of canonical_view] canonical view;

(resnetl) [resnetl, right=0.5cm of scene_img, fill=orange!20] 90ResNet;
(resnets) [resnets, right=0.5cm of canonical_view, fill=orange!20] 90ResNet;
(transformer) [cnn, right=0.5cm of resnetl, yshift=-0.6cm, fill=gray!20] 90Dense layer;
(t1) [t, right=0.5cm of transformer, yshift=0.6cm, fill=red!50] 90;
(t2) [t, right=0.5cm of transformer, yshift=0cm, fill=green!50] 90;
(t3) [t, right=0.5cm of transformer, yshift=-0.6cm, fill=cyan!50] 90;
(embedding_caption) [above=0.01cm of t1] obj embedding;
(readout) [cnn, right=0.5cm of t2, yshift=0cm, fill=green!10] 90Predicate Extractor;
(output_sigma) [output_sigma, right=0.2cm of readout] $\sigma$;

[arrow] (scene_img.east) – ([yshift=0cm,xshift=0cm]resnetl.west);
[arrow] (canonical_view.east) – ([yshift=0cm,xshift=0cm]resnets.west);
[arrow] (resnetl.east) – (resnetl-|transformer.west);
[arrow] (resnets.east) – (resnets-|transformer.west);
[arrow] ([yshift=0.6cm,xshift=0cm]transformer.east) – (t1.west);
[arrow] ([yshift=0cm,xshift=0cm]transformer.east) – (t2.west);
[arrow] ([yshift=-0.6cm,xshift=0cm]transformer.east) – (t3.west);

[arrow] (t1.east) – ([yshift=0.6cm,xshift=0cm]readout.west);
[arrow] (t2.east) – ([yshift=0cm,xshift=0cm]readout.west);
[arrow] (t3.east) – ([yshift=-0.6cm,xshift=0cm]readout.west);

[arrow] (readout.east) – (output_sigma.west);

[node distance = 0cm, auto]
(scene_img) [input] ;
(scene_img_caption) [above=-0.1cm of scene_img] input scene image;

(canonical_view) [input, below=0.5cm of scene_img, fill=red!20]

(canonical_view_caption) [above=-0.1cm of canonical_view] canonical view;

(scene_img_patches) [input, right=0.5cm of scene_img, fill=blue!20]
[width=.25,trim=.0 .66 .66 .0,clip]img/scene_img.png
[width=.25,trim=.33 .66 .33 .0,clip]img/scene_img.png
[width=.25,trim=.66 .66 .0 .0,clip]img/scene_img.png

[width=.25,trim=.0 .33 .66 .33,clip]img/scene_img.png
[width=.25,trim=.33 .33 .33 .33,clip]img/scene_img.png
[width=.25,trim=.66 .33 .0 .33,clip]img/scene_img.png

[width=.25,trim=.0 .0 .66 .66,clip]img/scene_img.png
[width=.25,trim=.33 .0 .33 .66,clip]img/scene_img.png
[width=.25,trim=.66 .0 .0 .66,clip]img/scene_img.png
(scene_img_patches_caption) [above=-0.1cm of scene_img_patches] scene image patches;
(stacking) [resnets, right=0.5cm of canonical_view, fill=orange!20] 90Stacking;

(resnet) [resnets, right=0.5cm of stacking, fill=orange!20] 903D-ResNet;
(transformer) [cnn, right=0.5cm of scene_img_patches, yshift=-1.2cm, fill=gray!20] 90Dense layer;
(t1) [t, right=0.5cm of transformer, yshift=0.6cm, fill=red!50] 90;
(t2) [t, right=0.5cm of transformer, yshift=-0cm, fill=green!50] 90;
(t3) [t, right=0.5cm of transformer, yshift=-0.6cm, fill=cyan!50] 90;
(embedding_caption) [above=0.01cm of t1] obj embedding;
(readout) [cnn, right=0.5cm of t2, yshift=0cm, fill=green!10] 90Predicate Extractor;
(output_sigma) [output_sigma, right=0.2cm of readout] $\sigma$;

[arrow] (scene_img.east) – ([yshift=0cm,xshift=0cm]scene_img_patches.west);
[arrow] (canonical_view.east) – ([yshift=0cm,xshift=0cm]stacking.west);
[arrow] (scene_img_patches.south -| stacking.north) – (stacking.north);
[arrow] (stacking.east) – (resnet.west);
[arrow] (resnet.east) – (resnet-|transformer.west);
[arrow] (transformer.east|-t1.west) – (t1.west);
[arrow] (transformer.east|-t2.west) – (t2.west);
[arrow] (transformer.east|-t3.west) – (t3.west);

[arrow] (t1.east) – (t1-|readout.west);
[arrow] (t2.east) – (t2-|readout.west);
[arrow] (t3.east) – (t3-|readout.west);

[arrow] (readout.east) – (output_sigma.west);

Proposed architectures based on ResNet: (left) a model based on a conventional ResNet, and (right) a model based on a 3D-ResNet.

§.§ Embedding Information Extraction

One option is to model PE to have the same output dimensionality as the number of the predicates that we are concerned with: for example, the predicates on-the-left, on-the-right, in-front-of and behind.

In our current architecture, there is a single PE for each predicate as shown in Fig. <ref>.

We trained the model with batch polling method.

This method is based on a polling approach (Fig. <ref>), where the data corresponding to different PEs are taken out in batches and used for the corresponding training, thus updating the corresponding PEs and the transformer they share.
With the polling method, the transformer is updated with data from multiple predicates simultaneously at each gradient step, leading to a more reliable optimization.

§.§ Details of the Models

The goal of this study is to learn a classifier $t$ predicting an object set ${O} \subset \mathcal{O}$ that comprises the objects manipulated by the actions $a_{1:K}$ generated by a feasible nonlinear program $P(g, S)$.
In order to achieve this objective, the key consideration is how to represent the scene $S$ and the desired configuration $\mathcal{S}_{goal}(g)$ as inputs to the classifier. This representation should not only facilitate the training of a precise classifier but also allow it to perform well on different scenes, including those with varying numbers of objects.
Due to the fixed-size input characteristics of the image sensor, we can argue that the image sensor meets our requirements.
Formally, we assume that there is a generative process to produce an image $I \in \mathbb{R}^{w \times h \times 3}$ from the scene $S$, either via a rendering engine in simulation or a camera in the real world.
We choose an image acquired by a tilted (virtual) camera to capture as much geometric information as possible.

Thanks to the attention mechanism, all we need to identify for our query is a canonical view per object $q \in \mathbb{R}^{32 \times 32}$ obtained in advance.
Note that the dimensions of the canonical views here are same as the patches of scene image, which are also 32 by 32.

The input to the classifier $t(I, p_{0,...,L}, q_i)$ is the image $I$ as three channel image $I \in \mathbb{R}^{w \times h \times 3}$ and canonical views for the objects, which, together with the goal predicate, describe the goal and the queried object, $p_{0,...,L} \in \mathbb{R}^{32 \times 32}$, $q_i$, respectively ($L=1$ for our PEs in this study).
Fig. <ref> visualizes these input images, and Fig. <ref> describes our architecture with a high level of abstraction.

The network is trained to output the probability that the queried objects will be included in the manipulation plan which is the result of the nonlinear program $P(g, S)$

\begin{equation} \label{eq:4}
t(I, p_{0,...,L}, q_{i})=p(Z_{S}(i)|F_{S}(a_{1:K})=1)_{o \in \mathcal{O}}
\end{equation}

using the standard weighted binary cross-entropy loss

\begin{equation} \label{eq:5}
\begin{aligned}
L(w) = \sum_{i \in \mathcal{O}} & \eta Z_{S}(i)\log(t(I, p_{0,...,L}, q_{i};w)) \\
& + (1 - Z_{S}(i))\log(1 - t(I, p_{0,...L}, q_{i};w)).
\end{aligned}
\end{equation}

The training data $\mathcal{D}={(I_{i}, \mathrm{P}_{i}, Q_{i})}_{i=1}^{d}$ consists of scene images $I_{i}$ with the canonical views of the predicate related objects $\mathrm{P}_{i}$ and query object $Q_{i}$.
In our experiments the majority of the objects are irrelevant to the goal state. Therefore, adjusting the value of weighting factor $\eta \leq 1$ in eq. <ref> becomes important to balance the loss, otherwise the classifier could indiscriminately predict the biased label to achieve high accuracies.
If the prediction $t(I, p_{0,...L}, q_{i})$ for object $i$ is lower than the threshold $\beta \in \mathbb{R}$, then we decide that $Z_{S}(i) = 1$, i.e., the object is relevant to the action sequence.
In the experiments, we discuss the influence of this threshold.

For the ViT and the 3D-ResNet models, the scene image is first cropped into small patches of a predefined size, we use 32 by 32 in our experiments.
For the ViT-based model, the patches of the scene images, along with the canonical views associated with the predicates and the query objects, are then fed into the first layer of the CNN network, which is used to transform them into an embedding with a size of 768 required by the transformer.
The embedding converted from the scene image will then be attached with the corresponding position label, while the embedding generated from the canonical views will only be attached with position label zero.
In essence, we only use the embeddings that capture the correlation between context patches and canonical views, and also between canonical views themselves, while discarding the correlation between context patches (for details, see [Dosovitskiy et al., 2020, Yuan et al., 2022]).
These embeddings with labels attached are then fed into the standard transformer.
Finally, we take out the output corresponding to the canonical views and feed them into the final fully connected layer, which has two layers and 512 neurons per layer (which we call the Predicate Extractor).

§ GUIDING MOTION PLANNING USING LEARNED
VISUAL HEURISTICS

The classifier defined in (<ref>) can be used not only to predict the relevance of a goal state depicted by a goal predicate within a scene, but also to serve as a heuristic for optimizing the mixed-integer program presented in (<ref>).
This is because it produces a score that can be normalized and interpreted as the probability of relevance.
The goal of this heuristic is to find a set of objects $O \subset \mathcal{O}(S)$ such that the NLP $P(g, S)$ is feasible with all actions $a_{1:K}$ of the objects that are operated on $o \in O$ (i.e., $\theta(a) \in O$).
In this way, the search space for finding a feasible NLP $P(g, S)$ could be reduced to a certain level, which improves the time complexity of the whole computation.

§.§ Admissible Heuristic

The classifier may serve as an admissible heuristic.
In this context, the term admissible suggests that errors occurring during classification cannot inhibit the discovery of a solution to a viable problem, denoted as $P(g, S)$.
This is realized by first evaluating $t(I, p_{0,...,L}, q_{i}) = p_{i}$ for all $o_{i} \in \mathcal{O}$, where $q_{i}$ is the canonical view of $o_{i}$.
The NLP $P(g, S)$ first searches for a feasible path with actions only considering the objects in $O \subset \mathcal{O}$.
If it is feasible, then a solution will be given, otherwise the whole object set will be used for searching for a feasible path.
In the instance of an infeasible problem, where no viable path leads to the goal state, all action-object tuples will be examined. Consequently, employing the learned network as an admissible heuristic offers no advantage in such scenarios.

§.§ Non-Admissible Heuristic

Similar to real-world settings, many objects are irrelevant for the goal state in our experiments.
In addition, as long as the object set predicted by the classifier does not contain all of the objects that the NLP has to consider in order to find a feasible path, the NLP will clearly be unable to find a feasible path.
In Sec. <ref> we mentioned the hyperparameter $\beta$ as the threshold for the binary output $Z_{S}(i)$ (Eq. <ref>).
In practice we can reduce $\beta$ to a certain amount to increase the size of object set, and to have a higher probability to fully cover the object set which is required by a feasible path finding procedure.
Because the predicted object set is still a subset of the whole object set $\mathcal{O}$, the search space remains limited, thus the time complexity is expected to be improved.
During the experiments, we assess and analyze both the false irrelevant rate and the time savings concerning the relevance threshold.
For $\beta = 0$, this method becomes identical to the admissible heuristic of Section <ref>.

An example of scene for the on top predicate.

2c    3caccuracy of classifier     3cperformance on LGP

data     predicate     $\eta$     2*true rel.     2*true irrel.     2*total accuracy      2*runtime(s)     2*our runtime(s)     2*speed-up

on-the-left       0.86     87%     94%     93%     $3.06 \pm 0.100$     $1.83 \pm 0.065$     40.2%

on-the-right      0.86     88%     94%     93%     $3.37 \pm 0.104$     $1.87 \pm 0.061$     44.8%

in-front          0.86     92%     96%     95%     $3.64 \pm 0.105$     $2.09 \pm 0.062$     42.6%

behind            0.86     92%     95%     95%     $3.36 \pm 0.104$     $1.92 \pm 0.065$     42.8%
total      0.86     90%     95%     94%     $4.73$     $1.93$     42.6%

on     0.66     76%     89%     87%     $4.47 \pm 0.103$     $3.70 \pm 0.101$     17.2%

Accuracies on all validation sets and Running time comparison of LGP using vision transformer based deep learning model and original LGP. We report 67% confidence intervals.

§ EXPERIMENTS

In our experiments, we consider two types of predicates: predicates describing objects on the same horizontal plane, such as on-the-left and on-the-right, and predictates describing objects in different planes, here we choose on top.

§.§ Data Generation

For all experiments, we place objects randomly on the table and a camera looking down 45 degrees to get the scene image.
We assume that each distinct object has a different color.
For the on-top predicate, we also randomly select a target object and then randomly place another object on top of the target object.
This is done in order to increase the likelihood that an object in the scene is task-relevant.
Randomly distributing all objects would in almost all cases result in all objects being on the same horizontal plane, which would mean that for the on-top predicate, no other objects apart from those being part of the goal are relevant.
Placing one object on top of the target ensures that the classifier is forced to make a decision on whether it needs to be removed.
This will depend on whether there is enough space for the goal object or not, and is therefore a nontrivial visual task.
Fig. <ref> depicts a typical scene for the same-plane case, where the task is to place object $c$ relative to object $d$ (e.g., on the left-hand side of object $d$). The classifier gets an extra input object $e$, and outputs the probability that object $e$ has to be interacted with in order to complete the mentioned task.
All tasks in such configurations would take between 2 and 6 steps to complete, with the number of steps required depending on the size of the area occupied by the target location, the size of the objects to be placed, and the spatial arrangement.

Fig <ref> shows an example of the different-plane case.

For generating labels, we rely on the outcome of LGP.
Given this initial scene and a goal predicate indicating a spatial relationship between two objects, the LGP agent has to figure out both the sequence of high-level decisions and the motion plans that can realize those tasks jointly.
If the agent successfully finds a solution, the objects which have been manipulated according to the motion plan will be labeled as relevant and recorded.
If no plan has been found, nothing is recorded.

As we mentioned before, for any given goal predicate, the majority of objects in the same scene $S$ are irrelevant.
Training a single model to output all $4$ predicates simultaneously would mean that there are $2^4=16$ different labels, some of which (for example the all true case) never occur, resulting in a very imbalanced dataset.
We therefore decide to train 4 different model heads, the PEs, on a single binary label that is independently balanced.
We still share weights of the model body, the details of which are described in section  <ref>.
This also means that the model scales gracefully to a large number of predicates.

[node distance = 0cm, auto]

(transformer) [transformer, above=1cm of pps] Transformer;

(readout) [readout, above=1.2cm of transformer] PEs;
(mlp1) [mlp, above=2.5cm of transformer, xshift=-5.9cm, fill=red!25] PE1;
(mlp2) [mlp, right=0.5cm of mlp1, fill=teal!25] PE2;
(mlp3) [mlp, right=0.5cm of mlp2, fill=blue!25] PE3;
(mlp4) [mlp, right=0.5cm of mlp3, fill=brown!25] PE4;

(dataset, minimum height=2.6cm) [dataset, above=1.2cm of readout] Datasets;
(data1) [mlp, above=2.5cm of readout, xshift=-5.9cm, fill=red!50] dataset1;
(data2) [mlp, right=0.5cm of data1, fill=teal!50] dataset2;
(data3) [mlp, right=0.5cm of data2, fill=blue!50] dataset3;
(data4) [mlp, right=0.5cm of data3, fill=brown!50] dataset4;

[arrow, line width=0.1mm, color=red] ([yshift=0cm,xshift=-1cm]mlp1.south) – ++(0,-1.6) -| ([yshift=0cm,xshift=-1cm]transformer.north);
[arrow, line width=0.1mm, color=teal] ([yshift=0cm,xshift=-1cm]mlp2.south) – ++(0,-0.7) -| ([yshift=0cm,xshift=-1cm]transformer.north);
[arrow, line width=0.1mm, color=blue] ([yshift=0cm,xshift=-1cm]mlp3.south) – ++(0,-0.9) -| ([yshift=0cm,xshift=-1cm]transformer.north);
[arrow, line width=0.1mm, color=brown] ([yshift=0cm,xshift=-1cm]mlp4.south) – ++(0,-1.1) -| ([yshift=0cm,xshift=-1cm]transformer.north);

[arrow, line width=0.1mm, color=red] ([yshift=0cm,xshift=0cm]mlp1.south) – ++(0,-1.6) -| ([yshift=0cm,xshift=0cm]transformer.north);
[arrow, line width=0.1mm, color=teal] ([yshift=0cm,xshift=0cm]mlp2.south) – ++(0,-0.7) -| ([yshift=0cm,xshift=0cm]transformer.north);
[arrow, line width=0.1mm, color=blue] ([yshift=0cm,xshift=0cm]mlp3.south) – ++(0,-0.9) -| ([yshift=0cm,xshift=0cm]transformer.north);
[arrow, line width=0.1mm, color=brown] ([yshift=0cm,xshift=0cm]mlp4.south) – ++(0,-1.1) -| ([yshift=0cm,xshift=0cm]transformer.north);

[arrow, line width=0.1mm, color=red] ([yshift=0cm,xshift=1cm]mlp1.south) – ++(0,-1.6) -| ([yshift=0cm,xshift=1cm]transformer.north);
[arrow, line width=0.1mm, color=teal] ([yshift=0cm,xshift=1cm]mlp2.south) – ++(0,-0.7) -| ([yshift=0cm,xshift=1cm]transformer.north);
[arrow, line width=0.1mm, color=blue] ([yshift=0cm,xshift=1cm]mlp3.south) – ++(0,-0.9) -| ([yshift=0cm,xshift=1cm]transformer.north);
[arrow, line width=0.1mm, color=brown] ([yshift=0cm,xshift=1cm]mlp4.south) – ++(0,-1.1) -| ([yshift=0cm,xshift=1cm]transformer.north);

[arrow, line width=0.1mm, color=red] ([yshift=0cm,xshift=-1cm]data1.south) – ++(0,-1.6) -| node [arrow_txt] 1. ([yshift=0cm,xshift=-1cm]mlp1.north);
[arrow, line width=0.1mm, color=teal] ([yshift=0cm,xshift=-1cm]data2.south) – ++(0,-1.6) -| node [arrow_txt] 2. ([yshift=0cm,xshift=-1cm]mlp2.north);
[arrow, line width=0.1mm, color=blue] ([yshift=0cm,xshift=-1cm]data3.south) – ++(0,-1.6) -| node [arrow_txt] 3. ([yshift=0cm,xshift=-1cm]mlp3.north);
[arrow, line width=0.1mm, color=brown] ([yshift=0cm,xshift=-1cm]data4.south) – ++(0,-1.6) -| node [arrow_txt] 4. ([yshift=0cm,xshift=-1cm]mlp4.north);

[arrow, line width=0.1mm, color=red] ([yshift=0cm,xshift=0cm]data1.south) – ++(0,-1.6) -| node [arrow_txt] 5. ([yshift=0cm,xshift=0cm]mlp1.north);
[arrow, line width=0.1mm, color=teal] ([yshift=0cm,xshift=0cm]data2.south) – ++(0,-1.6) -| node [arrow_txt] 6. ([yshift=0cm,xshift=0cm]mlp2.north);
[arrow, line width=0.1mm, color=blue] ([yshift=0cm,xshift=0cm]data3.south) – ++(0,-1.6) -| node [arrow_txt] 7. ([yshift=0cm,xshift=0cm]mlp3.north);
[arrow, line width=0.1mm, color=brown] ([yshift=0cm,xshift=0cm]data4.south) – ++(0,-1.6) -| node [arrow_txt] 8. ([yshift=0cm,xshift=0cm]mlp4.north);

[arrow, line width=0.1mm, color=red] ([yshift=0cm,xshift=1cm]data1.south) – ++(0,-1.6) -| node [arrow_txt] 9... ([yshift=0cm,xshift=1cm]mlp1.north);
[arrow, line width=0.1mm, color=teal] ([yshift=0cm,xshift=1cm]data2.south) – ++(0,-1.6) -| node [arrow_txt] 10... ([yshift=0cm,xshift=1cm]mlp2.north);
[arrow, line width=0.1mm, color=blue] ([yshift=0cm,xshift=1cm]data3.south) – ++(0,-1.6) -| node [arrow_txt] 11... ([yshift=0cm,xshift=1cm]mlp3.north);
[arrow, line width=0.1mm, color=brown] ([yshift=0cm,xshift=1cm]data4.south) – ++(0,-1.6) -| node [arrow_txt] 12... ([yshift=0cm,xshift=1cm]mlp4.north);

The batch polling training Procedure

§.§ Training the multi-head model

We use $4$ different model heads, one for each predicate.
The body of the model is shared between these.
The model is trained end-to-end by periodically taking a batch of data from each predicate's corresponding training set (Fig. <ref>).

In an alternative approach, we also tried training each predicate head to convergence, before moving on to the next one.
Unsurprisingly, the first approach we described works significantly better and used throughout the paper.

predicate     relevant (pct.)     irrelevant (pct.)
4*Training set
on-the-left         1235677 (13.7%)     7764323 (86.3%)

on-the-right        1239533 (13.8%)     7760467 (86.2%)

in-front            1239388 (13.8%)     7760612 (86.2%)

behind              1248629 (13.9%)     7751371 (86.1%)

4*Test set
on-the-left         12448 (13.8%)     77552 (86.2%)

on-the-right        12557 (14.0%)     77443 (86.0%)

in-front            12401 (13.8%)     77599 (86.2%)

behind              12495 (13.9%)     77505 (86.1%)

Training set     on       3650855 (30.4%)     8349145 (69.6%)

Test set         on       34689 (28.9%)     85311 (71.1%)

Distribution of labels in each dataset

§.§ Predicates in the same plane

For the first phase of the experiments, we consider datasets for on-the-left, on-the-right, in-front and behind.
Each dataset consists of 80,000 scenes in which there are between 3 and 10 objects in the scenes. Both of the training set and test set are generated with the same strategy.
As mentioned before, the data imbalance is quite severe, with only about 14% of the queries being relevant and about 86% of the queries being irrelevant.
Their data distribution is shown in the Table <ref>.
We use $\eta = 0.86$ (Eq. <ref>) to counteract the imbalance in the training dataset.
The network is trained with the Adam optimizer, using a learning rate of 0.00001 and a batch size of 50.

§.§ On-top predicate

The training data contains 160,000 scenes, and the scenes also have various number of objects.
30% of the objects were relevant, 70% were irrelevant.
The same imbalanced data issue arises here as well, and we use $\eta = 0.66$.
Again, the network is trained with the Adam optimizer, using a learning rate of 0.00001 and a batch size of 50.

[trim axis left, trim axis right]

scale only axis,
xmin=-0.05, xmax=0.7,
legend style=font=, nodes=scale=0.8, transform shape, fill=none,
y axis style/.style=
yticklabel style=#1,
ylabel style=#1,
y axis line style=#1,
ytick style=#1
every x tick label/.append style=font=,,
every y tick label/.append style=font=,,

tick label style=font=, /pgf/number format/fixed,
axis y line*=left,
x label style=at=(axis description cs:0.5,0.1),anchor=north,
y label style=at=(axis description cs:0.2,.5),anchor=south,
xlabel=Threshold $\beta$,

width= 0.8,
ymin=0, ymax=1.05,
xtick=0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
ytick=0, 0.2, 0.4, 0.6, 0.8, 1,
(0,   0.8)
(0.065, 0.81)
(0.13,  0.82)
(0.195, 0.83)
(0.26,  0.84)
(0.325, 0.85)
(0.39,  0.86)
(0.455, 0.87)
(0.52,  0.88)
(0.585, 0.89)
(0.65,  0.90)
(0,   0.91)
(0.065, 0.90)
(0.13,  0.89)
(0.195, 0.88)
(0.26,  0.86)
(0.325, 0.85)
(0.39,  0.84)
(0.455, 0.82)
(0.52,  0.81)
(0.585, 0.79)
(0.65,  0.77)
[smooth, color=black,]
(0,   0.88)
(0.065, 0.87)
(0.13,  0.87)
(0.195, 0.86)
(0.26,  0.86)
(0.325, 0.85)
(0.39,  0.84)
(0.455, 0.84)
(0.52,  0.83)
(0.585, 0.82)
(0.65,  0.81)

tick label style=font=,
axis y line*=right,
axis x line=none,
width= 0.8,
ymin=0, ymax=1.05,
ytick=0, 0.2, 0.4, 0.6, 0.8, 1,

legend style=at=(0.25,0.69), anchor=north west,
every y tick label/.append style=font=, xshift=-0.5ex,
scaled ticks=false,
tick label style=/pgf/number format/fixed

/pgfplots/refstyle=trplottrue relevant
/pgfplots/refstyle=tiplottrue irrelevant
/pgfplots/refstyle=total_acctotal accuracy
(0, 0.2)
(0.065, 0.19)
(0.13, 0.17)
(0.195, 0.17)
(0.26, 0.16)
(0.325, 0.15)
(0.39, 0.14)
(0.455, 0.13)
(0.52, 0.12)
(0.585, 0.11)
(0.65, 0.10)
; false irrelevant

Classifier accuracies

[trim axis left, trim axis right]

scale only axis,
xmin=-0.05, xmax=0.7,
legend style=font=, nodes=scale=0.8, transform shape,
y axis style/.style=
yticklabel style=#1,
ylabel style=#1,
y axis line style=#1,
ytick style=#1
every x tick label/.append style=font=,,
every y tick label/.append style=font=,,

tick label style=font=, /pgf/number format/fixed,
axis y line*=left,
ylabel near ticks, yticklabel pos=right,
x label style=at=(axis description cs:0.5,0.1),anchor=north,
y label style=at=(axis description cs:1.3,.5),anchor=south,

xlabel=Threshold $\beta$,

width= 0.8,
ymin=-5, ymax=40,
xtick=0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7,
ytick=0, 10, 20, 30, 40,

tick label style=font=, /pgf/number format/fixed,
axis y line*=right,
axis x line=none,
width= 0.8,
ymin=-5, ymax=40,

ytick=0, 10, 20, 30, 40,

every y tick label/.append style=font=, xshift=-0.5ex,

[smooth, color=Red1,]
(0.0, 17.2)
(0.065, 19.0)
(0.13,  31.0)
(0.195, 16.0)
(0.26,  23.3)
(0.325, 22.1)
(0.39,  23.1)
(0.455, 22.3)
(0.52,  22.1)
(0.585, 18.2)
(0.65,  20.4)

LGP speed-up

§.§ Efficiency Improvement on LGP

Finally, we compare the improvement brought by the model for sequential planning of actions using LGP.
We found out that the accuracies for different models are comparable (Tbl. <ref>).
Since the runtime of LGP is only dependent on the accuracy of the model, we will focus on the results obtained using the model based on the Transformer.
The results for both same-plane and different-plane datasets are shown in Tbl. <ref>.
We use the unmodified LGP as the baseline in our experiments.
We do not show the impact of false irrelevant predictions of the classifier on the performance of the LGP.
The influence of false irrelevant predictions on the runtime depends on other design choices like maximum runtime settings, which we factor out here.

The model performs better on the first dataset than on the second dataset, where we find an over 40% improvement on LGP.
The second dataset has only 76 percent of true relevant due to the complexity of the predicates and the occlusion problem of the visual sensors, but it still provides a 17% improvement.

As mentioned in section  <ref>, $\beta$ directly affects the percentage of false irrelevant predictions, and the presence of false irrelevant prediction means that LGP has a high probability of not finding a feasible manipulation sequence.

We analyse how $\beta$ would affect the runtime of LGP, and the results are shown in Fig. <ref> and Fig. <ref>.
The data in Fig. <ref> is obtained by using the classifier trained on the data for the $on-top$ predicate.
As $\beta$ increases, the model will adopt a more conservative strategy, so that more data will be marked as relevant.
This also leads to a decrease in the number of false irrelevant cases.
Since the action planning of LGP only considers the objects marked as "relevant", missing the objects that should be considered will result in LGP not finding a feasible plan.
By adjusting $\beta$, we can avoid this situation to some extent.

§.§ Comparison of base models

In the last part of the experiments we compare the training results based on different visual deep learning models.
All results in the experiments are obtained for the more complicated $on-top$ predicate dataset.

2c    3caccuracy of classifier

2*Base Model     $\eta$     2*true rel.     2*true irrel.     2*total acc.

Transformer     0.66     81%     89%     87%

3D-ResNet     0.66     79%     90%     87%

Default ResNet     0.66     81%     91%     88%

Comparison of accuracy rates using different base models

The best results are obtained using Default ResNet as the base model.
This is followed by the 3D-ResNet-based model, and finally the Transformer-based model.
It is worth noting that even though we used very different vision models, the final results are not significantly different.

§ CONCLUSION

As other visual models, our model can not circumvent problems due to occlusion and other ambiguities. Regardless of how the vision sensors are placed, it is possible for smaller objects to be obscured by larger objects, as is often the case for the on-top predicate.
Furthermore, our dataset is relatively homogeneous in terms of object color selection, and the images are not affected by ambient lighting.
While this can likely be addressed in a straightforward way by using photorealistic renderings and by randomizing lights and colors, it currently limits the number of scenarios the model can be used in.

Sequential manipulation planning becomes inefficient when applied to scenes where many non-relevant objects are present.
We address this in this study by using a learned visual heuristic to classify objects into non-relevant and relevant objects, given a symbolic goal.
We use a deep learning network that uses the scene image and canonical views describing the goal and a queried object as input, and makes predictions about the involvement of this object in the scene.
These predictions can be directly used in the process of sequential manipulation planning, thus directly improving the efficiency of planning.

We experimented on two different types of predicates: the first for spatial relations on the same plane and the second for spatial relations on different planes.
They both perform well in terms of prediction accuracy and reduce the time required for sequential manipulation planning.

[Wells et al., 2019]
A. M. Wells, N. T. Dantam, A. Shrivastava, and L. E. Kavraki, “Learning
feasibility for task and motion planning in tabletop environments,”
IEEE robotics and automation letters, vol. 4, no. 2, pp. 1255–1262,

[Toussaint et al., 2018]
M. A. Toussaint, K. R. Allen, K. A. Smith, and J. B. Tenenbaum,
“Differentiable physics and stable modes for tool-use and manipulation
planning,” 2018.

[Rodrıguez et al., 2019]
I. Rodrıguez, K. Nottensteiner, D. Leidner, M. Kaßecker, F. Stulp, and
A. Albu-Schäffer, “Iteratively refined feasibility checks in robotic
assembly sequence planning,” IEEE Robotics and Automation Letters,
vol. 4, no. 2, pp. 1416–1423, 2019.

[Driess et al., 2019]
D. Driess, O. Oguz, and M. Toussaint, “Hierarchical task and motion planning
using logic-geometric programming (hlgp),” in RSS Workshop on Robust
Task and Motion Planning, 2019.

[Kase et al., 2020]
K. Kase, C. Paxton, H. Mazhar, T. Ogata, and D. Fox, “Transferable task
execution from pixels through deep planning domain learning,” in 2020
IEEE International Conference on Robotics and Automation (ICRA).1em
plus 0.5em minus 0.4emIEEE, 2020, pp. 10 459–10 465.

[Driess et al., 2021]
D. Driess, J.-S. Ha, R. Tedrake, and M. Toussaint, “Learning geometric
reasoning and control for long-horizon tasks from visual input,” in
2021 IEEE International Conference on Robotics and Automation
(ICRA).1em plus 0.5em minus 0.4emIEEE, 2021, pp.
14 298–14 305.

[Driess et al., 2020]
D. Driess, O. Oguz, J.-S. Ha, and M. Toussaint, “Deep visual heuristics:
Learning feasibility of mixed-integer programs for manipulation planning,”
in 2020 IEEE International Conference on Robotics and Automation
(ICRA).1em plus 0.5em minus 0.4emIEEE, 2020, pp. 9563–9569.

[Yuan et al., 2022]
W. Yuan, C. Paxton, K. Desingh, and D. Fox, “Sornet: Spatial object-centric
representations for sequential manipulation,” in Conference on Robot
Learning.1em plus 0.5em minus 0.4emPMLR, 2022, pp. 148–157.

[Alili et al., 2010]
S. Alili, A. K. Pandey, E. A. Sisbot, and R. Alami, “Interleaving symbolic and
geometric reasoning for a robotic assistant,” in ICAPS Workshop on
Combining Action and Motion Planning, 2010.

[Srivastava et al., 2014]
S. Srivastava, E. Fang, L. Riano, R. Chitnis, S. Russell, and P. Abbeel,
“Combined task and motion planning through an extensible planner-independent
interface layer,” in 2014 IEEE international conference on robotics
and automation (ICRA).1em plus 0.5em minus 0.4emIEEE, 2014,
pp. 639–646.

[de Silva et al., 2013]
L. de Silva, A. K. Pandey, M. Gharbi, and R. Alami, “Towards combining htn
planning and geometric task planning,” arXiv preprint
arXiv:1307.1482, 2013.

[Dantam et al., 2018]
N. T. Dantam, Z. K. Kingston, S. Chaudhuri, and L. E. Kavraki, “An incremental
constraint-based framework for task and motion planning,” The
International Journal of Robotics Research, vol. 37, no. 10, pp. 1134–1151,

[Siméon et al., 2004]
T. Siméon, J.-P. Laumond, J. Cortés, and A. Sahbani, “Manipulation
planning with probabilistic roadmaps,” The International Journal of
Robotics Research, vol. 23, no. 7-8, pp. 729–746, 2004.

[Erdem et al., 2011]
E. Erdem, K. Haspalamutgil, C. Palaz, V. Patoglu, and T. Uras, “Combining
high-level causal reasoning with low-level geometric reasoning and motion
planning for robotic manipulation,” in 2011 IEEE International
Conference on Robotics and Automation.1em plus 0.5em minus
0.4emIEEE, 2011, pp. 4575–4581.

[Lagriffoul et al., 2012]
F. Lagriffoul, D. Dimitrov, A. Saffiotti, and L. Karlsson, “Constraint
propagation on interval bounds for dealing with geometric backtracking,” in
2012 IEEE/RSJ International Conference on Intelligent Robots and
Systems.1em plus 0.5em minus 0.4emIEEE, 2012, pp. 957–964.

[Lagriffoul et al., 2014]
F. Lagriffoul, D. Dimitrov, J. Bidot, A. Saffiotti, and L. Karlsson,
“Efficiently combining task and motion planning using geometric
constraints,” The International Journal of Robotics Research,
vol. 33, no. 14, pp. 1726–1747, 2014.

[Lozano-Pérez and Kaelbling, 2014]
T. Lozano-Pérez and L. P. Kaelbling, “A constraint-based method for
solving sequential manipulation planning problems,” in 2014 IEEE/RSJ
International Conference on Intelligent Robots and Systems.1em plus
0.5em minus 0.4emIEEE, 2014, pp. 3684–3691.

[Shoukry et al., 2016]
Y. Shoukry, P. Nuzzo, I. Saha, A. L. Sangiovanni-Vincentelli, S. A. Seshia,
G. J. Pappas, and P. Tabuada, “Scalable lazy smt-based motion planning,” in
2016 IEEE 55th Conference on Decision and Control (CDC).1em
plus 0.5em minus 0.4emIEEE, 2016, pp. 6683–6688.

[Hadfield-Menell et al., 2016]
D. Hadfield-Menell, C. Lin, R. Chitnis, S. Russell, and P. Abbeel, “Sequential
quadratic programming for task plan optimization,” in 2016 IEEE/RSJ
International Conference on Intelligent Robots and Systems (IROS).1em
plus 0.5em minus 0.4emIEEE, 2016, pp. 5040–5047.

[Toussaint, 2015]
M. Toussaint, “Logic-geometric programming: An optimization-based approach to
combined task and motion planning,” in Twenty-Fourth International
Joint Conference on Artificial Intelligence, 2015.

[Toussaint et al., 2020]
M. Toussaint, J.-S. Ha, and D. Driess, “Describing physics for physical
reasoning: Force-based sequential manipulation planning,” IEEE
Robotics and Automation Letters, vol. 5, no. 4, pp. 6209–6216, 2020.

[Ha et al., 2020]
J.-S. Ha, D. Driess, and M. Toussaint, “A probabilistic framework for
constrained manipulations and task and motion planning under uncertainty,”
in 2020 IEEE International Conference on Robotics and Automation
(ICRA).1em plus 0.5em minus 0.4emIEEE, 2020, pp. 6745–6751.

[Hogan et al., 2018]
F. R. Hogan, E. R. Grau, and A. Rodriguez, “Reactive planar manipulation with
convex hybrid mpc,” in 2018 IEEE International Conference on Robotics
and Automation (ICRA).1em plus 0.5em minus 0.4emIEEE, 2018,
pp. 247–253.

[Deits and Tedrake, 2014]
R. Deits and R. Tedrake, “Footstep planning on uneven terrain with
mixed-integer convex optimization,” in 2014 IEEE-RAS international
conference on humanoid robots.1em plus 0.5em minus 0.4emIEEE,
2014, pp. 279–286.

[Vaswani et al., 2017]
A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez,
Ł. Kaiser, and I. Polosukhin, “Attention is all you need,”
Advances in neural information processing systems, vol. 30, 2017.

[Dosovitskiy et al., 2020]
A. Dosovitskiy, L. Beyer, A. Kolesnikov, D. Weissenborn, X. Zhai,
T. Unterthiner, M. Dehghani, M. Minderer, G. Heigold, S. Gelly et al.,
“An image is worth 16x16 words: Transformers for image recognition at
scale,” arXiv preprint arXiv:2010.11929, 2020.

[Carion et al., 2020]
N. Carion, F. Massa, G. Synnaeve, N. Usunier, A. Kirillov, and S. Zagoruyko,
“End-to-end object detection with transformers,” in European
conference on computer vision.1em plus 0.5em minus 0.4emSpringer, 2020, pp. 213–229.

[Guo et al., 2022]
Y. Guo, S. Yang, W. Chen, L. Ma, D. Xie, and S. Pu, “2nd place solution for
iccv 2021 vipriors image classification challenge: An attract-and-repulse
learning approach,” ArXiv, vol. abs/2206.06168, 2022.

[Nguyen et al., 2020]
S. Nguyen, O. S. Oguz, V. N. Hartmann, and M. Toussaint, “Self-supervised
learning of scene-graph representations for robotic sequential manipulation
planning.” in CoRL, 2020, pp. 2104–2119.

[He et al., 2016]
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.

[Toussaint and Lopes, 2017]
M. Toussaint and M. Lopes, “Multi-bound tree search for logic-geometric
programming in cooperative manipulation domains,” in 2017 IEEE
International Conference on Robotics and Automation (ICRA).1em plus
0.5em minus 0.4emIEEE, 2017, pp. 4044–4051.

[Mason, 1985]
M. Mason, “The mechanics of manipulation,” in Proceedings. 1985 IEEE
International Conference on Robotics and Automation, vol. 2.1em plus
0.5em minus 0.4emIEEE, 1985, pp. 544–548.


# Augmenting MRI scan data with real-time predictions of glioblastoma brain
tumor evolution using exponential time integrators

Magdalena Pabisz1, Judit Muñoz-Matute2,3, Maciej Paszyński1 1 AGH University
of Krakow, Poland
2 The Basque Center for Applied Mathematics (BCAM), Bilbao, Spain
3 The Oden Institute for Computational Engineering and Sciences,
The University of Texas at Austin, USA

###### Abstract

We present an ultra-fast simulator to augment the MRI scan imaging of
glioblastoma brain tumors with predictions of future evolution. We consider
the glioblastoma tumor growth model based on the Fisher-Kolmogorov diffusion-
reaction equation with logistic growth. For the discretization we employ
finite differences in space coupled with a time integrator in time employing
the routines from [2] to compute the actions of the exponentials of the linear
operator. By combining these methods, we can perform the prediction of the
tumor evolution for several months forward within a couple of seconds on a
modern laptop. This method does not require HPC supercomputing centers, and it
can be performed on the fly using a laptop with Windows 10, Octave
simulations, and ParaView visualization. We illustrate our simulations by
predicting the tumor growth evolution based on three-dimensional MRI scan
data.

###### keywords:

Glioblastoma brain tumor , MRI scan data , Prediction of tumor evolution ,
Exponential time integrators

##  1 Introduction

Glioblastoma is a malignant brain tumor with a high mortality rate [7]. This
tumor is highly aggressive, and it generates microvascular proliferation that
is not visible on MRI scans [20]. Thus, computer-based simulation predictions
of the evolution of the brain tumor are important in planning treatment and
surgery. But highly efficient simulators of the tumor growth take several
minutes to compute a single time step [12, 24, 23, 11].The Physics Informed
Neural Networks originally proposed by Karniadakis et al. [18, 4, 3]
instantiated for the brain tumor simulation in [22] can predict the future
behavior of the glioblastoma tumor evolution within one hour for the patient-
specific case.

In this work, we present a computational tool to perform ultra-fast
simulations of three dimensional progression of glioblastoma brain tumor
models. We consider the Fisher-Kolmogorov diffusion-reaction equation with
logistic growth [21] describing the evolution of glioblastoma tumor in a non-
homogenous environment of the human head. Our simulations employs the MRI scan
data consisting of 29 bitmaps, with pixel’s intensity varying from 0 to 255.
We introduce materials such as skull, skin, white matter, gray matter, and
air, based on the pixels intensity.

For the discretization of the transient semilinear model we combine two
methods. First, we semidiscretize the space variable employing finite
diffences. Then, we discretize the resulting semilinear system of ordinary
differential equations by an exponential integrator method [8]. For
simplicity, we focus on the exponential Euler method which is first order but
the extension to higher order exponential-based methods in time is
straightforward. We employ the routine presented in [2] for computing the
action of the corresponding $\varphi$-functions over vector. Then, we post-
process the results and we use ParaView for visualization. We test this
strategy on both 2D and 3D non-stationary glioblastoma models employing MRI
data for the diffusion coefficient. The numerical results show that we can
perform $100$ iterations over $32\times 32\times 32$ computational mesh in
less than 3 seconds on a laptop. Thus, this simulator can be employed "on the
fly".

The efficiency comes from the fact that the exponential routines are extremely
efficient on operators coming from a semidiscretization in space employing
finite differences. Moreover, the exponential time integrators [9, 16, 5] are
designed to solve semilinear problems that the one presented in this article.
They are suitable for long-time simulations and they are unconditionally
stable. From the authors understanding, there is only one work in the
literature [13] where the authors employ exponential time integrators for the
simulation of cancer models. This article aims to underscore the potential of
this tool in the simulation of glioblastomas in the human brain.

Our future wrok will include development of the parallel version of the
simulator, targeting GPGPU cards. This will allow us to perform simulations on
high resolution data of the brain. We also plan to incorporate the drug
delivery into our model [15]. Having the chemotherapy included into our
simulator, we plan to attack the data assimilation problem using the
supermodeling technique [17, 19].

##  2 Glioblastoma tumor model

Let $\Omega\subset{\mathbb{R}}^{3}$ an open set and $I=(0,T)$ with $T>0$, we
consider the following Fisher-Kolmogorov diffusion-reaction equation with
logistic growth [21] describing the brain tumor dynamics.

$\displaystyle{\left\\{\begin{aligned} \frac{\partial u}{\partial
t}=\underbrace{\nabla\cdot(D(\mathbf{x})\nabla u)}_{\textrm{Tumor cell
diffusion}}&+\underbrace{\rho u(1-u)}_{\textrm{Tumor cell
proliferation}},&\mbox{in}&\;\Omega\times I,\\\ \nabla u\cdot
n&=0,&\mbox{on}&\;\partial\Omega\times I,\\\
u(\mathbf{x},0)&=u_{0},&\mbox{on}&\;\Omega\times\\{0\\}.\\\
\end{aligned}\right.}$ (1)

Here, $u(x,y,z;t)$ represents the tumor cell density scalar field, and
$D(x,y,z)$ represents the diffusion coefficient, describing how the tumor
cells expand in a tissue, similarly to the diffusion phenomena. The $\rho$
coefficient describes the proliferation rate of the tumor cells (how fast they
multiply). Modeling the tumor cell dynamics by this kind of PDE results in a
progression of the tumor "traveling wave" with the velocity $2\sqrt{D}\rho$,
according to [22].

The diffusion coefficient usually patient-specific. Following [6], it
expressed as a linear combination of gray and white matter coefficients,
namely $D=p_{w}D_{w}+p_{g}D_{g}$, where $p_{w},p_{g}$ corresponds to the
proportions of the white and gray matter, obtained from the MRI scans. Here,
we considered the case of isotropic diffusion in both the white and the grey
matter [10]. Following [14] and [22], we assume that $D_{w}=10D_{g}$, and we
select $D_{w}=0.13$ [mm2days-1], $D_{g}=0.013$ [mm2days-1]. Moreover,
following [6] we set $\rho=0.025$ [days]-1. With this setup, assuming 70
percent of the gray matter and 30 percent of the white matter, we expect the
tumor growth to progress with a velocity of
$2\sqrt{D}\rho=2\sqrt{p_{w}D_{w}+p_{g}D_{g}}\rho=0.11$[mm days-1]. Thus, we
expect the tumor to progress to 11 mm per 100 days.

##  3 Finite differences

We first semi-discretize (1) in space employing finite differences. Namely, we
introduce a regular grid with equidistant points in each espatial direction:

$\\{x_{i,j,k}=((i-1)h,(j-1)h,(k-1)h)\\}_{i=1,...,N_{x};j=1,...,N_{y};k=1,...,N_{z}},$
(2)

and we represent the values of the tumor cell densities at these point at time
moment $t$ as

$\\{u^{t}_{i,j,k}=u(x_{i,j,k};t)\\}_{i=1,...,N_{x};j=1,...,N_{y};k=1,...,N_{z}}.$
(3)

We discretize the equation at time moment $t$ using finite differences as
follows

$\displaystyle\begin{aligned} \frac{\partial u^{t}_{i,j,k}}{\partial
t}=&\frac{\partial D(x_{i,j,k})}{\partial x_{1}}\frac{\partial
u^{t}_{i,j,k}}{\partial
x_{1}}+D(x_{i,j,k})\frac{\partial^{2}u^{t}_{i,j,k}}{\partial x_{1}^{2}}+\\\
&\frac{\partial D(x_{i,j,k})}{\partial x_{2}}\frac{\partial
u^{t}_{i,j,k}}{\partial
x_{2}}+D(x_{i,j,k})\frac{\partial^{2}u^{t}_{i,j,k}}{\partial x_{2}^{2}}\\\
&\frac{\partial D(x_{i,j,k})}{\partial x_{3}}\frac{\partial
u^{t}_{i,j,k}}{\partial
x_{3}}+D(x_{i,j,k})\frac{\partial^{2}u^{t}_{i,j,k}}{\partial x_{3}^{2}}+\rho
u^{t}_{i,j,k}(1-u^{t}_{i,j,k}),\end{aligned}$ (4)

with

$\displaystyle\begin{aligned} \frac{\partial u^{t}_{i,j,k}}{\partial
x_{1}}&=\frac{u^{t}_{i+1,j,k}-u^{t}_{i,j,k}}{h},\\\ \frac{\partial
u^{t}_{i,j,k}}{\partial x_{2}}&=\frac{u^{t}_{i,j+1,k}-u^{t}_{i,j,k}}{h},\\\
\frac{\partial u^{t}_{i,j,k}}{\partial
x_{3}}&=\frac{u^{t}_{i,j,k+1}-u^{t}_{i,j,k}}{h},\\\
\frac{\partial^{2}u^{t}_{i,j,k}}{\partial
x_{1}^{2}}=&\frac{u^{t}_{i+1,j,k}-2u^{t}_{i,j,k}+u^{t}_{i-1,j,k}}{h^{2}},\\\
\frac{\partial^{2}u^{t}_{i,j,k}}{\partial
x_{2}^{2}}=&\frac{u^{t}_{i,j+1,k}-2u^{t}_{i,j,k}+u^{t}_{i,j-1,k}}{h^{2}},\\\
\frac{\partial^{2}u^{t}_{i,j,k}}{\partial
x_{3}^{2}}=&\frac{u^{t}_{i,j,k+1}-2u^{t}_{i,j,k}+u^{t}_{i,j,k-1}}{h^{2}},\\\
\end{aligned}$ (5)

and we obtain the following system of semilinear Ordinary Differential
Equations

$\displaystyle{\left\\{\begin{aligned}
\dot{U}(t)&=AU(t)+F(U(t)),&\mbox{in}&\;I,\\\ U(0)&=U_{0},\\\
\end{aligned}\right.}$ (6)

where $U(t)=\\{u^{t}_{i,j,k}\\}_{i=1,...,N_{x};j=1,...,N_{y};k=1,...,N_{z}}$
is the time-dependent vector of the degrees of freedom in space. To derive the
$A$ operator, we simplify the derivation, assuming $\frac{\partial
D_{i,j,k}}{\partial x_{i}}=0$,

$\displaystyle\begin{aligned} \frac{\partial u^{t}_{i,j,k}}{\partial
t}&=D(x_{i,j,k})\frac{u^{t}_{i+1,j,k}-2u^{t}_{i,j,k}+u^{t}_{i-1,j,k}}{h^{2}}\\\
&+D(x_{i,j,k})\frac{u^{t}_{i,j+1,k}-2u^{t}_{i,j,k}+u^{t}_{i,j-1,k}}{h^{2}}\\\
&+D(x_{i,j,k})\frac{u^{t}_{i,j,k+1}-2u^{t}_{i,j,k}+u^{t}_{i,j,k-1}}{h^{2}}+\rho
u^{t}_{i,j,k}(1-u^{t}_{i,j,k}),\end{aligned}$ (7)

and we group the terms

$\displaystyle\begin{aligned} h^{2}\frac{\partial u^{t}_{i,j,k}}{\partial
t}&=-6D(x_{i,j,k})u^{t}_{i,j,k}+D(x_{i,j,k})u^{t}_{i-1,j,k}\\\
&+D(x_{i,j,k})u^{t}_{i,j-1,k}+D(x_{i,j,k})u^{t}_{i,j,k-1}\\\
&+D(x_{i,j,k})u^{t}_{i+1,j,k}+D(x_{i,j,k})u^{t}_{i,j+1,k}\\\
&+D(x_{i,j,k})u^{t}_{i,j,k+1}+h^{2}\rho
u^{t}_{i,j,k}(1-u^{t}_{i,j,k}),\end{aligned}$ (8)

The entries of matrix $A$ are

$\displaystyle A_{i,j,k;l,m,n}=\displaystyle{\left\\{\begin{aligned}
-6D(x_{i,j,k}),&\quad(i,j,k)==(l,m,n),\\\
D(x_{i-1,j,k}),&\quad(l,m,n)\in\\{(i-1,j,k),(i+1,j,k),\\\
&\quad\qquad\qquad\quad(i,j-1,k),(i-1,j+1,k),\\\
&\quad\qquad\qquad\quad(i-1,j,k-1),(i-1,j,k+1)\\},\\\
0,&\quad\textrm{otherwise}.\end{aligned}\right.}$ (9)

We employ
$\\{1,...,N_{x}\\}\times\\{1,...,N_{y}\\}\times\\{1,...,N_{z}\\}\ni(i,j,k)\rightarrow
global(i,j,k)=i+(j-1)N_{y}+(j-1)(k-1)N_{y}N_{z}$ is the mapping from the
integer coordinates into the global rows / columns numbering. The $F$ operator
is given by $F(U(t))=\rho U(t)(1-U(t))$.

##  4 Exponential integrators

For the time discretization, we consider a uniform partition of the time
interval as

$0=t_{0}<t_{1}<\ldots<t_{N-1}<t_{N}=T,$ (10)

we define $I_{n}=(t_{n},t_{n+1})$ and $\tau=t_{n+1}-t_{n},\;\forall
n=0,\ldots,N-1$. Let $U_{n}$ be the numerical approximation to the solution of
(6) at $t_{n}$, we know that the integral representation of the solution of
(6) at $t_{n+1}$, also known as the variation-of-constants formula, reads

$U_{n+1}=e^{\tau A}U_{n}+h\int_{0}^{1}e^{(1-\theta)\tau
A}F(U(t_{n}+\tau\theta))d\theta.$ (11)

Different approximations of the nonlinear term in (11) lead to different
exponential time integration methods [9]. All these methods are expressed in
terms of the so-called $\varphi$-functions defined as

$\displaystyle{\left\\{\begin{aligned} \varphi_{0}(z)&=e^{z},\\\
\varphi_{p}(z)&=\int_{0}^{1}e^{(1-\theta)z}\frac{\theta^{p-1}}{(p-1)!}d\theta,\;\forall
p\geq 1,\end{aligned}\right.}$ (12)

which satisfy the following recurrence relation

$\varphi_{p+1}(z)=\frac{1}{z}\left(\varphi_{p}(z)-\frac{1}{p!}\right).$ (13)

Here, we will focus on the simplest exponential integrator method, the
Exponential Euler method, which is first order in time. For higher-order
methods we refer to [9, 16]. For that, we approximate in (11) the non-linear
term with its value at $t_{n}$ that it is known, i.e.,
$F(U(t_{n}+\tau\theta))\approx F(U_{n})$. Integrating exactly in (11), we
obtain

$U_{n+1}=\varphi_{0}(\tau A)U_{n}+\tau\varphi_{1}(\tau A)F(U_{n}).$ (14)

which is given in terms of the $\varphi$-functions (12). Finally, employing
the recurrence formula (13), we rewrite (14) as

$U_{n+1}=U_{n}+\tau\varphi_{1}(\tau A)(F(U_{n})+\tau AU_{n}).$ (15)

For the numerical results, we employ the MATLAB routines from [2] for
computing the action of $\varphi$-functions over vectors. In this routine the
authors employ the scaling and squaring method together with a truncated
Taylor series approximation to the exponential of a matrix. As we show in the
numerical results, this routine applied to operator $\tau A$ coming from
finite difference semidiscretization in space is extremely efficient.
Moreover, exponential integrators are suitable for long-time simulations as
they are unconditionally stable.

##  5 Numerical results

###  5.1 MRI scan of the human brain

Figure 1: MRI scans of the head of Maciej Paszyński.

Our simulations are based on digital data, the MRI scan with 29 two-
dimensional slices, each one with 532 times 565 pixels. Each pixel’s intensity
is a value from the range of [0, 255], and it’s proportional to the material’s
(skull, skin, white matter or gray matter, and air) normalized density.
Exemplary slices of the human head from the MRI scan of the head of one of the
authors (Maciej Paszyński) are presented in Figure 1. Next, according to the
MRI scan data, we employ material data changing on the skull, skin, white
matter, gray matter, and air. We assume air (MRI scan data $\leq$ 1), white
matter (1 $\leq$ approximation $\leq$ 230), grey matter (230 $<$ approximation
$\leq$ 240) and skull (approximation $>$ 240).

###  5.2 Two-dimensional glioblastoma brain tumor simulation

We consider one central slice of the MRI scan data. Let
$\Omega\subset[0,200]^{2}$ [mm2] be an open set and $I=(150,3500)$ [days] the
simulated time interval. We solve the two-dimensional version of the Fisher-
Kolmogorov diffusion-reaction equation with logistic growth describing the
brain tumor dynamics (1) [21]. We prescribe the following initial state

$u(x_{1},x_{2};t_{0})=0.1\exp(-10((x_{1}-110)^{2}+(x_{2}-140)^{2})).$ (16)

We select $\rho=2.5\times 0.01$ [days-1] [6]. The diffusion coefficient for
the white matter is defined as $D_{w}=1.3\times 0.1$ [mm2days-1] [6]. The
diffusion coefficient for the gray matter is defined as $D_{w}=0.13\times 0.1$
[mm2days-1] [10.1007/978-3-642-22092-0_60].

With this parameters setup, we start from $t_{0}=150$ [days] with $T=3500$,
and we perform 100 time steps of the exponential integrator simulation. Thus,
one time step represents $\frac{3500-150}{100}=33.5$ [days]. The spatial mesh
consists of $193\times 193$ points. The 100 time steps exponential integrators
method executed on a laptop with Win10, using Octave, with 11th Gen Intel(R)
Core(TM) i5-11500H @ 2.90GHz, 2.92 GHz, and 32 GB of RAM, takes less than 6
seconds. The selected snapshots from the simulation are presented in Figures
2-3. In these pictures, we do not present the brain tissue, just the density
of the tumor cells.

This single two-dimensional slice of the human head has around 70 percent of
the gray matter, and 30 percent of the white matter. According to our
evaluation, we have the theoretical predicted tumor velocity as of
$2\sqrt{D}\rho=0.11$ [mm days -1]. This means that we have 11mm progression
per 100 days. The progression from day 150 (our initial state, where the tumor
had a diameter of 1mm) to day 518, presented in Figures 2-3 shows the
progression from 1mm to 42mm (as measuring the tumor radius in one dimension).
This implies the tumor velocity of $\frac{42}{518-150}=0.11$ [mm days -1].
This confirms the agreement between the numerical results and the model setup.

Figure 2: Sequence of snapshoots from the two-dimensional tumor brain
simulations.

Figure 3: Sequence of snapshoots from the two-dimensional tumor brain
simulations. Continuation.

###  5.3 Three dimensional glioblastoma brain tumor simulation

We construct 3D computational grid out of 28 slices of the MRI scan. We want
to solve the Fisher-Kolmogorov diffusion-reaction equation with logistic
growth (1) describing the brain tumor dynamics [21] over
$\Omega\subset[0,200]^{3}$ [mm3], and time interval $I=(150,3500)$ [days] and
prescribed initial state

${u(x_{1},x_{2},x_{3};t_{0})=0.1\exp(-10((x_{1}-102)^{2}+(x_{2}-138)^{2}+(x_{3}-96)^{2}))}$
(17)

Following [6, 10.1007/978-3-642-22092-0_60] we consider $\rho=2.5\times 0.01$
[days-1]. $D_{w}=1.3\times 0.1$ [mm2days-1] [6]. $D_{g}=0.13\times 0.1$
[mm2days-1]. We perform 100 time steps, starting from $t_{0}=150$ [days], with
time step $=33.5$ [days]. The spatial mesh consists of $50\times 50\times 50$
points. The 100 time steps exponential integrators method executed on a laptop
with Win10, using Octave, with 11th Gen Intel(R) Core(TM) i5-11500H @ 2.90GHz,
2.92 GHz and 32 GB of RAM, takes less than 3 seconds. The selected snapshots
from the simulation are presented in Figure 4-6. We show the brain tissue and
the growing tumor cells. The visualization has been performed through the
interface with ParaView [1].

The percentage of the gray matter and white matter in the full three-
dimensional MRI scan of the human head is around 85 percent of the gray
matter, and 15 percent of the white matter. This implies the following
theoretical predicted tumor velocity as of $2\sqrt{D}\rho=$
$2\sqrt{p_{w}D_{w}+p_{g}D_{g}}\rho=0.087$ [mm days -1]. This means that we
have 9mm progression per 100 days.

The progression from day 150 (our initial state, where the tumor had a
diameter of 1mm) to day 485, presented in Figure 4 shows the progression from
1mm to 29mm (as measuring the tumor radius in one dimension) . This implies
the tumor velocity of $\frac{28}{485-150}=0.083$ [mm days -1].

(a) Initial condition. Day 150.

(b) Snapshot at day 217.

(c) Snapshot at day 284.

(d) Snapshot at day 351.

(e) Snapshot at day 418.

(f) Snapshot at day 485.

Figure 4: Sequence of snapshots from the 3D tumor brain simulation. Top view.
Snapshots from day 150 (initial state), one time step is 33 days. Each
snapshot is two time steps.

(a) Initial condition. Day 552.

(b) Snapshot at day 619.

(c) Snapshot at day 686.

(d) Snapshot at day 719.

(e) Snapshot at day 786.

(f) Snapshot at day 853.

Figure 5: Sequence of snapshots from the 3D tumor brain simulation. Top view.
Continuation from day Snapshots from day 552, one time step is 33 days. Each
snapshot is two time steps.

(a) Initial condition. Day 920.

(b) Snapshot at day 986.

(c) Snapshot at day 1052.

(d) Snapshot at day 1118.

(e) Snapshot at day 1184.

(f) Snapshot at day 1250.

Figure 6: Sequence of snapshots from the 3D tumor brain simulation. Top view.
Continuation from day Snapshots from day 920, one time step is 33 days. Each
snapshot is two time steps.

##  6 Conclusions

Based on three-dimensional MRI scans of the human head, we build an ultra-fast
simulator of brain tumor growth. We employed the Fisher-Kolmogorov diffusion-
reaction equation with logistic growth to model the glioblastoma tumor. We
used the exponential integrators method, employing the routine from [2] to
compute the action of the exponentials, with a finite difference
discretization in space. Our two-dimensional simulator written in Octave can
predict a one-year forward tumor evolution within 6 seconds using the
resolution of 193$\times$193 pixels of a single two-dimensional slice of the
human head. We also developed a three-dimensional simulator with Octave and
ParaView visualization. The three-dimensional simulator on a laptop can
perform a one-year tumor growth prediction with a resolution of $32\times
32\times 32$ pixels within 3 seconds. For both two-dimensional and three-
dimensional simulations, we set the parameters of the simulator based on the
average percentage of white and gray matter. With good agreement, we have
compared the average velocity of the tumor growth resulting from the numerical
experiments with the velocity predicted by the theory. As future work we plan
analysis of the concurrency of the exponential integrators algorithm and
development of a parallel version of the simulator, targeting GPGPU computing
cards. Having the parallel simulator we will run large simulations on high-
resolution human brain data. Another future direction of our research is the
incorporation of the drug delivery term into our model [15]. In particular we
plan to focus on the chemotherapy modeling. We also plan to adress the data
assimilation problem using the supermodeling technique [17, 19].

## Acknowledgements

This work was partially supported by the program “Excellence initiative -
research university" for the AGH University of Kraków. Judit Muñoz-Matute has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie individual fellowship
grant agreement No. 101017984 (GEODPG).

## References

  * [1] Ahrens, J., Geveci, B. & Law, C. ParaView: An End-User Tool for Large Data Visualization. Visualization Handbook. (2005), ISBN 978-0123875822
  * [2] Al-Mohy, A. & Higham, N. Computing the action of the matrix exponential, with an application to exponential integrators. SIAM Journal On Scientific Computing. 33, 488-511 (2011)
  * [3] Cai, S., Mao, Z., Wang, Z., Yin, M. & Karniadakis, G. Physics-informed neural networks (PINNs) for fluid mechanics: A review. Acta Mechanica Sinica. 37, 1727-1738 (2021)
  * [4] Chen, Y., Lu, L., Karniadakis, G. & Dal Negro, L. Physics-informed neural networks for inverse problems in nano-optics and metamaterials. Optics Express. 28, 11618-11633 (2020)
  * [5] Croci, M. & Muñoz-Matute, J. Exploiting Kronecker structure in exponential integrators: Fast approximation of the action of $\phi$-functions of matrices via quadrature. Journal Of Computational Science. 67 pp. 101966 (2023), https://www.sciencedirect.com/science/article/pii/S1877750323000261
  * [6] Harpold, H., Alvord, E. & Swanson, K. The evolution of mathematical modeling of glioma proliferation and invasion. Journal Of Neuropathology And Experimental Neurology. 66, 1-9 (2007)
  * [7] Hertler, C., Felsberg, J., Gramatzki, D., Le Rhun, E., Clarke, J., Soffietti, R., Wick, W., Chinot, O., Ducray, F., Roth, P., McDonald, K., Hau, P., Hottinger, A., Jaap, R., Schnell, O., Marosi, C., Glantz, M., Darlix, A., Lombardi, G., Krex, D., Glas, M., Reardon, D., Bent, M., Lefranc, F., Herrlinger, U., Razis, E., Carpentier, A., Phillips, S., Rudà, R., Wick, A., Tabouret, E., Meyronet, D., Maurage, C., Rushing, E., Rapkins, R., Bumes, E., Hegi, M., Weyerbrock, A., Aregawi, D., Gonzalez-Gomez, C., Pellerino, A., Klein, M., Preusser, M., Bendszus, M., Golfinopoulos, V., Deimling, A., Gorlia, T., Y. Wen, P., Reifenberger, G. & Weller, M. Long-term survival with IDH wildtype glioblastoma: first results from the ETERNITY Brain Tumor Funders’ Collaborative Consortium (EORTC 1419). European Journal Of Cancer Research. 189 (2023)
  * [8] Hochbruck, M. & Ostermann, A. Explicit exponential Runge–Kutta methods for semilinear parabolic problems. SIAM Journal On Numerical Analysis. 43, 1069-1090 (2005)
  * [9] Hochbruck, M. & Ostermann, A. Exponential integrators. Acta Numerica. 19 pp. 209-286 (2010)
  * [10] Hogea, C., Davatzikos, C. & Biros, G. An image-driven parameter estimation problem for a reaction–diffusion glioma growth model with mass effects. Journal Of Mathematical Biology. 56 pp. 793-825 (2008)
  * [11] Kłusek, A., Łoś, M., Paszyński, M. & Dzwinel, W. Efficient model of tumor dynamics simulated in multi-GPU environment. The International Journal Of High Performance Computing Applications. 33, 489-506 (2019)
  * [12] Lipková, J., Angelikopoulos, P., Wu, S., Alberts, E., Wiestler, B., Diehl, C., Preibisch, C., Pyka, T., Combs, S., Hadjidoukas, P., Leemput, K., Koumoutsakos, P., Lowengrub, J. & Menze, B. Personalized Radiotherapy Design for Glioblastoma: Integrating Mathematical Tumor Models, Multimodal Scans, and Bayesian Inference. IEEE Transactions On Medical Imaging. 38 pp. 1875-1884 (2018)
  * [13] Maddalena, L. & Ragni, S. Existence of solutions and numerical approximation of a non-local tumor growth model. Mathematical Medicine And Biology: A Journal Of The IMA. 37, 58-82 (2020)
  * [14] Menze, B. H., Van Leemput, K., Honkela, A., Konukoglu, E., Weber, M.-A. , Ayache, N. & Golland, P. A generative approach for image-based modeling of tumor growth. In G. Székely and H. K. Hahn, editors, Information Processing in Medical Imaging, 735–747, Berlin, Heidelberg (2011)
  * [15] McDaniel, J., Kostelich, E., Kuang, Y., Nagy, J., Preul, M., Moore, N. & Matirosyan, N. Data Assimilation in Brain Tumor Models. Mathematical Methods And Models In Biomedicine. pp. 233-262 (2013)
  * [16] Muñoz-Matute, J. & Demkowicz, L. Multistage DPG time-marching scheme for nonlinear problems. ArXiv Preprint ArXiv:2309.00069. (2023)
  * [17] Paszyński, M., Siwik, L., Dzwinel, W. & Pingali, K. Supermodeling, a convergent data assimilation meta-procedure used in simulation of tumor progression. Computers & Mathematics With Applications. 113 pp. 214-224 (2022)
  * [18] Raissi, M., Perdikaris, P. & Karniadakis, G. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal Of Computational Physics. 378 pp. 686-707 (2019)
  * [19] Siwik, L., Łoś, M., Kłusek, A., Paszyńska, A., Pingali, K., Dzwinel, W. & Paszyński, M. Tuning three-dimensional tumor progression simulations on a cluster of GPGPUs. Journal Of Computational And Applied Mathematics. 412 pp. 114308 (2022)
  * [20] Stupp, R., Brada, M., Van den Bent, M., Tonn, J., Pentheroudakis, G. & ESMO Guidelines Working Group High-grade glioma: ESMO clinical practice guidelines for diagnosis, treatment and follow-up. Annals Of Oncology. 25 pp. 93-101 (2014,9)
  * [21] Swanson, K., Alvord Jr, E. & Murray, J. A quantitative model for differential motility of gliomas in grey and white matter. Cell Proliferation. 33, 317-329 (2000)
  * [22] Zhu, A. Accelerating Parameter Inference in Diffusion-Reaction Models of Glioblastoma Using Physics-Informed Neural Networks. (SIAM Undergraduate Resources Online,2022)
  * [23] Łoś, M., Kłusek, A., Hassaan, M., Pingali, K., Dzwinel, W. & Paszyński, M. Parallel fast isogeometric L2 projection solver with GALOIS system for 3D tumor growth simulations. Computer Methods In Applied Mechanics And Engineering. 343 pp. 1-22 (2019)
  * [24] Łoś, M., Paszyński, M., Kłusek, A. & Dzwinel, W. Application of fast isogeometric L2 projection solver for tumor growth simulations. Computer Methods In Applied Mechanics And Engineering. 316 pp. 1257-1269 (2017), Special Issue on Isogeometric Analysis: Progress and Challenges