Simulation of Spiking Neural P Systems with Sparse Matrix-Vector Operations

processes Article Simulation of Spiking Neural P Systems with Sparse Matrix-Vector Operations Miguel Ángel Martínez-del-Amor 1,2,* , David Orellana-Martín 1, Ignacio Pérez-Hurtado 1, Francis George C. Cabarle 3and Henry N. Adorna 3   Citation: Martínez-del-Amor, M.Á.; Orellana-Martín, D.; Pérez-Hurtado, I.; Cabarle, F.G.C.; Adorna, H.N. Simulation of Spiking Neural P Systems with Sparse Matrix-Vector Operations. Processes 2021,9, 690. https://doi.org/10.3390/ pr9040690 Academic Editors: Mengchu Zhou and Zhiwei Gao Received: 28 February 2021 Accepted: 12 April 2021 Published: 14 April 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: c 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1Research Group on Natural Computing, Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, 41012 Seville, Spain; [email protected] (D.O.-M.); per[email protected] (I.P.-H.) 2Smart Computer Systems Research and Engineering Lab (SCORE), Research Institute of Computer Engineering (I3US), Universidad de Sevilla, 41012 Seville, Spain 3Algorithms and Complexity Laboratory, Department of Computer Science, University of the Philippines Diliman, Quezon City 1101, Philippines; [email protected] (F.G.C.C.); [email protected] (H.N.A.) *Correspondence: [email protected] Abstract: To date, parallel simulation algorithms for spiking neural P (SNP) systems are based on a matrix representation. This way, the simulation is implemented with linear algebra operations, which can be easily parallelized on high performance computing platforms such as GPUs. Although it has been convenient for the first generation of GPU-based simulators, such as CuSNP, there are some bottlenecks to sort out. For example, the proposed matrix representations of SNP systems lead to very sparse matrices, where the majority of values are zero. It is known that sparse matrices can compromise the performance of algorithms since they involve a waste of memory and time. This problem has been extensively studied in the literature of parallel computing. In this paper, we analyze some of these ideas and apply them to represent some variants of SNP systems. We also provide a new simulation algorithm based on a novel compressed representation for sparse matrices. We also conclude which SNP system variant better suits our new compressed matrix representation. Keywords: spiking neural P systems; simulation algorithm; sparse matrix-vector operations; compressed matrix representation; GPU computing 1. Introduction Membrane computing [ 1 , 2 ] is an interdisciplinary research area in the intersection of computer science and cellular biology mainly [ 3 ], but also with many other fields such as engineering, neuroscience, systems biology, chemistry, etc. The aim is to study computational devices called P systems, taking inspiration from how living cells process information. Spiking neural P (SNP) systems [ 4 ] are a type of P system composed of a directed graph inspired by how neurons are interconnected by axons and synapses in the brain. Neurons communicate through spikes, and the time difference between them plays an important role in the computation. Therefore, this model belongs to the known third generation of artificial neural networks, i.e., based on spikes. Aside from computing numbers, SNP systems can also compute strings, and hence, languages. More general ways to provide the input or receive the output include the use of spike trains, i.e., a stream or sequence of spikes entering or leaving the system. Further results and details on computability, complexity, and applications of spiking neural P systems are detailed in [ 5 – 7 ], a dedicated chapter in the Handbook in [ 8 ], and an extensive bibliography until February 2016 in [9]. Moreover, there is a wide range of SNP system variants: with delays, with weights [ 10 ], with astrocytes [ 11 ], with anti-spikes [12] , dendrites [ 13 ], rules on synapses [ 14 ], scheduled synapses [ 15 ], stochastic firing [ 16 ], numerical [17], etc. Processes 2021,9, 690. https://doi.org/10.3390/pr9040690 https://www.mdpi.com/journal/processes Processes 2021,9, 690 2 of 30 The research on applications and variants of SNP systems has required the development of simulators. The simulation of SNP systems was initially carried out through sequential simulators such as pLinguaCore [ 18 ]. In 2010, a matrix representation of SNP systems was introduced [ 19 ]. Since then, most simulation algorithms are based on matrices and vector representations, and consists of a set of linear algebra operations. This way, parallel simulators can be efficiently implemented, since matrix-vector multiplications are easy to parallelize. Moreover, there are efficient algebra libraries that can be used out-of-the-box, although they have not been explored yet for this purpose. For instance, GPUs are parallel devices optimized for certain matrix operations [ 20 ], and can handle matrix operations efficiently. We can say, without loss of generality, that these matrix representations of SNP systems fit well to the highly parallel architecture of these devices. This have been harnessed already by introducing CuSNP, a set of simulators for SNP systems implemented with CUDA [ 21 – 24 ]. Simulators for specific solutions have been also defined in the literature [ 5 , 25 ]. Moreover, this is not unique for SNP systems, many simulators for other P system variants have been accelerated on GPUs [26–28]. However, this matrix representation can be sparse (i.e., having a majority of zero values) because the directed graph of SNP systems is not usually fully connected. A first approach to tackle this problem was presented in [ 29 ], where some of the ideas described in this work were described. Following these ideas, in [ 30 ], the transition matrix was split to reduce the memory footprint of the SNP representation. In many disciplines, sparse vector-matrix operations are very usual, and hence, many solutions based on compressed implementations have been proposed in the literature [31]. In this paper, we introduce compressed representations for the simulation of SNP systems based on sparse vector-matrix operations. First, we provide two approaches to compress the transition matrix for the simulation of SNP systems with static graph. Second, we extend these algorithms and data structures for SNP systems with dynamic graphs (division, budding, and plasticity). Finally, we make a complexity analysis and comparison of the algorithms to draw some conclusions. The paper is structured as follows: Section 2provides required concepts for the methods and algorithms here defined; Section 3defines the designs of the representations; Section 4 contains the detailed algorithms based on the compressed representations; Section 5 shows the results on complexity analyses of the algorithms; Section 6provides final conclusions, remarks, and plans of future work. 2. Preliminaries In this section we briefly introduce the concepts employed in this work. Firstly, we define the standard model of spiking neural P systems and three variants. Second, a matrix-based simulation algorithm for this model is revisited. Third, the fundamentals of compressed formats for sparse matrix-vector operations are given. 2.1. Spiking Neural P Systems Let us first formally introduce the definition of spiking neural P system. This model was first introduced in [4]. Definition 1. Aspiking neural P system of degree q ≥1is a tuple Π= (O,syn,σ1, . . . , σq,iout) where: •O={a}is the singleton alphabet (a is called spike); •syn ={(i , j)| ,1 ≤i , j≤q , i6=j} represents the arcs of a directed graph G= (V , syn) whose nodes are V ={1, . . . , q}; •σ1, . . . , σqare neurons of the form σi= (ni,Ri),1 ≤i≤q, Processes 2021,9, 690 3 of 30 where: –ni≥0is the initial number of spikes within neuron labeled by i; and –Riis a finite set of rules associated to neuron labeled by i, of one of the following forms: ( 1 )E/ac→ap, being E a regular expression over {a}, c ≥p≥1(firing rules); ( 2 )as→λ for some s≥ 1, with the restriction that for each rule E/ac→ap of type of type (1)from Ri, we have as6∈ L(E)(forgetting rules). •iout ∈ {1, 2, . . . , q}such that (iout,j)/∈syn, for any 1≤j≤q. A spiking neural P system of degree q≥ 1 can be viewed as a set of q neurons {σ1 , . . . , σq} interconnected by the arcs of a directed graph syn , called synapse graph. There is a distinguished neuron label iout , called output neuron ( σiout ), which communicates with the environment. If a neuron σi contains k spikes at an instant t , and ak∈L(E) , k≥c , then the rule E/ac→ap can be applied. By the application of that rule, c spikes are removed from neuron σi and the neuron fires producing p spikes immediately. Thus, each neuron σj such that (σi , σj)∈G receives p spikes. For σiout , the output neuron iout , the spikes are sent to the environment. The rules of type ( 2 ) are forgetting rules, and they are applied as follows: If neuron σi contains exactly s spikes, then the rule as→λ from Ri can be applied. By the application of this rule all sspikes are removed from σi. In spiking neural P systems, a global clock is assumed, marking the time for the whole system. Only one rule can be executed in each neuron at step t . As models of computation, spiking neural P systems are Turing complete, i.e., as powerful as Turing machines. On one hand, a common way to introduce the input (instance of the problem to solve) to the system is to encode it into some or all of the initial spikes ni ’s (inside each neuron i ). On the other hand, a common way to obtain the output is by observing neuron iout : either by getting the interval t2−t1=n , where σiout sent its first two spikes at times t1 and t2 (we say n is computed or generated by the system), or by counting all the spikes sent by σiout to the environment until the system halts. For the rest of the paper, we call this model spiking neural P systems with static structure, or just static SNP, given that the graph associated with it does not change along the computation. Next, we briefly introduce and focus on three variants with a dynamic graph: division, budding, and plasticity. A broader explanation of them and more variants are provided at [32–35]. Finally, let us introduce some notations and definitions: •pres(i): for a neuron σi, the presynapses of this neuron is pres(i) = {j|(i,j)∈syn}. •outdegree(i): for a neuron σi, the out degree of this neuron is: outdegree(i) = |pres(i)|. •ines(i): for a neuron σi, the insynapses of this neuron is ines(i) = {j|(j,i)∈syn}. •indegree(i): for a neuron σi, the in degree of this neuron is: outdegree(i) = |ines(i)|. 2.1.1. Spiking Neural P Systems with Budding Rules Based on the idea of neuronal budding, where a cell is divided into two new cells, we can abstract it to budding rules. In this process, the new neurons can differ in some aspects: their connections, contents, and shape. A budding rule has the following form: [E]i→[ ]i/[ ]j, where Eis a regular expression and i,j∈ {1, . . . , q}. If a neuron σi contains s spikes, as∈L(E) , and there is no neuron σj such that there exists a synapse (i , j) in the system, then the rule [E]i→[ ]i/[ ]j is enabled and it can be executed. A new neuron σj is created, and both neurons σi and σj are empty after the execution of the rule. This neuron σi keeps all the synapses that were going in, and this σj inherits all the synapses that were going out of σi in the previous configuration. There is Processes 2021,9, 690 4 of 30 also a synapse (i , j) between neurons σi and σj , and the rest of synapses of σj are given to the neuron depending on the synapses of syn. 2.1.2. Spiking Neural P Systems with Division Rules Inspired by the process of mitosis, division rules have been widely used within the field of membrane computing. In SNP systems, a division rule can be defined as follows: [E]i→[ ]j||[ ]k, where Eis a regular expression and i,j,k∈ {1, . . . , q}. If a neuron σi contains s spikes, as∈L(E) , and there is no neuron σg such that the synapse (g , i) or (i , g) exists in the system, g∈ {j , k} , then the rule [E]i→[ ]j||[ ]k is enabled and it can be executed. Neuron σi is then divided into two new cells, σj and σk . The new cells are empty at the time of their creation. The new neurons keep the synapses previously associated to the original neuron σi , that is, if there was a synapse from σg to σi , then a new synapse from σg to σj and a new one to σk are created, and if there was a synapse from σi to σg , then a new synapse from σj to σg and a new one from σk to σg are created. The rest of synapses of σjand σkare given by the ones defined in syn. 2.1.3. Spiking Neural P Systems with Plasticity Rules It is known that new synapses can be created in the brain thanks to the process of synaptogenesis. We can recreate this process in the framework of spiking neural P systems defining plasticity rules in the following form: E/ac→αk(i,Nj), where E is a regular expression, c≥ 1, α∈ {+ , − , ± , ∓} , k≥ 1 and Nj⊆ { 1, . . . , q} (a.k.a. neuron set). Recall that pres(i)is the set of presynapses of neuron σi. If a neuron σi contains s spikes, as∈L(E) , then the rule E/ac→αk(i , Nj) is enabled and can be executed. The rule consumes c spikes and, depending on the value of α , it performs one of the following: • If α= + and Nj−pres(i) = ∅ , or if α=− and pres(i) = ∅ , then there is nothing more to do. • If α= + and |Nj−pres(i)| ≤ k , deterministically create a synapse to every σg , g∈Nj−pres(i) . Otherwise, if |Nj−pres(i)|>k , then non-deterministically select k neurons in Nj−pres(i)and create one synapse to each selected neuron. • If α=− and |pres(i)| ≤ k , deterministically delete all synapses in pres(i) . Otherwise, if |pres(i)|>k , then non-deterministically select k neurons in pres(i) and delete each synapse to the selected neurons. • If α={± , ∓} , create (respectively, delete) synapses at time t and then delete (resp., create) synapses at time t+ 1. Even when this rule is applied, neurons are still open, that is, they can continue receiving spikes. Let us notice that if, for some σi , we apply a plasticity rule with α∈ {+ , ± , ∓} , when a synapse is created, a spike is sent from σi to the neuron that has been connected. That is, when σi attaches to σj through this method, we have immediately transferring one spike to σj. 2.2. Matrix Representation for SNP Systems Usually, parallel, P system simulators make use of ad-hoc representations, tailored for a certain variant [ 26 – 28 ]. In order to ease the simulation of static SNP system and its deployment to parallel environments, a matrix representation was introduced [ 19 ]. By using a set of algebraic operations, it is possible to reproduce the transitions of a computation. Although the baseline representation only involves SNP systems without delays and static structure, many extensions have followed such as for enabling delays [21,22] , Processes 2021,9, 690 5 of 30 handling non-determinism [ 24 ], plasticity rules [ 36 ], rules on synapses [ 37 ], and dendrite P systems [38] . In this section we briefly introduce the definitions for the matrix representation of the basic model of spiking neural P systems without delays, as defined above. We also provide the pseudocode to simulate just one computation of any P system of the variant using this matrix representation. In our notation, we use capital letters for vectors and matrices, and [] for accessing values: V[i] is the value at position i of the vector V , and M[i,j]is the value at row iand column jof matrix M. For an SNP system Π of degree (q , m) ( q neurons and m rules, where m=∑q i=1|Ri| ), we define the following vectors and matrices: Configuration vector : Ck is the vector containing all spikes in every neuron on the kth computation step/time, where C0 denotes the initial configuration; i.e., C0[i] = ni , for neuron σi= (ni,Ri). It contains qelements. Spiking vector : Sk shows if a rule is going to fire at the transition step k (having value 1 ) or not (having value 0). Given the non-determinism nature of SNP systems, it would be possible to have a set of valid spiking vectors, which is denoted as SVk . However, for the computation of the next configuration vector, only a spiking vector is used. It contains melements. Transition matrix:MΠis a matrix comprised of m·qelements given as M[i,j] =        −c, rule riis in σjand is applied consuming cspikes; p, rule riis in σs, with s6=jand (s,j)∈syn, and is applied producing pspikes; 0, rule riis in σswith s6=jand (s,j)/∈syn. In this representation, rows represent rules and columns represent neurons in the spiking transition matrix. Note also that a negative entry corresponds to the consumption of spikes. Thus, it is easy to observe that each row has exactly one negative entry, and each column has at least one negative entry [19]. Hence, to compute the transition k , it is enough to select a spiking vector Sk from all possibilities SVkand calculate: Ck=Sk·MΠ+Ck−1. The pseudocode to simulate a computation of an SNP system is as described in Algorithm 1. The selection of valid spiking vectors can be done in different ways, as in [21,22] . This returns a set of valid spiking vectors. In this work, we focus on just one computation, but non-determinism can be tackled by maintaining a queue of generated configurations [24]. Algorithm 1 MAIN PROCEDURE: simulating one computation for static spiking neural P systems. Require: A SNP system Πof degree (q,m), and a limit Lof time steps. Ensure: A computation of the system 1: (MΠ,C0)←INIT(Π) 2: k←0 3: repeat 4: SVk←SPIKING_VECTORS(Ck,Π).Calculate all possible spiking vectors 5: Sk←GET_ONE_RANDOMLY(SVk).Pick one spiking vector randomly 6: Ck+1←Ck+Sk·MΠ.Compute next configuration 7: k←k+1 8: until k≥L∨SVk=∅.Stop condition: maximum steps or no more applicable rules 9: return C0. . . Ck−1 In this work we focus on compressing the representation, specifically the transition matrix, so the determination of the spiking vector is not affecting these designs. Therefore, we use a straightforward approach and select just one valid spiking vector randomly. The representations here depicted only affect how the computation of the next configuration is done (matrix-vector multiplication at line 6in Algorithm 1). Processes 2021,9, 690 6 of 30 2.3. Sparse Matrix-Vector Operations Algebraic operations have been studied deeply in parallel computing solutions. Specifically, GPU computing provides large speedups when accelerating such kind of operations. This technology allows us to run scientific computations in parallel on the GPU, given that a GPU device typically contains thousands of cores and high memory bandwidth [ 39 ]. However, parallel computing on a GPU has more constraints than on a CPU: threads have to run in an SPMD fashion while accessing data in a coalesced way; that is, best performance is achieved when the execution of threads is synchronized and accessing contiguous data from memory. In fact, GPUs have been employed for P system simulations since the introduction of CUDA. Matrix computation is a highly optimized operation in CUDA [ 40 ], and there are many efficient libraries for algebra computations like cuBLAS. It is usual that when working with large matrices, these are almost “empty”, or with a majority of zero values. This is known as sparse matrix, and this downgrades the runtime in two ways: lot of memory is wasted, and lot of operations are redundant. Given the importance of linear algebra in many computational disciplines, sparse vector-matrix operations (SpMV) have been subject of study in parallel computing (and so, on GPUs). Today there exists many approaches to tackle this problem [ 41 ]. Let us focus on two formats to represent sparse matrices in a compressed way, assuming that threads access rows in parallel: • CSR format. Only non-null values are represented by using three arrays: row pointers, non-zero values, and columns (see Figure 1for an example). First, the row-pointers array is accessed, which contains a position per row of the original matrix. Each position says the index where the row start in the non-zero values and columns arrays. The non-zero values and the columns arrays can be seen as a single array of pairs, since every entry has to be accessed at the same time. Once a row is indexed, then a loop over the values in that row has to be performed, so that the corresponding column is found, and therefore, the value. If the column is not present, then the value is assumed to be zero, since this data structure contains all non-zero values. The main advantage is that it is a full-compressed format if NumNonZeroValues · 2 >NumZeroValues , where NumNonZeroValues and NumZeroValues are the number of non-zero and zero values in the original matrix, respectively. However, the drawbacks is that the search of elements in the non-zero values and columns arrays is not coalesced when using parallelism per row. Moreover, since it is a full-compressed format, there is no room for modifying the values, such as introducing new non-zero values; • ELL format. This representation aims at increasing the memory coalescing access of threads in CUDA. This is achieved by using a matrix of pairs, containing a trimmed, transposed version of the original matrix (see Figure 2for an example). Each column of the ELL matrix is devoted to each row of the matrix, even though the row is empty (all elements are zero). Every element is a pair, where the first position denotes the column and the second is the value, of only the non-zero elements in the corresponding row. However, the size of the matrix is fixed, so the number of columns equals the number of rows of the original matrix, but the number of rows is the maximum length of a row in terms of non-zero values; in other words, the maximum amount of non-zero elements in a row of the original matrix. Rows containing fewer elements pad the difference with null elements. The main advantage of this format is that threads always access the elements of all rows in coalesced way, and the null elements padded by short rows can be utilized to incorporate new data. However, there is a waste of memory, which is worst when the rows are unbalanced in terms of number of zeros. Processes 2021,9, 690 7 of 30 3 0 1 0 0 2 4 1 0 0 0 0 8 1 5 1 1 3 6 7 3 1 2 4 1 8 1 5 1 1 3 2 3 4 1 2 3 4 Row pointers: Non-zero val: Columns: Figure 1. CSR format example. Non-zero val array stores the non-null values, columns array stores the column indexes, and row pointers are the positions where each row starts in the previous arrays. 3 0 1 0 0 2 4 1 0 0 0 0 8 1 5 1 (1,3) (2,2) X (1,8) (3,1) (3,4) X (2,1) X (4,1) X (3,5) X X X (4,1) Maximum number of non-null values in a row of the original matrix Column Value transpose Figure 2. ELL format example. Note that it represents the transpose of the original matrix to increase the coalesced access in GPU devices. It includes pairs of column and value for every row. The compressed matrix has a number of columns equals to the number of original rows, and a number of rows equals to the maximum amount of non-null values in an original row. 3. Methods SNP systems in the literature typically do not contain fully connected graphs. This means that the transition matrix gets very sparse and, therefore, both computing time and memory are wasted. However, further optimizations based on SpMV can be conveyed. In the following subsections we discuss some approaches. Of course, if the graph inherent to an SNP system leads to a compressed transition matrix, then a normal (sparse) format can be employed, because using compressed formats will increase the memory footprint. In this work, we focus on the basic model of spiking neural P systems without delays as defined above, as well as three variants with dynamic network: budding, division and plasticity. The set of algorithms defined next are designed to take advantage of data parallelism, what is convenient for GPUs and vector co-processors. Their pseudocodes are detailed in Section 4of this paper. In Algorithm 2we generalize the pseudocode disposed in Algorithm 1to be able to handle both static and dynamic networks. This way, each variant requires to re-define just some functions from the ones defined for the static SNP system variant using vector and matrices without compression (i.e., sparse representation). In order to understand the algorithms, we will present in this section the main new data structures and their behavior. The detailed pseudocodes are available in Section 4. Algorithm 2 MAIN PROCEDURE: simulating one computation for spiking neural P systems. Require: An SNP system Πof degree (q,m), and a limit Lof time steps. Ensure: A computation of the system 1: (C0,M0)←INIT(Π) 2: k←0 3: repeat 4: SVk←SPIKING_VECTORS(Ck,Mk).Calculate all possible spiking vectors 5: Sk←GET_ONE_RANDOMLY(SVk).Pick one spiking vector randomly 6: (Ck+1,Mk+1)←COMPUTE_NEXT(Ck,Mk,SK).Compute next configuration 7: k←k+1 8: until k≥L∨SVk=∅.Stop condition: maximum steps or no more applicable rules 9: return C0. . . Ck−1 As a convention, those vectors and matrices using subindex k are dynamic and can change during the simulation time, while those with Π subindex are constructed at the Processes 2021,9, 690 8 of 30 beginning and are invariant. Capital letters refer to vectors and matrices, and small letters are scalar numbers. A total order over the rules defined in the system is assumed, which is denoted as R=Vq i=1Ri . For the sake of simplicity, we represent each rule rj≡Ej/acj→apj , with 1 ≤j≤m , as a tuple (i , Ej , cj , pj) , where i is the subindex of the set Ri where rj belongs (i.e., the neuron where it is contained). Specifically, forgetting rules just have pj=0. For static SNP systems using sparse representation, we use the following vectors and matrices: •Preconditions vector PΠ is a vector storing the preconditions of the rules; that is, both the regular expression and the consumed spikes. Initially, PΠ[j] = (Ej , cj) , for each rj∈R,rj= (ij,Ej,cj,pj), 1 ≤j≤m. •Neuron-rule map vector NΠ is a vector that maps each neuron index with its rules indexes. Specifically, NΠ[i] is the index of the first rule in the neuron. Given that rules have been ordered in R as mentioned above, rules belonging to the same neuron have contiguous indexes. Thus, it is enough to store just the first index. In this sense, the first rule in neuron i is NΠ[i] and the last one is NΠ[i+ 1 ]− 1. In other words, NΠ contains q+ 1 elements, and it is initialized as follows: NΠ[i] = ∑i−1 h=1 1 +|Rh| . Specifically, NΠ[1] = 1 and NΠ[i] = N0[i−1] + |Ri−1|, for 2 ≤i≤q+1 •Mk is the transition tuple , where Mk= (MΠ , PΠ , NΠ) . If the variant has a dynamic network, the transition matrix needs to be modified. Therefore, we start with M0 . The following algorithms show how they are constructed. Algorithm 2can be easily transformed to Algorithm 1by defining the INIT and COMPUTE_NEXT functions as in Algorithm 3. They work exactly as already specified in Section 2.2; that is, using the usual vector-matrix multiplication operation to calculate the next configuration vector. We will also detail how the selection of spiking vectors can be done. This is defined in Algorithm 4, and it is based on previous ideas already presented in [21,22] . First, SPIKING_VECTORS function calculates the set of all possible spiking vectors by using a recursive function over neuron index i . It gathers all spiking vectors that can be generated for neurons i0>i and then. If neuron i contains applicable rules, it populates a spiking vector for each of these rules, and from each of the generated spiking vectors form neurons i0>i . Finally, neuron i propagates these spiking vectors to the next neuron i−1. 3.1. Approach with ELL Format Our first approach to compress the representation of the transition matrix, MΠ , is to use the ELL format (see Figure 3for an example). The reason for using ELL and not other compressed formats for sparse matrices (CSR, COO, BSR, ...) is to enable extensions for dynamic networks, as seen later. ELL can give some room for modifications without much memory re-allocations, while CSR requires us to modify the whole matrix to add new elements. ELL format leads to the new compressed matrix Ms Π . The following aspects have been taken into consideration: • The ELL format represents the transpose of the original matrix, so now rows correspond to neurons and columns to rules. This is convenient for SIMD processors such as GPUs. • The number of rows of Ms Π equals the maximum amount of non-zero values in a row of MΠ , denoted by z0 . It can be shown that z0=z+ 1, where z is the maximum output degree found in the neurons of the SNP system. Specifically, z=max{outdegree(i)| 1 ≤ i≤q} (see definition in Section 2.1). z0 can be derived from the composition of the transition matrix, where row j devoted for rule rj≡(ij , Ej , cj , pj) contains the values +pj for every neuron i (columns) connected though an output synapse with the neuron where the rule belongs to (i.e., i∈pres(ij) ), and a value −cj for consuming the spikes in the neuron the rule belongs to (i.e., ij). Processes 2021,9, 690 9 of 30 • The values inside columns can be sorted, so that the consumption of spikes ( −c values) are placed at the first row. In this way, if implemented in parallel, all threads can start by doing the same task: consuming spikes. •Every position of Ms Πis a pair (as illustrated in Figure 3), where the first element is a neuron label, and the second is the number of spikes produced (+p). A parallel code can be implemented with this design by assigning a thread to each rule, and so, one per column of the spiking vector Sk and one per column of Ms Π (rows of the original transition matrix). For the vector-matrix multiplication, it is enough to have a loop of z0 steps at maximum through the columns. In the loop of each column j , the corresponding value in the spiking vector Sk[j] (either 0 or 1) is multiplied to the value xj in the pair (ij , xj) , and added to the neuron id nj in the configuration vector Ck[nj] . In case the SNP network contains hubs (nodes with high amount of input synapses), then we can opt for a parallel reduction per column. Since some threads might write to same positions in the configuration vector at the same time, a solution would be to use atomic adding operations, which are available on devices such as GPUs. Algorithm 3 Functions for static SNP systems with sparse matrix representation. 1: procedure INIT(Π) 2: (C0,NΠ)←INIT_NEURON_VECTORS(Π).Initialize vectors only related to neurons. 3: (MΠ,PΠ)←INIT_RULE_MATRICES(Π).Initialize matrices related to rules. 4: M0←(MΠ,PΠ,NΠ) 5: return (C0,M0) 6: end procedure 7: procedure INIT_NEURON_VECTORS(Π) 8: (q,σ1, . . . , σq)←Π.Get information from Π 9: C0←EMPTY_VECTOR(q).Create initial configuration 10: NΠ←EMPTY_VECTOR(q+1) .Create neuron-rule vector 11: NΠ[1]←1 12: for all i←1 . . . qdo .For each neuron 13: (ni,Ri)←σi.Get info of the neuron from Π 14: C0[i]←ni.Initial configuration 15: NΠ[i+1]←NΠ[i] + |Ri|.Neuron-rule map vector initialization 16: end for 17: return (C0,NΠ) 18: end procedure 19: procedure INIT_RULE_MATRICES(Π) 20: (R,m,q,pres)←Π.Get information from Π 21: PΠ←EMPTY_VECTOR(m).Create preconditions vector 22: MΠ←EMPTY_MATRIX(m,q).Create transition matrix 23: for all rj∈R,j←1 . . . mdo .For each rule (column). This loop is parallelizable. 24: rj≡(ij,Ej,cj,pj).Get info of the rule 25: PΠ[j]←(Ej,cj).Store it in precondition vector 26: MΠ[j,ij]← −cj.Construct transition matrix 27: for all i∈pres(ij)do .For each connected neuron to ij 28: MΠ[j,i]←pj.Construct transition matrix 29: end for 30: end for 31: return (MΠ,PΠ) 32: end procedure 33: procedure COMPUTE_NEXT(Ck,Mk,SK) 34: (MΠ,_,_)←Mk.Get some content of transition tuple 35: return (Ck+Sk·MΠ,Mk).Only the configuration is updated. 36: end procedure Processes 2021,9, 690 16 of 30 for optimized matrix representation, is slightly modified (just two lines) to support the neuron-id map vector. • APPLICABLE function as in Algorithm 8. This function, when dealing with division rules, has to search if there are existing synapses for the neurons involved. If they exist, the division rule does not apply. • COMPUTE_NEXT function as in Algorithm 9, to include the operations described above. It now needs to expand the synapse matrix Sykeither by columns (when new neurons are created) or by rows if there is a neuron from which we need to create a synapse to the new neuron and it has already the maximum out degree z . In this case, we need to re-allocate the synapse matrix in order to extend it by one row (this is written in the pseudocode with the function EXPAND_MATRIX). Finally, let us remark that we can easily detect if type of a rule rj at it’s associated cj value: if 0, it is a budding rule, if it is positive number, a spiking rule, otherwise (negative value) it is a division rule. Algorithm 7 Initialization functions for dynamic SNP systems with budding and division rules over optimized compressed matrix representation. 1: procedure INIT(Π) 2: (C0,NΠ)←INIT_NEURON_VECTORS(Π).Initialize vectors only related to neurons. 3: (M0 0,Q0,PΠ)←INIT_RULE_MATRICES(Π).Initialize matrices related to rules. 4: M0←(M0 0,Q0,PΠ,NΠ) 5: return (C0,M0) 6: end procedure 7: procedure INIT_NEURON_VECTORS(Π) 8: (q,qmax,σ1, . . . , σq)←Π.Get information from Π 9: C0←EMPTY_VECTOR(qmax).Create initial configuration 10: NΠ←EMPTY_VECTOR(q+1) .Create neuron-rule vector 11: NΠ[1]←1 12: for all i←1 . . . qdo .For each neuron 13: (ni,Ri)←σi.Get info of the neuron from Π 14: C0[i]←ni.Initial configuration 15: NΠ[i+1]←NΠ[i] + |Ri|.Neuron-rule map vector initialization 16: end for 17: return (C0,NΠ) 18: end procedure 19: procedure INIT_RULE_MATRICES(Π) 20: (R,m,q,qmax,z,pres)←Π.Get information from Π 21: PΠ←EMPTY_VECTOR(m).Create preconditions vector 22: RuΠ←EMPTY_VECTOR(m) 23: for all rj∈R(j←1 . . . m)do .For each rule (column). This loop is parallelizable. 24: rj≡(ij,Ej,cj,pj).Get info of the rule 25: PΠ[j]←(Ej,cj).Store it in precondition vector 26: RuΠ[j]←(cj,pj).Store it in rule vector 27: end for 28: Q0←EMPTY_VECTOR(qmax) 29: Sy0←EMPTY_MATRIX(z,qmax) 30: for i←1 . . . qdo .For each neuron label 31: Q0[i]←i 32: k←1.kis our iterator, compacting the rows 33: for all h∈pres(i)do .For each out synapse 34: Sy0[k,i]←h 35: k←k+1.We know that k≤z 36: end for 37: end for 38: M0 0←(RuΠ,Sy0) 39: return (M0 0,Q0,PΠ) 40: end procedure Processes 2021,9, 690 17 of 30 Algorithm 8 Applicable functions for dynamic SNP systems with budding rules over optimized compressed matrix representation. 1: procedure APPLICABLE(i,j,Ck,Mk) 2: (M0 k,Qk,PΠ,_)←Mk.Get some content of transition tuple 3: (RuΠ,Syk)←M0 k.Get some content of transition matrix 4: (lj,hj)←RuΠ[j].Preconditions of the rule 5: (Ej,cj)←PΠ[j].Preconditions of the rule 6: if cj>0then .If a spiking or forgetting rule 7: return Ck[i]∈L(Ej)∧Ck[i]≥cj. If rule j is applicable in neuron icj= 0 . If a budding rule 8: b←False .Check if synapse (i,hj)exists 9: w←1 10: while ¬b∧Syk[w,i]6=0∧w≤zdo .Until an empty value or the maximum 11: b←Qk[Syk[w,i]] = hj.If synapse exists 12: w←w+1 13: end while 14: return Ck[i]∈L(Ej)∧ ¬b 15: else .If a division rule 16: b←False .Check if synapse (i,hj)or (i,−lj)exists 17: w←1 18: while ¬b∧Syk[w,i]6=0∧w≤zdo .Until an empty value 19: b←Qk[Syk[w,i]] = hj∨Qk[Syk[w,i]] = −lj.If either synapse exists 20: w←w+1 21: end while 22: for x←1 . . . qkdo .Search for neurons with label hjor lj 23: b0←Qk[x] = hj∨Qk[x] = −lj.Either hjor lj 24: w←1 25: while b0∧ ¬b∧Syk[w,x]6=0∧w≤zdo .Combination of conditions 26: b←Qk[Syk[w,x]] = i.If the synapse exists 27: w←w+1 28: end while 29: end for 30: return Ck[i]∈L(Ej)∧ ¬b 31: end if 32: end procedure 3.3.2. Plasticity Rules For dynamic SNP systems with plasticity rules, the synapse matrix can be allocated in advance to the exact size q , since no new neurons are created. Thus, there is no need of using a neuron-id map vector as before. However, enough rows (value z ) in the synapse matrix have to be pre-established to support the maximum amount of synapses. Fortunately, this can be pre-computed by looking to the initial out degrees of the neurons and the size of the neuron sets in the plasticity rules adding synapses. We encode a plasticity rule rj≡Ej/ac j→αjkj(ij , Nj) , with αj= +/−/±/∓ as follows: rj= (ij , Ej , cj , αj , kj , Nj) . Next, we define the value of zp for SNP systems with plasticity rules: zp=max{|pres(i)∪ Nti| ,1 ≤i≤q} , where Nti=Sj=m j=1Nj , rj∈Ri , rj= (ij , Ej , cj , αj , kj , Nj) , αj∈ {+ , ± , ∓} . In other words, zp is the maximum out degree ( z ) that a neuron can have initially plus those new connections that can be created with plasticity rules inside that neuron. This result can be refined for plasticity rules having α∈ {± , ∓} , because we know up to k new synapses can be created at a time. However, for simplicity, we will use the formula above. Processes 2021,9, 690 18 of 30 Algorithm 9 Compute next function for dynamic SNP systems with budding and division rules using optimized compressed matrix representation. 1: procedure COMPUTE_NEXT(Ck,Mk,Sk) 2: Ck+1←Ck.Create a copy of Ck 3: (M0 k,Qk,_, _)←Mk.Extract info from transition tuple 4: (RuΠ,Syk)←M0 k.Extract info from transition matrix 5: qk← |Ck|.Get current amount of neurons. 6: for i←1 . . . qkdo .For each neuron. This loop is parallelizable. 7: j←Sk[i].Index of rule to fire in the neuron 8: if j6=0then .Only if there is a rule. 9: (cj,pj)←RuΠ[j].Get rule info 10: if cj>0then .Execution of a spiking or forgetting rule 11: Ck+1[i]←Ck+1[i]−cj.Consume spikes in firing neuron 12: w←1.Next while stops if pj=0, i.e., a firing rule 13: while pj>0∧Syk[w,i]6=0∧w≤zdo .Until an empty value 14: h←Syk[w,i].Get connected neuron by a synapse 15: Ck+1[h]←Ck+1[h] + pj.Produce spikes in connected neuron 16: w←w+1 17: end while 18: else if cj=0then .Execution of a budding rule 19: qk←qk+1.Increment counter of neurons 20: Ck+1[i]←0.Empty the neuron i, neuron qkis 0 already 21: Qk[qk]←pj.pjis the label of the new neuron 22: for w←1 . . . zdo 23: Syk[w,qk]←Syk[w,i].Copy column ito the new one 24: if w=1then .Update out synapses of i 25: Syk[w,i]←qk.The only new out synapse of i 26: else 27: Syk[w,i]←0.No more out synapses for i 28: end if 29: end for 30: else .Execution of a division rule 31: (hj,lj)←(cj,−pj).Get new neurons labels 32: qk←qk+1.Increment counter of neurons 33: Ck+1[i]←0.Empty the neuron i, neuron qkis 0 already 34: Qk[i]←hj.The new label of the neuron 35: Qk[qk]←lj.The label of the new neuron 36: for w←1 . . . zdo .Copy out synapses to new neuron 37: Syk[w,qk]←Syk[w,i].Copy column ito the new one 38: end for 39: for x←1 . . . qk−1do .Search for in synapses of neuron i 40: b←False .Boolean saying the synapse was found 41: w←1 42: while Syk[w,x]6=0∧w≤zdo .Search the end of the column 43: b←b∨Syk[w,x] = i.Search for neuron i 44: w←w+1 45: end while 46: if bthen .If synapse was found, add new neuron at the end of the column 47: if w=zthen 48: z←z+1.The neuron xhas a larger out degree than z 49: Syk←EXPAND_MATRIX(Syk,z,qk).Extend with one more row 50: end if 51: Syk[w,x]←qk.This can lead to overflows if w=z 52: end if 53: end for 54: end if 55: end if 56: end for 57: return (Ck+1,Mk) 58: end procedure Processes 2021,9, 690 19 of 30 First, we need to represent plasticity rules into vectors. We assume that np is the total amount of plasticity rules in the system, and that there is a total order between these rules. Given a plasticity rule, we can initialize the neuron-map and the precondition vector as with spiking rules. But in this case, we need a couple of new vectors and modify existing ones in order to represent all plasticity rules rj= (ij , Ej , cj , αj , kj , Nj) , with j∈ { 1 . . . np} (following the imposed total order): • Rule vector RuΠ stores the following pair for a plasticity rule rj : (cj , −j) , that is, the consumed spikes cj and the unique index of the plasticity rule j . This index is used to access the following vector, and it is stored as a negative value in order to detect that this is a plasticity rule. •Plasticity rule vector PrΠ , of size np , contains a tuple for each plasticity rule rj of the form (αj , kj , nij , nej) . The values nij and nej are used as indexes, from nij (start) to nej (end) for the following vector. •Plasticity neuron vector PnΠ , of size npr =∑j=np j=1|Nj| , represents all neuron sets of plasticity rules. Thus, the elements of Nj are stored, in an ordered way, between PnΠ[nij]to PnΠ[nej]. •Time vector Tk is used to prevent neurons from applying rules during one step if the plasticity rule applied was of the type αj∈ {±,∓}. It contains binary (0 or 1) values. • Transition matrix is therefore M0 k= (RuΠ , Syk , Tk , PrΠ , PnΠ) . Note that the Synapse matrix can be modified at each step, so we use sub-index k. The following operations have to be performed to reproduce the behavior plasticity rules (see Figure 7for an illustration): 1. For each column in the synapse matrix executing a plasticity rule deleting x synapses: (a) If the intersection of the rule’s neuron set and the current synapses in Syk is larger than x, then randomly select xsynapses. (b) Loop through the rows (up to zp iterations) to search the selected neurons and set them to null. Given that holes might appear in the column, its values can be sorted (or compacted). 2. For each column in the synapse matrix executing a plasticity rule adding x synapses: (a) If the difference of the rule’s neuron set and the current synapses in Syk is larger than x, then randomly select xneurons. (b) Loop through the rows (up to zp iterations) to insert the selected new synapses while keeping the order. Checking the applicability of plasticity rules is much simpler than for division rules, given that the preconditions only affect to the local neuron and we do not need to know if there are existing synapses. However, for a plasticity rule r in a neuron i , and in order to create new or delete existing synapses, we need to check which neurons declared in r are already in the column i in the synapse matrix. This search can be O(zp·nr) , being nr the length of the neuron set in r . Nevertheless, by maintaining always the order in the column, this search can be done easily in O(zp+nr). Given that it is not usual to have budding and division rules together with plasticity rules, the pseudocode is based on the optimized matrix representation for static SNP systems (and not for division and budding) in Section 3.2. Algorithm 10 shows the redefinition of INIT_RULE_MATRICES and COMPUTE_NEXT functions, replacing those from Algorithm 3. For COMPUTE_NEXT, the implementation is very similar to the original one, but it just call to a new function, PLASTICITY, which actually modify the synapses of the neuron (by just modifying its corresponding column in the synapse matrix). This function and its auxiliaries are defined in Algorithms 11 and 12, respectively. Processes 2021,9, 690 20 of 30 I4I2 I3 I4 I5 I6 I1I1 I4 I1 I5 I6 I1 α=+,k=3,i=2,N={Ik=3,k=3,i=2,N={Ii=2,k=3,i=2,N={IN={I2,k=3,i=2,N={II3,k=3,i=2,N={II5}I2 I2 I4 I3 I5 zp α=-,k=3,i=2,N={Ik=2,k=3,i=2,N={Ii=4,k=3,i=2,N={IN={I2,k=3,i=2,N={II3,k=3,i=2,N={II4} I3 I4 I1 I2 I3 I5 I5 I6 I1 I3 sort sort Figure 7. Illustration of application of a plasticity rule in the synapse matrix with compressed representation. Light blue cells in the synapse matrix are empty values (0), dark cells are positions with values greater than 0 (i.e., with neuron label). Two examples are give, in case of adding new synapses ( top ) and in case of deleting synapses ( bottom ). We sort the synapses per column for more efficiency. 4. Algorithms In this section we define the algorithms implementing the methods described in Section 3. Let us first define a generic function to create a new, empty (all values to 0) vector of size s as follows: EMPTY_VECTOR( s ). In order to create an empty matrix with f rows and c columns, we will use the following function: EMPTY_MATRIX( f , c ). Next, the pseudocodes for simulating static SNP systems with sparse representation are given. Algorithm 3shows the INIT and COMPUTE_NEXT functions, while Algorithm 4shows the selection of spiking vectors. For ELL-based matrix representation for static SNP systems , we need to re-define only two functions (INIT_RULE_MATRICES and COMPUTE_NEXT) from Algorithm 3 (static SNP systems with sparse representation) as shown in Algorithm 5. For our optimized matrix representation for static SNP systems , we need to redefine only two functions (INIT_RULE_MATRICES and COMPUTE_NEXT) from Algorithm 3 (static SNP systems with sparse representation) as shown in Algorithm 6 . Moreover, the following two lines in the definition of COMBINATIONS function at Algorithm 4are required, in order to support a spiking vector of size q: •Line 13 at Algorithm 4:S←EMPTY_VECTOR(q) •Line 22 at Algorithm 4:S0[i]←j For dynamic SNP systems with budding and division rules , the following functions are redefined: INIT functions as in Algorithm 7, APPLICABLE function as in Algorithm 8, and COMPUTE_NEXT function as in Algorithm 9. The SPIKING_VECTORS function, as defined in Algorithm 4and modified in Section 3.2 for optimized matrix representation, is slightly modified (just two lines) to support the neuron-id map vector as follows: •Line 6at Algorithm 4:(_, Qk,PΠ,NΠ)←Mk •Line 18 at Algorithm 4:for j←NΠ[Qk[i]] . . . NΠ[Qk[i] + 1]−1 Processes 2021,9, 690 21 of 30 Algorithm 10 Functions for dynamic SNP systems with plasticity rules using optimized compressed matrix representation. 1: procedure INIT_RULE_MATRICES(Π) 2: (R,m,q,zp,np,npr,pres)←Π.Get information from Π,zpis different for plasticity. 3: PΠ←EMPTY_VECTOR(m).Create preconditions vector 4: RuΠ←EMPTY_VECTOR(m).Create rule vector 5: PrΠ←EMPTY_VECTOR(np).Create plasticity rule vector 6: PnΠ←EMPTY_VECTOR(npr).Create plasticity neuron vector 7: pk ←1.Counter for plasticity rule vector 8: ni ←1.Counter for plasticity neuron vector 9: for all rj∈R,j←1 . . . mdo .For each rule (column). This loop is parallelizable. 10: if rj= (ij,Ej,cj,pj)then .If spiking rule 11: PΠ[j]←(Ej,cj).Store it in precondition vector 12: RuΠ[j]←(cj,pj).Store it in rule vector 13: else if rj= (ij,Ej,cj,αj,kj,Nj)then .If plasticity rule 14: PΠ[j]←(Ej,cj).Store it in precondition vector 15: RuΠ[j]←(cj,−pk).Store it in rule vector 16: PrΠ[pk]←(αj,kj,ni,ni +|Nj|).Store it in rule vector 17: PnΠ[ni]←SORT(Nj).Sort and store Njin Pn after position ni 18: pk ←pk +1 19: ni ←ni +|Nj| 20: end if 21: end for 22: T0←EMPTY_VECTOR(q).Create time vector 23: Sy0←EMPTY_MATRIX(zp,q).Create synapse matrix 24: for all i←1 . . . qdo .For each neuron (column in synapse matrix) 25: k←1.kis our iterator, compacting the rows 26: for all h∈pres(i)do .For each out synapse 27: Sy0[k,i]←h 28: k←k+1.We know that k≤zp 29: end for 30: end for 31: M0 0←(RuΠ,Sy0,T0,PrΠ,PnΠ).New transition matrix 32: return (M0 0,PΠ) 33: end procedure 34: procedure COMPUTE_NEXT(Ck,Mk,Sk) 35: Ck+1←Ck.Create a copy of Ck 36: (M0 k,PΠ,NΠ)←Mk.Get some content of transition tuple 37: (RuΠ,Syk,PrΠ,PnΠ)←M0 k 38: for i←1 . . . qdo .For each neuron. This loop is parallelizable. 39: j←Sk[i].Index of rule to fire in the neuron 40: if j6=0∨Tk[i] = 1then .Only if there is a rule or blocked neuron 41: (cj,pj)←RuΠ[j].Get rule info 42: Ck+1[i]←Ck+1[i]−cj.Consume spikes in firing neuron 43: if pj>0then .If a spiking rule 44: w←1 45: while SyΠ[w,i]6=0∧w≤zpdo .Until an empty value or the maximum 46: h←SyΠ[w,i].Get connected neuron by a synapse 47: Ck+1[h]←Ck+1[h] + pj.Produce spikes in connected neuron 48: w←w+1 49: end while 50: else if pj<0then .If a plasticity rule 51: (A,t)←PLASTICITY(Syk[,i],−pj,PrΠ,PnΠ).Modify only column i 52: Syk[,i]←A Tk[i]←t 53: end if 54: end if 55: Tk[i]←0.Reset time vector 56: end for 57: M0 k+1←(RuΠ,Syk,Tk,PrΠ,PnΠ).Next transition matrix 58: return (Ck+1,(M0 Π,PΠ,NΠ) 59: end procedure Processes 2021,9, 690 22 of 30 For dynamic SNP systems with plasticity rules , the pseudocode is based on the optimized matrix representation for static SNP systems (and not for division and budding) in Section 3.2. Algorithm 10 shows the re-definition of INIT_RULE_MATRICES and COMPUTE_NEXT functions, replacing those from Algorithm 3. As for line 17, we assume that the function SORT exists, which takes a set of neurons, sorts them by id, and generates a vector. Moreover, we can copy vectors directly from one position by just one assignation. The new PLASTICITY function is defined in Algorithm 11, and its auxiliaries are defined in Algorithm 12. Algorithm 11 Function for plasticity mechanism using optimized compressed matrix representation. 1: procedure PLASTICITY(Ai,j,PrΠ,PnΠ) 2: (αj,kj,nij,nej)gets PrΠ[j].Get info of plasticity rule 3: nprj←nej−nij.Number of neurons in the neuron set 4: Nj←EMPTY_VECTOR(nprj).Create a vector with the neuron set 5: for x←1 . . . nprjdo 6: Nj[x]←PnΠ[nij+x].Copy the contents of the neuron set 7: end for 8: t←0.Vale for time vector, only 1 for mp,pm 9: if αj=−then .Delete synapses 10: Ai←DEL_SYNAPSES(Ai,kj,Nj) 11: else if αj= + then .Add synapses 12: Ai←ADD_SYNAPSES(Ai,kj,Nj) 13: else if αj=±then .Add and delete synapses 14: Ai←ADD_SYNAPSES(Ai,kj,Nj) 15: Ai←DEL_SYNAPSES(Ai,kj,Nj) 16: t←1 17: else if αj=∓then .Delete and add synapses 18: Ai←DEL_SYNAPSES(Ai,kj,Nj) 19: Ai←ADD_SYNAPSES(Ai,kj,Nj) 20: t←1 21: end if 22: return (Ai,t) 23: end procedure In order to keep Algorithm 12 simple, we assume that the functions INTERSEC, DIFF, and DELETE_RANDOM are already defined. As mentioned above, INTERSEC and DIFF can be implemented with algorithms of complexity O(zp+nr) , given that the vectors (a column of synapse matrix and a chunk of plasticity neuron vector) are already sorted. We also assume that DELETE_RANDOM is a function that randomly select k elements from a total of n while keeping the order between elements. This can be done with an algorithm of complexity O(k2). 5. Results In this section we conduct a complexity analysis (for both time and memory) of the algorithms. In order to define the formulas, we need to introduce a set of descriptors for a spiking neural P system Π . These are described in Table 1. Moreover, Table 2 summarizes the vectors and matrices employed by each representation, and their corresponding sizes defined according to the descriptors. We use the following short names for the representations: Sparse (original sparse representation as Section 3), ELL (ELL compressed representation as in Section 3.1), optimized static (optimized static compressed representation as in Section 3.2), division and budding (optimized dynamic compressed representation for division and budding as in Section 3.3.1), and plasticity (optimized dynamic compressed representation for plasticity as in Section 3.3.2). Processes 2021,9, 690 23 of 30 Algorithm 12 Auxiliary functions for plasticity mechanism using optimized compressed matrix representation. 1: procedure DEL_SYNAPSES(A,k,N) 2: N0←INTERSEC(A,N).Calculate A∩N(involved synapses to be deleted) 3: if |N0|>kthen .If more than kneurons, select randomly 4: N0←DELETE_RANDOM(N0,|N0| − k).A random set of kneurons 5: end if 6: k← |N0|.The new amount of synapses to delete 7: w←p←s←1.Initialize iterators 8: while w≤zp∧A[w]6=0do .Loop over the column 9: if A[w] = N0[s]then .Synapse to delete 10: A[w]←0.Delete the synapse 11: s←s+1.Advance in N0vector 12: else 13: if p<wthen .Need to compact the vector 14: A[p]←A[w].pis the last compacted position 15: A[w]←0 16: end if 17: p←p+1.Advance the last compacted position p 18: end if 19: w←w+1 20: end while 21: return A 22: end procedure 23: procedure ADD_SYNAPSES(A,k,N,n) 24: N0←DIFF(N,A).Calculate N\A(not involved synapses to create) 25: if |N0|>kthen .If more than kneurons, select randomly 26: N0←DELETE_RANDOM(N0,|N0| − k).A random set of kneurons 27: end if 28: k← |N0|.The new amount of synapses to delete 29: B←EMPTY_VECTOR(zp).Create the output 30: w←p←s←1.Initialize iterators 31: while w≤zp∧ ¬(A[w] = 0∧s≤k)do .Loop over the column 32: if A[w]>N0[s]then .Synapse to add 33: B[p]←N0[s].Add the synapse 34: s←s+1 35: else 36: B[p]←A[w].Keep the synapse 37: w←w+1 38: end if 39: p←p+1 40: end while 41: return B 42: end procedure Table 1. Descriptors of an SNP system. Descriptor Description qNumber of initial neurons mTotal number of rules zNumber of rows (column size) for optimized matrix representation. Also, maximum out degree of a neuron z0Number of rows (column size) for ELL matrix representation. z0=z+1 qmax Maximum amount of neurons to handle during simulation for division and budding. q0≥q zpNumber of rows (column size) for optimized matrix representation for plasticity npMaximum size of a neuron set in plasticity rules. npr Sum of neuron set sizes of all plasticity rules. kpMaximum value of kin plasticity rules. Processes 2021,9, 690 24 of 30 Table 2. Size of matrices employed in the representations for SNP systems. Those whose name are in bold were used for the total calculation, which assumes just one spiking vector. Notation/Name Sparse ELL Optimized Static Division and Budding Plasticity CkConfiguration Vector q q q qmax q SkSpiking vector m m q qmax q PΠPreconditions vector 2m2m2m2m2m NΠNeuron-rule map vector q+1q+1q+1q+1q+1 SykSynapse matrix q·z qmax ·zq·zp RuΠRule vector 2m2m2m QkNeuron-id map vector qmax PrΠRule vector 4m PnΠRule vector npr TkRule vector q M0 ΠTransition matrix m·q2·m·z0q·z+2m qmax ·z+2mq(zp+1) + 6m+npr MkTransition tuple m·q+2m+q+1m(2z0+2) + q+1q(z+1) + 4m+1qmax (z+1) + 4m+q+1q(zp+2) + 8m+npr +1 TOTAL m·q+3m+2q+1m(2z0+3) + 2q+1q(z+3) + 4m+1qmax (z+2) + 4m+2q+1q(zp+4) + 8m+npr +1 According to Table 2, we can limit the value of qmax for dynamic SNP systems with division and budding with the following formula: qmax =bMT−4m+2q+1 z+2c , where MT is the maximum amount of memory in the system (measured in the word size employed to encode the elements of all the matrices and vectors; e.g., 4 Bytes). Moreover, we can infer when the matrix representation will be smaller for static SNP systems: ELL is better than sparse when z<q−2 2 ; optimized is better than ELL when z>q−m 2m−q ; optimized is better than sparse when z<m− 2. In other words, our optimized compressed representation is worth when the the maximum out degree of the neurons is less than the total number of rules minus 2. For dynamic SNP systems, given can say that a solution to a problem using an SNP with plasticity rules is better than a solution based on division and budding, if qmax >q(zp+2)+4m+npr z+2 ; in other words, if we can know the maximum amount of neurons to generate, and this number is greater than a formula based on number of initial neurons, number of rules, and number of elements in the neuron set and the max out degree, then the solution will need less memory using plasticity. Finally, Table 3shows the order of complexity of each function as defined for each representation. We can see that COMPUTE_NEXT gets reduced in complexity as well when using optimized static representation against ELL and sparse, given that we expect that m≤q , and also z<m− 2. However, we can see that implementing division and budding explodes the complexity of the algorithms, since they need to loop over all the neurons checking for in-synapses. This also depends on the total amount of generated neurons in a given step. This is also the case for the generation of spiking vectors, because the applicability function also needs to loop over all existing neurons. However, for dynamic networks, plasticity keeps the complexity with the amount of neurons, the value z , and the descriptors of plasticity rules (max value of kand amount of neurons in a neuron set np). Table 3. Algorithmic order of complexity of the main functions employed in the simulation loop (i.e., excluding init functions) for each representation. Function Sparse ELL Optimized Static Division & Budding Plasticity APPLICABLE O(1)O(1)O(1)O(qmax ·z)O(1) SPIKING_VECTORS O(m)O(m)O(m)O(m·qmax ·z)O(m) PLASTICITY O(np+k2 p+zp) COMPUTE_NEXT O(q·m)O(z0·m)O(z·q)O(qmax ·z)O(q·(np+k2 p+zp)) Processes 2021,9, 690 25 of 30 Therefore, we can see that using our compressed representations, both the memory footprint of the simulators and their complexity are reduced, as long as the maximum out degree of neurons is a low number. Furthermore, we can see that for dynamic networks, plasticity is an option that keeps the complexity balanced, since we know in advance the amount of neurons and synapses. Let us make an easy example of comparison with an example from the literature. For example, if we take the SNP system for sorting natural numbers as defined in [ 42 ], then we have that q= 3 n , m=n+n2 and z=n , where n is the amount of natural numbers to sort. Thus: • The size of the sparse representation is 3 n3+ 6 n2+ 5 n+ 1 and the complexity of COMPUTE_NEXT is O(n3). • The size of the ELL representation is 2 n3+ 7 n2+ 11 n+ 1 and the complexity of COMPUTE_NEXT is O(n3). • The size of the optimized representation is 7 n2+ 13 n+ 1 and the complexity of COMPUTE_NEXT is O(n2). The optimized representation drastically decreases the order of complexity and amount of memory spent for the algorithms, going from orders of n3 to n2 . ELL has a similar order of complexity to that of sparse, but the amount of memory is just a bit decreased. Figure 8shows that the reduction of the memory footprint achieved with the compressed representations takes effect after n> 3. Figure 9shows that the optimized representation scales better than ELL and sparse. ELL is only a bit better than the sparse representation, demonstrating the need for using the optimized one, which significantly scales much better. Finally, we also analyze a uniform solution to 3SAT with SNP systems without delays as in [ 43 ] (Figure 10). We can see that q= 8 n3+ 3 n+ 3, m= 64 n3+ 6 n+ 3 and z= 8 n3 , where n is the amount of variables in the 3SAT instance. We can see that z<m− 3, so our optimized implementation will be able to save some memory. Therefore: • The size of the sparse representation is 512 n6+ 240 n4+ 424 n3+ 18 n2+ 51 n+ 25 and the complexity of COMPUTE_NEXT is O(n6). • The size of the ELL representation is 1024 n6+ 96 n4+ 384 n3+ 36 n+ 22 and the complexity of COMPUTE_NEXT is O(n6). • The size of the optimized representation is 64 n6+ 24 n4+ 304 n3+ 33 n+ 22 and the complexity of COMPUTE_NEXT is O(n6). Figure 8. Memory size of the matrix representation (Y-axis) depending on the amount of natural numbers ( n , X-axis) for the model of natural number sorting in [ 42 ], using sparse, ELL and optimized representation, for only n=1,2, 3,4.