GPU implementation of evolving spiking neural P systems

GPU implementation of evolving spiking neural P systems Rogelio V. Gungon a , Katreen Kyle M. Hernandez a , Francis George C. Cabarle a , Ren Tristan A. de la Cruz a , Henry N. Adorna a , Miguel Á. Martínez-del-Amor b,c, ⇑ , David Orellana-Martín b,c , Ignacio Pérez-Hurtado b,c a Department of Computer Science, University of the Philippines Diliman, Quezon City 1101, Philippines b Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, Seville 41012, Spain c Smart Computer Systems Research and Engineering Lab (SCORE), Research Institute of Computer Engineering (I3US), Universidad de Sevilla, Seville 41012, Spain article info Article history: Received 9 August 2021 Revised 30 May 2022 Accepted 27 June 2022 Available online 30 June 2022 Communicated by Zidong Wang Keywords: Membrane computing Spiking neural P systems Genetic algorithm Evolutionary computing GPU computing CUDA abstract Methods for optimizing and evolving spiking neural P systems (in short, SN P systems) have been previously developed with the use of a genetic algorithm framework. So far, these computations, both evolving and simulating, were done only sequentially. Due to the non-deterministic and parallel nature of SN P systems, it is natural to harness parallel processors in implementing its evolution and simulation. In this work, a parallel framework for the evolution of SN P Systems is presented. This is the result of extending our previous work by implementing it on a CUDA-enabled graphics processing unit and adapting CuSNP design in simulations. Using binary addition and binary subtraction with 3 different categories each as initial SN P systems, the GPU-based evolution runs up to 9x faster with respect to its CPU-based evolution counterparts. Overall, when considering the whole process, the GPU framework is up to 3 times faster than the CPU version. Ó2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction Evolution is evident in nature. In biology, cells display features of self-healing, adaptation, and learning through evolution in a parallel manner, with a degree of non-determinism [1]. The paradigm of adapting and evolving to reach approximate solutions that nature presents brought life to natural computing. This paper is an interplay between two branches of natural computing: Membrane computing explores new computational models inspired by biological cells, specifically cellular membranes. This branch deals with distributed and parallel computing models that perform calculations using a biologically-inspired process. The process of computation is abstracted from chemical reactions when objects inside the cell cross cell membranes. The various types of membrane systems are known as P systems after Gheorghe Pa ˘un first formulated the model in [2]. Evolutionary computing covers optimization techniques inspired by principles of biological evolution. This branch situates an initial population of candidate solutions in an environment where less desired solutions are more likely to be removed. As a result, the population will gradually evolve to increase its fitness. Incorporating two different research fields is a method used by researchers in an attempt to discover new solutions to problems. As such, the interactions between membrane computing and evolutionary computing, called evolutionary membrane computing (EMC) has been explored in [3]. Three new research topics have emerged from the study [4]: membrane-inspired evolutionary algorithms (MIEAs), automated design of membrane computing models (ADMCM), and membrane evolutionary algorithms (MEA). So far, designing a P system for performing a specific task has been a challenge for researchers. P systems are carefully designed by experts and could not be automatically created by programs. Automatically designing a P system by using programsis an ongoing research direction under ADMCM [5–7]. From these works, it has been shown that afeasible way to do this is touse evolutionary algorithms to evolve a population of P systems until it is well designed. In this paper, a new framework which can be categorized under ADMCM is proposed. We use one of the main paradigms of evolutionary computing called genetic algorithm (GA) to evolve spiking neural P systems or SN P systems. An SN P System is a https://doi.org/10.1016/j.neucom.2022.06.094 0925-2312/Ó2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author at: Departamento de Ciencias de la Computación e Inteligencia Artificial, E.T.S. Ingeniería Informática, Avda. Reina Mercedes S/N, 41012 Sevilla, Spain. E-mail addresses: [email protected] (K.K.M. Hernandez), fccabarle@up. edu.ph (F.G.C. Cabarle), [email protected] (R.T.A. de la Cruz), [email protected]. ph (H.N. Adorna), [email protected] (M.Á. Martínez-del-Amor), [email protected] (D. Orellana-Martín), [email protected] (I. Pérez-Hurtado). Neurocomputing 503 (2022) 140–161 Contents lists available at ScienceDirect Neurocomputing journal homepage: www.elsevier.com/locate/neucom neural-like P system that is introduced in [8], where neurons and synapses are represented as nodes of a graph connected by edges. The objects present in the neurons, called spikes, are handled by specific sets of rules. This is discussed further in Section 2.1. In evolving SN P Systems, it is necessary to develop ways to help researchers compute, analyze, and extract results faster. Evolution has to be as efficient and fast as possible to handle complex and large-sized instances. Because of its nature as a population-based search method, GA has great potential to improve its performance by utilizing parallel computations [9]. Moreover, the inherent parallelism of SN P systems also leads us to looking for simulators developed over highly parallel platforms [10–15]. The main reason for using Graphics Processing Units (GPUs) for this purpose is that they offer a cheap platform with a high degree of data parallelism. GPUs include thousands of computing cores that can be used for general purpose computations. CUDA is a technology that exposee the GPU as a highly parallel processor [16]. This opens the avenue for massively parallel computing. The CUDA platform is designed to work with multiple programming languages such as C and Python [17]. This accessibility makes it easier for researchers to use GPU resources. In previous work [18], SN P Systems are evolved and simulated sequentially. In this paper, we aim to produce a more time-efficient system by maximizing the highly parallel nature of SN P Systems on a CUDA-enabled GPU. This was done by introducing the following improvements: 1. The GA operators, mutation and crossover, are implemented to run partially parallel on a CUDA-enabled GPU by using threads to facilitate arithmetic processes such as synapse modification. This is further discussed in Section 3.3. 2. The framework in [18] simulates the systems in the population sequentially. In this paper, the framework is modified so that it can simulate all SN P systems in the population in parallel. 3. The CuSNP design is adapted from [15] to use threads in modifying neuron states of an SN P system in order to speed up individual simulations. 4. The aforementioned developments were able to accomplish an evolution speedup of up to 9x with respect to the framework in [18]. Overall, considering also data transfers, the system is up to 3x faster. Similarworks totheone wepresent herealso appearin [19,20].A recent study [20] developed a simulator, a representation for SN P systems with polarization (PSN P system), and a GA framework for evolving them. PSN P systems are a variant of SN P systems in which the charge of the neuron determines whether a rule can be applied. The framework modifies the synapses and rules to evolve the PSN P systems. The framework was tested on two functions, AND and OR, over three categories: Perfect, Extra, and Mutated. Another recent study [19] proposed a genetic algorithm framework for evolving Spiking Neural P systems where the framework is limited to modifying only the synapses and rules. The framework was tested on two initial SN P systems, i.e., bitwise ADD and SUB, over three types. Each type has varying fitness, number of neurons, and number of rules. [18–20]eachpresentedageneticalgorithmframeworkforevolving their respective SN P system variant. They designed their framework to run sequentially on a CPU. In contrast, this work aims to accelerate the GA framework in [18] using CUDA-enabled GPU. Using GPU, this work evolves and simulates SN P systems in parallel using the canonical GAoperators in the framework from [18].Like in [18],theGAframeworkinthispaperis constrainedtoevolveadeterministic initial SN P system only through its topology. The rules inside the neurons must stay constant, while the synapses and neurons may vary. This is different from [19,20] where the framework is limited only to modifying the synapsesand the rules of the given SN P system while the number of neurons is fixed. This paper is organized as follows: In Section 2, we provide the preliminaries for the study. Section 3discusses the chromosome representation of SN P systems, GA framework, and GA function. The experimentation setup and results are discussed in Sections 4 and 5, respectively. Lastly, conclusions and direction of future works are discussed in Section 6. 2. Preliminaries In this section, SN P systems are formally defined, together with a matrix representation of the model. SubSections 2.1 and 2.2 are devoted to discussing them respectively. The matrix representation of an SN P system is the basis for the simulation algorithm used in this paper. The nature of SN P systems allowed their representation as matrices, which is a crucial step in simulating them on highly parallel devices such as GPUs. Using a CUDA-enabled GPU in simulating SN P systems have been studied in related works [12,11,13]; this paper will primarily focus on CuSNP [14,15] as its simulation design is the one adapted in SubSection 3.3.2. CuSNP is given a brief introduction in subSection 2.5. Some preliminaries about Genetic Algorithmare discussed in 2.3 and 2.4. 2.1. Spiking Neural P Systems Pa ˘un in [8] introduced the SN P systems [21] in a general form, namely, an extended (i.e., with the rules able to produce more than one spike) computing (i.e., able to take an input and provide an output) version. It consists of a set of neurons placed in the nodes of a directed graph and sending signals (spikes, denoted by the symbol a) along synapses (graph edges). It is further defined as follows: Definition 1 (SN P System). A spiking neural P system, of degree mP1, is a construct of the form P ¼ðO; r 1 ;...; r m ;syn;in;outÞ; where: 1. O¼fagis the singleton alphabet (ais called spike); 2. r 1 ;...; r m are neurons, of the form r i ¼ðn i ;R i Þ;16i6m, where: (a) n i P0 is the initial number of spikes contained in r i ; (b) R i is a finite set of rules of the following two forms: (1) (Spiking Rule) E=a c !a p ;d, where Eis a regular expression over aand cPpP1;dP0; (2) (Forgetting Rule) a s !k, for sP1, with the restriction that for each rule E=a c !a p ;dof type (1) from R i , we have a s RLðEÞ; 3. syn #f1;2;...;mgf1;2;...;mgwith i–jfor all ði;jÞ2syn;16i;j6m(synapses between neurons); 4. in;out 2f1;2;...;mgindicate the input and output neurons, respectively. The rules of type (1) are called spiking rules, while those of type (2) are called forgetting rules. In the case of LðEÞ¼fa c g, spiking and forgetting rules can be shortened to a c !a p and a c !k, respectively; this is applied to the rules in Fig. 1. In each time unit, A rule E=a c !a b 2R i must be applied if neuron r i contains kspikes, where a k 2LðEÞand kPc.IfpP1, applying the rule consumes or removes cspikes from r i , so that kcspikes remain in r i . The neuron sends pspikes to every r j such that ði;jÞ2syn. The output neuron does not have a synapse directed to any other neuron, R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 141 only to the environment. The moments of time when a spike is emitted by the output neuron are marked with 1, the other moments are marked with 0. The binary sequence obtained in this way is called the output spike train of the system. A computation in a system starts with the initial configuration. The initial configuration of a system is a record of the spikes each neuron holds before the feeding on input spikes or spike trains. In order to compute a function f:N k !N, we introduce knatural numbers n 1 ;...;n k in the system by reading from the environment a binary sequence z¼10 n 1 1 10 n 2 1 1...10 n k 1 1. This means that the input neuron of P receives a spike in each step corresponding to a digit 1 from the string zand no spike otherwise. Note that we input exactly kþ1 spikes, i.e., after the last spike we assume that no further spike is coming to the input neuron. 2.2. Matrix Representation Matrix manipulation plays a significant role in evolution and simulation of SN P systems. Matrix algorithms are numerically intensive computations, these are commonly implemented in parallel [22]. We can offload these computations to GPU. In this section, we discuss a matrix representation for SN P systems. Matrix representation of SN P systems without delay was first introduced in [23].In[14], modifications were made for the representation to include delays. This representation is used for the numerical computations in simulating SN P systems on a GPU. Let P be an SN P system with mneurons and nrules. We use the following definitions from [14] to represent P . Definition 2 (Configuration Vector). The configuration vector C ðkÞ ¼hc 1 ;...;c m i, where c i is the amount of spikes in r i at time k. Definition 3 (Spiking Vector). A spiking vector S ðkÞ ¼hs 1 ;...;s n i where s ðkÞ i ¼1ifE i is satisfied and r i is applied; 0 otherwise  Definition 4 (Status Vector). The kth status vector is denoted by St ðkÞ ¼hst 1 ;...;st m iwhere for each i2f1;2;...;mg, st i ¼1 if neuron mis open; 0 if neuron mis closed  Definition 5 (Rule Representation). Rules are represented by R ðkÞ ¼hr 1 ;...;r n iwhere for each i2f1;2;...;ng, a rule is represented by vector r i ¼ðE;j;d 0 ;cÞwhere Eis the regular expression for rule i;jis the neuron that contains the rule r i , d 0 ¼ 1 if the rule is not fired; 0 if the rule is fired; P1 if the rule is currently on delay 8 > < > : and cis the number of spikes in neuron r j consumed if r i is applied. Definition 6 (Delay Vector). The delay vector D ¼hd 1 ;...;d n icontains the delay count for each rule r i ;i2f1;2;...;ngin P . Definition 7 (Loss Vector). The loss vector LV ðkÞ ¼hl v 1 ;...;l v m i where for i2f1;2;...;mg;l v i is the number of consumed spikes in r i at the step k. Definition 8 (Gain Vector). The kth gain vector is denoted by GV ðkÞ ¼hg v 1 ;...;g v m iwhere for i2f1;2;...;mg;g v i is the number of spikes sent by connected neurons to neuron r i at the step k. Definition 9 (Transition Matrix). The transition matrix 1 of Pis an ordered set of vectors TV defined as TV ¼ftv1;...;tvngwhere for each i2f1;2;...;ng;tvi¼hp1;...;pmisuch that if ri2rswhere p j ¼P0 number of spikes produced by r i given ðs;jÞ2 syn; 0 otherwise  Definition 10 (Indicator Vector). The indicator vector IV k ¼hi v 1 ;...;i v m iindicated which rule will produce spikes at time k. Definition 11 (Removing Matrix). The removing matrix is denoted by RM ¼frm 1 ;...;rm n gwhere for each i2f1;2;...;ng; rm i ¼ht i ;...;t m isuch that if r i 2 r s where t j ¼P0 number of spikes consumed by r i given s¼j; 0 otherwise  Definition 12 (Net Gain Vector). The net gain vector at step kis denoted by NG ðkÞ ¼C ðkþ1Þ C ðkÞ . 2.3. Genetic Algorithm Solutions to optimization problems are represented by chromosomes. Genetic algorithm (in short, GA) is used to search through a pool of chromosomes, in which we call population, an approximate solution to optimization problems. There are different genetic operators as introduced in [24] but we will be using mutation, selection, and crossover. A modified canonical genetic algorithm is presented in Fig. 2 from [18]. The evolution starts by generating a pool of candidate solutions, which we call a population comprising of chromosomes. These chromosomes are mutations generated from an original solution or system. A fitness function defined prior to the start of the evolution is then used to assess the fitness of the chromosome. These fitness values are used in the selection process in which a subset of chromosomes from the population is chosen to proceed to the next process. This subset of chromosomes now comprise the new population. The population undergoes the crossover process in which pairs Fig. 1. An example of an SN P system. It simulates the binary addition function. r 1 and r 2 are the input neurons, while r 5 is the output neuron. See Table 1 for. a sample run. 1 Transition matrix differs from the spiking transition matrix in [23], the former only contains spikes gained by each neuron while the latter contains both spikes consumed and produced by each neuron at step k. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 142 of chromosomes are randomly chosen to become parents and create offspring. The offspring will either replace their parent chromosomes in the population or are simply appended. The population then undergoes mutation in which bits of randomly chosen chromosomes are randomly flipped. This is done to introduce new genetic material to the next population to escape a local minimum [18]. One iteration of evolution is called a generation. The entire process is repeated until one of the predefined halting conditions is met: in this case, convergence condition and maximum number of generations. 2.4. Fitness and Precision We will be using the definition of the terms fitness and precision as used in [18].Precision is the output of an evaluation function, which computes the output for each chromosome independently from its peers. The evaluation function works by simulating a set of input spike trains on the chromosome and compares its outputs to the corresponding desired outputs, using longest common substring (LCS) rate to determine its precision. This is further discussed in Section 3.3.2.Fitness of a chromosome is computed by getting its precision with respect to its peer chromosomes’ precision in the population. 2.5. CuSNP Given that this implementation is designed to run on a CUDA- enabled GPU, to label this implementation as parallel would be a misnomer. It is only partially parallel due to the computing architecture of CUDA where the GPU works with the CPU. The CPU is considered as the master node which calls for the execution of multiple compute kernels on the GPU when needed [26]. In this case, the CPU is also used to prepare the variables for computing and control the flow of evolution. Different functions are defined in the kernel which are called by the CPU to execute. As with the current technologies involving CuSNP [14,15], the parallel implementation involved assigning threads to neurons or rules depending on the kernel function called. Aside from facilitating the flow of compute kernel execution and evolution, the CPU also initializes, pre-processes, and postprocesses the data [14]. Pre-processing of data includes converting it from P-Lingua format to its chromosome representation (13), extracting and preparing specific arrays from the available data that will fit the needs of the called kernel function, and copying the data from the host memory to the device memory. Postprocessing includes gathering the processed data from device memory to host memory, printing the data for storage and further inspection, and formatting the chromosomes to P-Lingua format. 3. Evolution Design Genetic algorithm is used to evolve an initial SN P system into one with more desirable properties: in this case, topology and precision. SubSection 3.2 discusses the Genetic Algorithmframework which tells how the evolution is facilitated. SubSection 3.3 discusses the Genetic AlgorithmFunction or how evolution is executed. Chromosomes are used to represent a solution, i.e. an SN P system, within a population. An SN P system P can be encoded into a chromosome, formally defined in SubSection 3.1. 3.1. Chromosome Representation SN P Systems are represented by chromosomes as stated in Section 2.3. A chromosome encoding was introduced by [18] which is defined as follows. Definition 13 (Chromosome). An SN P system P of mneurons is represented by chromosome chrom, a construct of the form chrom= (num_nrns, spikes, rules, input_nrns, output_nrn, syn_matrix, precision, fitness, active_nrns, connected_nrns, connected_to_in- put_nrns, connected_to_output_nrns), where: 1. num nrns ¼mþ1, the environment is also counted as a neuron; 2. spikes ¼fn 0 1 ;...;n 0 mþ1 g, where n 0 i P0 is the number of initial spikes contained in r i , where r i is a neuron in P ; 3. rules ¼fR 0 1 ;...;R 0 mþ1 g, where R 0 i is a finite set of rules contained in r i of the same form as R i ; 4. input nrns ¼fi 1 ;...;i jinj g, where r i is an input neuron, i2jinj; 5. output nrn ¼i, where r i is the output neuron; 6. syn matrix ¼½a ij  mm , where a ij ¼1 if there exists a synapse from r i to r j ; 0 otherwise  7. precision ¼p, where 0 6p61; 8. fitness ¼f, where 0 6f61; 9. acti v e nrns ¼fi 1 ;...;i k gwith 0 6k6ðmþ1Þ, where r i is an active neuron; 10. connected nrns ¼fi 1 ;...;i k gwith 0 6k6ðmþ1Þ, where r i is a connected neuron; 11. connected to input nrns ¼fi 1 ;...;i k gwith 0 6k6ðmþ1Þ, where there exists a path from at least one of the input neurons to r i ; 12. connected to output nrns ¼fi 1 ;...;i k gwith 0 6k6ðmþ1Þ, where there exists a path from r i to the environment neuron. Neurons have two classifications: active and connected. A neuron is considered active if it is within a path from an input neuron to the environment. If not, a neuron can still be considered active given that it is connected to the SN P System by at least one synapse, it has at least one initial spike, and it has a corresponding rule to use. These conditions ensure that the neuron have a possibility of affecting the resulting spike train. If these conditions are not met, the neuron has no effect on the resulting spike train is considered inactive. Fig. 2. Canonical genetic algorithm flowchart. This is a slightly modified version of the binary genetic algorithm flowchart in [25], wherein we replaced ‘‘cost”, ‘‘select mates”, and ‘‘mating” with ‘‘fitness”, ‘‘selection”, and ‘‘crossover”, respectively [18]. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 143 A neuron is considered connected if there exists a path connecting it to the environment or at least one input neuron. Otherwise, the neuron is disconnected. 3.2. Genetic AlgorithmFramework In this section, we discuss the main framework of the evolution algorithm presented in Fig. 3. We have adapted the work of [18]. This framework is designed to transform an initial SN P system P init that approximates the function fða;b;c;...Þ¼zinto a smaller and more precise system P final that approximates the same output zgiven the same input/s fa;b;c;...g. The SN P systems are encoded in P-Lingua format [27]. The framework starts by transforming P init from P-Lingua format into chromosome chrom init ; this is denoted by the star in Fig. 3. The GA framework contains two prominent loops: the inner loop and outer loop. The inner loop represents the fundamental characteristic of genetic algorithms wherein it searches the best solution from a generated pool of candidate solutions, which we know to be the population. The best solution is defined as the SN P System that has the highest precision or smaller topology in comparison to the initial SN P System. The variable iter counts the number of iterations made inside the inner loop. The inner loop iterates if the GA function failed to produce a chrom final better than chrom init in terms of topology and precision, and iter has not reached the user-defined maximum number of iterations an inner loop can undergo which is denoted by variable q . There are two ways of exiting the loop. First is when proceeding to start a new outer loop; this is triggered when the produced chrom final is better than chrom init in terms of topology. Second is when halting altogether which is triggered when chrom final has a precision greater than or equal to chrom init ’s precision, or when the maximum number of iterations have been made. The outer loop automates the process of continuous evolution wherein the output of the inner loop, denoted by chrom final is used as the new chrom init . The implementation of the GA function’s Selection operator ensures that if its produced chrom final is not equal to its initial input chrom init ;chrom final is a better SN P System compared to chrom init . 3.3. Genetic AlgorithmFunction The GA function is illustrated in Fig. 4 which shows that it runs partially sequentially and partially parallel, these are further discussed in the following subsections. The function takes in three different sets of input. The first input is the chrom init which is the initial chromosome or SN P System to evolve. Given the set of fitness cases S¼fða 0 ;b 0 Þ;ða 1 ;b 1 Þ;...;ða k ;b k gwhere kis the number of fitness cases and kP1;a i are tuples consisting of input spike trains f a 0 ; a 1 gwhile b i is the corresponding ideal output spike train. From this, we can obtain the two other sets of input which are the set of input spike trains denoted by A¼fa 0 ;a 1 ;...;a k g and the set of ideal ouput spike trains denoted by B¼fb 0 ;b 1 ;...;b k g. The function takes chrom init as the initial input and forms the initial population, population 0 ¼fchrom 0 ;chrom 1 ;...;chrom p gwith pas the defined population size, which consists of pnumber of chrom init mutations. The initial population is then passed to the SN P System simulator along with Awhich produces a set of output spike trains B 0 . The sets Band B 0 are passed through the evaluation function to compute the precision of each chromosome of the population by comparing the value of b 0 i to the corresponding ideal output spike train b i . The population is sorted by precision from greatest to least. Based on the defined selection rate, a subset of chromosomes from the population are selected and the rest are discarded. Before the new population undergoes crossover, it first identifies if there exists a chromosome fit to be a chrom final . The selection process and conditions to be met in choosing a fit chrom final are further discussed in Section 3.3.2. Otherwise, it proceeds to the crossover process where pairs of chromosomes are randomly chosen as parents to produce offspring chromosomes which either replace the parents or are added to the current population, provided that the size of the population will still be within or exactly p. This represents one cycle of evolution or a generation. The evolution cycles until the maximum number of generations is made or the selection process selects a chrom final . 3.3.1. Mutation In the first generation, the population is built by mutating the chrom init ptimes. For the succeeding generations, each chromosome of the population has a chance of being mutated according to the user defined mutation parameter of range ½0;1. If a chromosome is to be mutated, one mutation method is chosen to executed on the chromosome. There are three methods which have equal chances of being chosen: 1. Disconnecting a neuron. This method is done by randomly choosing a neuron excluding input neurons, the environment neuron, and disconnected neurons. Synapses are created to connect the chosen neuron’s pre-synaptic neurons of its in-going synapses to all post-synaptic neurons of its outgoing synapses. The chosen neuron’s synapses are then deleted. 2. Deleting a synapse. This method is done by randomly choosing a neuron. From its synapses, one is randomly chosen to be deleted. A self loop can arise after execution of this method. 3. Adding a synapse. This method is done by randomly choosing a neuron excluding the environment neuron. This condition is to ensure that the environment neuron does not supply spikes to a neuron in the system. A synapse is then added to connect it to another randomly chosen neuron. Disconnected neurons are possibly reconnected to the system after execution of this method. As with deleting a synapse, a self loop can arise after execution of this method. The mutation process partially runs sequentially and in parallel. The preparation of variables in the form of arrays are done in the CPU while the execution of the methods are done in the GPU. There are three variables prepared: 1. A matrix of form ½syn matrix 0 ;...;syn matrix c with cbeing the number of chromosomes randomly chosen to be mutated with 06c6pwith respect to the mutation rate parameter, given that pis the number of chromosomes in the population; 2. An array of form ½task 0 ;...;task c containing tasks to be performed on each chromosome with mnumber of neurons, including the environment neuron. Tasks are denoted by task i where for each i2f0;1;...;cg, taski¼ 06j6miif neuron rjis to be disconnected; 0>jPðmiþ1Þif a synapse from rjis to be added or deleted 8 > < > : This shows that positive values denote that the method to disconnect a neuron is executed, and that negative values denote that the method to add or delete a synapse is to be executed. This design combines the mutation method decision process and choosing neurons. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 144 3. An array of form ½syn 0 ;...;syn c containing synapses to be added or deleted. These synapses are denoted by syn i 2f0;...;mgnfegwhere i2f0;1;...;cgand r e is the environment neuron. Due to the nature of synapse matrices, the addition or deletion of a synapse can be made by flipping the corresponding bit. These are copied from host to device for the parallel segment of the mutation process. In the device, there are pðmmÞthreads where ðmmÞrepresents a synapse matrix. The task array is checked per chromosome con which method to execute. If task c P0, the neuron r task c is disconnected from the SN P System. If task c <0;a ij from the syn matrix c is flipped, with i¼ðtask c 1Þ and j¼syn c .Ifa ij ’s value is flipped from 0 to 1, a synapse is added from r i to r j . On the other hand, if a ij ’s value is flipped from 1 to 0, the existing synapse from r i to r j is deleted. The output of the parallel mutation execution is an updated list of the population’s synapse matrices which is copied from device back to host. On CPU, the chromosomes are updated with its corresponding new synapse matrix. 3.3.2. Simulation, Evaluation Function and Selection In order to obtain set B 0 in Fig. 4, the chromosomes in the population must be simulated with a set of input spike train pairs or A. The simulation process runs partially parallel. As with the mutation process, the preparation of variables are done in the CPU and the simulation execution is done in the GPU. Algorithm1 gives an overview of the implementation of the simulator in pseudocode form. The algorithm for the actual simulation used in step 5 is from [14] (Algorithm1). Algorithm1 Overview of SN P system simulation algorithm. SEQ and PAR indicate which step is done sequentially or in parallel, respectively. Input: population g ¼fchrom 0 ;chrom 1 ;...;chrom p g, A¼fa 0 ;a 1 ;...;a k gwhere population g is the population of chromosomes at generation gwith pas the defined population size and Ais the set of input spike trains where k is the number of fitness cases. Output: B0¼fb0 0 ;b0 1 ;...;b0 k gwhere B0is the set of output spike trains. 1: (SEQ) Parse each chrom i from population g ; 2: (SEQ) Load parsed data from population g into variables; 3: (SEQ) Initialize and allocate memory for vectors; 4: (SEQ) Transfer data from host to device for simulation; 5: (PAR) Perform simulation on each chrom i for all input spike trains a i The longest common substring (LCS) [28] rate between the ideal and actual output spike trains, b i and b 0 i respectively in Fig. 4,is used to evaluate the precision of each chromosome. The longest common substring is defined as the longest string that is a substring of both ideal and actual output spike trains. The LCS rate can be computed by dividing the length of the longest common substring by the length of the ideal output spike train. This method is used so that computation errors from padding in the actual output spike trains can be avoided. Fig. 3. Baseline Genetic AlgorithmFramework. qis a user-defined parameter which defines the number of maximum iterations of inner loops within one iteration of the outer loop. Fig. 4. Genetic AlgorithmFunction. Modified from [18] to have certain sections run in a parallel manner. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 145 The selection operator first checks for convergence in every generation. Initially, chrom init is set as chrom best . The selection operator halts the generation when the following conditions are met: 1. there exists a chromosome in the current population with precision greater than chrom best 2. there exists a chromosome in the current population with the same precision as chrom best but with fewer neurons 3. the maximum number of generations is met Otherwise, the chromosomes are sorted from highest to lowest precision and a subset starting from the highest ranking chromosomes are selected to remain in the population that will proceed to the crossover process. The subset is chosen based on the selection rate. If the rate were set to 1:00, all chromosomes in the current population remain. This design was adapted from [15] which sequentially traverses each bit of the input spike train while simulating in parallel the firing of spikes or changes in matrices. In order simulate all fitness cases which have varying input spike train lengths, a maximum value MAX SIZE was used in order to make the lengths equal. MAX SIZE was also used to pad the matrices involved in simulation in order to maintain thread and block allocation uniformity. This uniformity is helpful in referencing the correct threads and blocks that handled the specific matrix indices in each fitness case and each chromosome of the population. 3.3.3. Crossover The crossover process runs partially sequentially and in parallel. The list of parent chromosome pairs Pand corresponding neurons eligible for swapping Sare prepared in the CPU. The process of preparing the list of parent chromosomes starts with creating two copies of the current population defined as pop 0 0 and pop 0 1 with the chromosomes shuffled. In the event that the size of the current population is less than the set population size, the difference between the values is taken and defined as d. For every chromosome chrom 0;i in pop 0 0 , the chromosomes in pop 0 1 are traversed and checked if there exists a neuron r 1;j that can be swapped with the neuron of chrom 0;i denoted by r 0;j . Neurons are eligible for swapping given that r 0;j and r 1;j are not disconnected or output neurons in chrom 0;i and chrom 1;i respectively. If the condition is met, parent chromosomes chrom 0;i and chrom 1;i are appended to P, and neuron index jis appended to S. This process is done until there are dnumber of parent chromosome pairs. After the crossover process, the generated offspring chromosomes are appended to the current population which is now be the passed as the next population. Otherwise, in the event that the size of the current population is equal to the set population size, the process of generating Pand S are done until the set maximum population size is met. However, not all chromosomes of pop 0 0 and pop 0 1 undergo the crossover process. The chance of being chosen to be parent chromosomes is determined by the set crossover rate. The chromosomes chosen to be parent chromosomes undergo the crossover process and its corresponding offspring chromosomes are appended to the next population. Those that aren’t chosen are added to the next population as is. The process of crossing over is executed in the GPU. There are three lists copied from the host to the device: two lists of parents P 0 ¼fchrom 0;0 ;...;chrom 0;p gand P 1 ¼fchrom 1;0 ;...;chrom 1;p gwith pbeing the number of parent chromosome pairs, and S¼ r 0 ;...; r p which contains the indices of the neurons to be swapped per parent chromosome pair of the same index. Crossing over is done by swapping the values of the out-going synapses of the pair’s chosen neuron r . 4. Experimentation This section discusses the technicalities of experimentation. We have adapted two variants of the GA function design from [18] in our experiments which are further discussed in subSection 4.2. Input used for the experiments are SN P Systems that simulate binary addition and binary subtraction; these are divided into different categories each with its own P init , namely Baseline, Original, and Adversarial. This is in order to measure evolution variant effectiveness. These categories are discussed in subSection 4.1. Finally, each experiment run has parameters that are set beforehand. These parameters include the population size, maximum number of inner loops, maximum number of generations per inner loop, mutation rate, crossover rate, and selection rate. These are enumerated in subSection 4.3. Source code files used in experimentation and results can be accessed inhttps://github.com/katkmh/ESNP-GPU. 4.1. Input Design Due to the design of the simulator, the framework only accepts deterministic SN P systems. The initial SN P systems P init used are limited to binary addition and binary subtraction. However, these operations are represented by three categories: Baseline,Original, and Adversarial. The SN P systems under the Baseline category have the closest output approximation while having the smallest designs in terms of number of neurons, synapses, and rules. The SN P systems under the Original category are significantly larger and less precise than those under the baseline category. There is one system under this category for binary addition, in which its set of neurons is a superset of its baseline counterpart. There are two systems for binary subtraction, referred to as Original 1 and Original 2 in this study. The decision for having two such systems are discussed in [18] in which Original 2 contains all the neurons from the baseline while Original 1 does not. The SN P systems under the Adversarial category follow the designs of its baseline counterparts but with madditional neurons, where mis the degree of its corresponding baseline counterpart. The process of generating adversarial versions from its baseline counterparts is further discussed in Section 3.3 of [18]. See Cfor visualizations of the input SN P systems. 4.2. Genetic AlgorithmVariants There are four variants of the genetic algorithm design that differs on two variables: whether the formation of self-loops is allowed during evolution and whether inactive or disconnected neurons are deleted. From the work of [18], the variants that are set to delete inactive neurons yielded higher average precision values compared to the variants that are set to delete disconnected neurons. Due to the aforementioned results from the previous study, the objective of this study, and our current technical limitations, we have only chosen two among the available variants: variant 1 and variant 3 in [18] which will be referenced as GA Variant 1 and GA Variant 2 in this paper respectively. In these two variants, the following validations remain consistent: 1. In chromosome comparison, only active neurons are counted. 2. After the execution of the mutation or crossover process, the resulting chromosome must always have a path from at least one input neuron to the environment neuron. If it does not meet this condition, the chromosome pre-mutation or pre-crossover is kept. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 146 3. If the resulting chromosome from the mutation or crossover process is the same as the input to the current outer loop, it is invalid. The chromosome pre-mutation or pre-crossover is kept. 4. The final SN P System P final must never have a self-loop. A method is called before declaring it as final. The method removes self-loops by transforming the neurons as shown in Fig. 5. 5. Inactive neurons are deleted after the end of each outer loop. The two variants only differ in allowing the formation of selfloops during evolution. Self-loops are neurons in which there exists an outgoing synapse connected to itself. GA Variant 1 does not allow self-loops during evolution while GA Variant 2 does. 4.3. Experiment Setup 4.3.1. GA Framework For each category under each operation sans Original 1 of binary subtraction, the GA framework is ran 5 times with the corresponding initial population set. Original 1 of binary subtraction was run 10 times; this decision is discussed in Section 5. There are two versions of the GA framework that is executed during the experimentation: the first one is from [18] which runs wholly sequentially and the second one uses GA function in Fig. 4 that runs partially parallel; these will be identified in this paper as GA-CPU and GA-GPU respectively. This is to provide comparison on the improvements of the SN P System evolution and simulation. Both GA framework versions have the 2 variants in SubSection 4.2 which will be identified in this paper as CPU Variant 1 and CPU Variant 2 for GA-CPU, and GPU Variant 1 and GPU Variant 2 for the GA-GPU. Given 7 P init designs used on 2 GA frameworks with 2 variants, with 6 P init ’s executed for 5 times and 1 P init executed for 10 times, the total number of runs during experimentation is 160. Two sets of data are prepared prior to executing the GA Frameworks: 4.3.2. Set of fitness cases There are two sets of fitness cases Sgenerated, one for each operation (binary addition and binary subtraction). Each set of fitness cases contain 100 pairs of input spike trains and their corresponding output spike train according to the operation. The maximum length of the input spike trains is 50. Table 1. 4.3.3. Set of initial populations With multiple experiment runs, a minor modification to the genetic algorithm function was made in order to maintain a uniform starting ground across all evolutions under each P init . The first step of the genetic algorithm function described in subSection 3.3, where a population 0 is generated from chrom init by mutating it ptimes, was run prior to evolution. By doing so, it generates an initial population population init for each of the 7 initial SN P systems P init . With the first step of the GA Function processed, all experiment runs take in population init as the input instead of chrom init before proceeding to the simulation process. By maintaining the same initial population across all experiment runs of each input SN P system P init , bias in the results due to one run having a more advantageous initial population can be avoided. Table 2. The values of por sizes of each generated initial population vary according to the operation and its category as shown in Table 3. Aside from the population sizes, the framework utilizes the same default parameters in every run: 1. Maximum number of inner loops: 10 2. Maximum number of generations per inner loop: 75 3. Mutation rate: 0.5 4. Crossover rate: 0.3 5. Selection rate: 0.6 5. Results and Discussion All experiments were run using Computing and Archiving Research Environment (COARE)’s High-Performance Computing (HPC) service [29] with the specifications cited in Table 4. Two types of experiments were conducted. First, we averaged the runtimes of the parallel stage (mutation, simulation and crossover, as shown in Fig. 4). The aim is to analyze how effective is the parallelization of the bottlenecks in the evolution process. Moreover, we use the average to take into account the variance of runs with operations such as binary substraction. Secondly, we performed overall executions, which include parts of the evolutionary algorithm and data transfer, to check how effective can be this preliminary version of the framework. 5.1. Average runtimes As mentioned in Section 4.3, the same parameters and corresponding initial population was used for all the experiments so there is expected similar results with respect to precision, number of generations in evolution, number of outer and inner loops, and resulting SN P system topology. The focus of this discussion will revolve around the runtime of the GA function, where an initial SN P system is evolved to a system with a more desirable precision and topology. Figs. 6a and 6b show the average evolution runtime (only parallel stages) of the CPU and GPU Variants using P init for binary addition and binary subtraction respectively. Average Evolution Runtime of each GA framework variant is taken by computing the average amount of time that the GA function took to execute per inner loop across all corresponding runs. The data visualized in the figures is recorded in Tables 5 and 6. These tabulated data show the average evolution runtime, standard deviation (Std Dev), and coefficient of variance (CV). The averages are taken using the evolution runtime values of all experiment runs made by each variant of each GA framework. GA-GPU yielded faster evolution runtimes in comparison to its GA-CPU counterparts of up to 9.14 speedup as shown in the figures, supplemented by the speedup data recorded in Table 7. These results show the optimization, in terms of runtime, made on the evolution of SN P systems when executing the GA operators on the CPU and GPU. Variant 2 precision values are higher than Variant 1 in both GACPU and GA-GPU which is similar to the results of [18]. This supports the hypothesis mentioned wherein allowing self-loops dur- Fig. 5. Removing a self-loop. Taken from [18], this figure illustrates how a chromosome with a self-loop is transformed. The system for computing a Boolean function of three variables in [8] utilizes a similar mechanism in the form of the auxiliary neurons. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 147 ing evolution may help in escaping a local maximum of the SN P system. On the other hand, Variant 2 performed better than Variant 1 in both GA-CPU and GA-GPU except on Original 1 P init and Original 2 P init for binary subtraction in terms of evolution runtime. Due to the stochastic nature of genetic algorithms, the SN P systems during evolution began to have varying topology, i.e., number of neurons and synapses. It is speculated that GPU Variant 2 performed slower than GPU Variant 1 for binary subtraction because the former had more runs where the SN P system contained more neurons and synapses. Thus, slowing down the average evolution runtime. The input data with greatly varied results from the majority is Original 1 P init for binary subtraction. Originally, there were only 5 runs made for this input data, however with large standard deviation values as shown in Table 8 consistent in GA-CPU and GAGPU, we have decided to execute 5 additional runs for each variant. With the updated values considering the additional runs (Table 6, row of Original 1), there is an inferred unpredictability with the runtime of evolution of initial chromosomes with low precision. There is also an inferred unpredictability in the final chromosome’s precision and topology as shown in B. While the precision values of the final chromosomes generated by GA-CPU and GA-GPU played between 58%to 65%, GPU Variant 2 made a notable run that was able to produce a final chromosome with precision value 93%shown in Fig. D.19. The additional factors that may contribute to the unpredictability of the Original category’s evolution’s performance cannot be inferred with the results at hand. However, it is suspected that this is due to the observation that the GA framework is most effective in evolving systems with a subgraph that already approximates function f, defined in subSection 3.2. As mentioned in Section 4.1, the Adversarial category was built from the Baseline category but with additional neurons thus there exists a subgraph that approximates function f. Given that the Baseline category’s precision is already 100%, if the evolution process is successful, the P init under Adversarial category could be evolved into a P final with 100%precision and topology similar to or exactly like the P init under Baseline category it was built from. This observation is supported by results shown in A. On the other hand, Original category was not generated from the Baseline category. It may have sets of neurons that are supersets of its baseline counterparts’, but this does not guarantee a subgraph that approximates function f. This observation was also present in [18] wherein the mentioned theory was that during evolution, it was getting trapped in its local maximum. In the event that precision could no longer be improved, the evolution exiting criteria that it will be able to satisfy is reduced to either decreased topology or maximum number of generations. It is shown in Athat the resulting topology of Binary Subtraction’s Original 2 is mostly smaller in comparison to Original 1 despite Original 2 having a set of neurons that is a superset of the Baseline category’s set of neurons. 5.2. Overall runtimes Table 9 shows the overall time of a single run for all categories and variants, and compares the GPU vs CPU performance. The table focuses on not only the parallel part of the algorithm, but on the overall process (see Fig. 4). This includes data transfer from and to the GPU, what becomes a bottleneck. GA-GPU is slower than GA-CPU in overall process for variant 1. However, for variant 2, Table 1 A sample run of the SN P system presented in Fig. 1. This system simulates the binary addition function: arepresents 1 and krepresents 0. The input spike trains in1and in2 correspond to what is received by the input neurons r1and r2, respectively. In this example, in1¼111 and in2¼101, thus the output spike train must be out ¼1100. time in 1 in 2 r 1 r 2 r 3 r 4 r 5 out 0aaa aka----- k 1aa akaa--- kk 2aaaka 2 a 2 a 2 kkk 3--aa a 2 a 2 a 2 k kkk 4---- a 3 a 3 a 3 kk kkk 5----aaa akk kkk 6 - - -----aakk kkk Table 2 Characteristics of the initial SN P Systems used in the experiments. Operation Category Precision Number of Neurons Number of Synapses Binary Addition Baseline 100% 5 10 Original 75.138375% 11 14 Adversarial 31.05496% 8 19 Binary Subtraction Baseline 100% 8 15 Original 1 38.60768% 32 38 Original 2 26.04025% 33 41 Adversarial 60.45469% 16 59 Table 3 Values of population sizes for each operation and category used in the experiments. Operation Category Population Size Binary Addition Baseline 15 Original 40 Adversarial 40 Binary Subtraction Baseline 20 Original 1 80 Original 2 80 Adversarial 40 Table 4 Specifications of HPC system where experiments were run. System Specifications Values Operating System CentOS Linux 7 (Core) Processor Intel(R) Xeon(R) CPU E5-2699 v4 @ 2.2 GHz Cores 240 (allocated per user) RAM 24 GB GPU NVIDIA Tesla P40 @ 1.3 GHz CUDA Cores/GPU 30 (multiprocessors) 12 (cores/processor) = 360 CUDA version 10.1 Device memory size 24 GB R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 148 Fig. C.13. SN P systems that simulates binary subtraction under Original 2 category. Fig. C.12. SN P systems that simulates binary subtraction under Original 1 category. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 155 Fig. C.14. SN P systems that simulates binary subtraction under Adversarial category. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 156 Appendix D. Noteworthy Final SN P Systems from experiments Figs. D.15,D.16,D.17, and D.18 challenge the theory that evolving an input SN P System under Adversarial category will yield the system under Baseline category from which it was built from. GPU Variant 1 was able to design a system that differed from its corresponding Baseline category SN P System in terms of topology while still yielding 100%precision when simulated on the same fitness cases. Fig. D.19 shows the outlying final SN P System from GPU Variant 2 in Fig. 7b that has 93%precision. Fig. D.15. P final from a GPU Variant 1 run on input SN P system simulating binary addition under Adversarial category (Fig. C.10). This system has 100%precision but differed in topology from the baseline system it was built from, now with 7 neurons and 14 synapses. The framework was able to generate this system in 228 generations. Fig. D.16. P final from a GPU Variant 1 run on input SN P system simulating binary addition under Adversarial category (Fig. C.10). This system has 100%precision but differed in topology from the baseline system it was built from, now with 7 neurons and 20 synapses. The framework was able to generate this system in 161 generations. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 157 Fig. D.17. P final from a GPU Variant 1 run on input SN P system simulating binary subtraction under Adversarial category (Fig. C.14). This system has 100%precision but differed in topology from the baseline system it was built from, now with 9 neurons and 17 synapses. The framework was able to generate this system in 719 generations. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 158 Fig. D.19. P final from a GPU Variant 2 run on input SN P system simulating binary subtraction under Original 1 category (Fig. C.12). It is mentioned in Section 5that category Original 1 had the most varied results. As shown in the box plot in Fig. 7b, there is an extreme outlier wherein the resulting system’s precision is 93%. The system illustrated above is the outlying resulting system. Fig. D.18. P final from a GPU Variant 1 run on input SN P system simulating binary subtraction under Adversarial category (Fig. C.14). This system has 100%precision but differed in topology from the baseline system it was built from, now with 9 neurons and 18 synapses. The framework was able to generate this system in 624 generations. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 159 References [1] G. Pa ˘un, From Cells to (Silicon) Computers, and Back, Springer New York, New York, NY, 2008, pp. 343–371, https://doi.org/10.1007/978-0-387-68546-5_15, URL:https://doi.org/10.1007/978-0-387-68546-5_15. [2] G. Pa ˘un, Computing with membranes, J. Comput. Syst. Sci. 61 (1) (2000) 108– 143. [3] G. Zhang, M. Gheorghe, L. Pan, M.J. Perez-Jimenez, Evolutionary membrane computing: a comprehensive survey and new results, Inf. Sci. 279 (2014) 528– 551. [4] M. Gheorghe, G. Pa ˘un, M. Pérez-Jiménez, G. Rozenberg, Frontiers of membrane computing: Open problems and research topics, Intern. J. Found. Computer Sci. (2013) 171–249. [5] X. Huang, G. Zhang, H. Rong, F. Ipate, Evolutionary design of a simple membrane system, in: M. Gheorghe, G. Pa ˘un, G. Rozenberg, A. Salomaa, S. Verlan (Eds.), CMC 2011: Membrane Computing, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 203–214. [6] Z. Ou, G. Zhang, T. Wang, X. Huang, Automatic design of cell-like p systems through tuning membrane structures, initial objects and evolution rules, Int. J. Unconv. Comput. 9 (5–6) (2013) 425–443. [7] G. Zhang, H. Rong, Z. Ou, M. Pérez-Jiménez, M. Gheorghe, Automatic design of deterministic and non-halting membrane systems by tuning syntactical ingredients, IEEE transactions on nanobioscience 13. doi:10.1109/ TNB.2014.2341618. [8] G. Pa ˘un, Spiking neural p systems. a tutorial, Bulletin of the European Association for Theoretical Computer Science EATCS. [9] J.R. Cheng, M. Gen, Parallel genetic algorithms with gpu computing, in: Industry 4.0-Impact on Intelligent Logistics and Manufacturing, IntechOpen, London, United Kingdom, 2020, https://doi.org/10.5772/intechopen.89152. [10] F.G.C. Cabarle, H. Adorna, M.Á. Martínez-del-Amor, A spiking neural p system simulator based on cuda, in: M. Gheorghe, G. Pa ˘un, G. Rozenberg, A. Salomaa, S. Verlan (Eds.), CMC 2011: Membrane Computing, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 87–103. [11] F.G.C. Cabarle, H. Adorna, M. Á. Martínez-del-Amor, An improved gpu simulator for spiking neural p systems, in: 2011 Sixth International Conference on Bio-Inspired Computing: Theories and Applications, IEEE, IEEE Computer Society, Washington, DC, United States, 2011, pp. 262–267. [12] F. Cabarle, H. Adorna, M.Á. Martínez-del-Amor, Simulating spiking neural p systems without delays using gpus, Int. J. Natural Computing Res. (IJNCR) 2 (2) (2011) 19–31. [13] F.G.C. Cabarle, H.N. Adorna, M.Á. Martínez-del-Amor, M.J. Pérez-Jiménez, Improving gpu simulations of spiking neural p systems, Romanian, J. Inform. Sci. Technol. 15 (2012) 5–20. [14] J. Carandang, J.M.B. Villaflores, F.G.C. Cabarle, H.N. Adorna, M.Á. Martínez-del- Amor, Cusnp: Spiking neural p systems simulators in cuda, Romanian, J. Inform. Sci. Technol. 20 (1) (2017) 57–70. [15] B.C.D. Aboy, E.J.A. Bariring, J.P. Carandang, F.G.C. Cabarle, R.T. De La Cruz, H.N. Adorna, M.Á. Martínez-del-Amor, Optimizations in cusnp simulator for spiking neural p systems on cuda gpus, in: 2019 International Conference on High Performance Computing Simulation (HPCS), IEEE, 2019, pp. 535–542, https:// doi.org/10.1109/HPCS48598.2019.9188174. [16] J. Sanders, E. Kandrot, CUDA by Example: An Introduction to General-Purpose GPU Programming, Pearson Education Inc, Boston, MA, United States, 2011. [17] A. Klöckner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov, A. Fasih, PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation, Parallel Comput. 38 (3) (2012) 157–174, https://doi.org/10.1016/ j.parco.2011.09.001. [18] L.J. Casauay, F.G.C. Cabarle, I.C.H. Macababayao, R.T.A. de la Cruz, H.N. Adorna, X. Zeng, M.Á. Martínez-del-Amor, A framework for evolving spiking neural p systems, Int. J. Unconventional Computing 16 (2–3) (2021) 121–139. [19] C.C.R. Zarate, F.G.C. Cabarle, I.C. Macababayao, R.T. De la Cruz, Evolving spiking neural p systems by fixing neurons, and varying rules and synapses, Philippine Computing Journal. (Special Issue on P systems) 14 (2) (2020) 21–30. [20] J.G.E. Juico, J.L. Silapan, F.G.C. Cabarle, I.C. Macababayao, R.T. De la Cruz, Evolving spiking neural p systems with polarization, Philippine Computing Journal. (Special Issue on P systems) 14 (2) (2020) 11–20. [21] M. Ionescu, G. Pa ˘un, T. Yokomori, Spiking neural P systems, Fundamenta Informaticae 71 (2, 3) (2006) 279–308. [22] K. Fatahalian, J. Sugerman, P. Hanrahan, Understanding the efficiency of gpu algorithms for matrix-matrix multiplication, in: Proceedings of the ACM SIGGRAPH/EUROGRAPHICS conference on Graphics hardware, Association for Computing Machinery, New York, United States, 2004, pp. 133–137. [23] X. Zeng, H. Adorna, M.Á. Martínez-del-Amor, L. Pan, M.J. Pérez-Jiménez, Matrix representation of spiking neural p systems, in: M. Gheorghe, T. Hinze, G. Pa ˘un, G. Rozenberg, A. Salomaa (Eds.), Membrane Computing, Springer Berlin Heidelberg, Berlin, Heidelberg, 2011, pp. 377–391. [24] J. Holland, Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control and Artificial Intelligence, MIT Press, Cambridge, MA, USA, 1992. [25] R.L. Haupt, S.E. Haupt, The Binary Genetic Algorithm, Ch. 2, John Wiley & Sons Ltd, Hoboken, New Jersey, 2003, pp. 27–50, https://doi.org/10.1002/ 0471671746.ch2, arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/ 0471671746.ch2, URL:https://onlinelibrary.wiley.com/doi/abs/10.1002/ 0471671746.ch2. [26] D. Luebke, Data parallel computing, in: D.B. Kirk, W.H. Wen-Mei (Eds.), Programming Massively Parallel Processors (Third Edition), Ch. 2, Morgan Kaufmann, Cambridge, MA, USA, 2017, pp. 19–41, https://doi.org/10.1016/ B978-0-12-811986-0.00002-9, URL:http://www.sciencedirect.com/science/ article/pii/B9780128119860000029. [27] L.F. Macías-Ramos, I. Pérez-Hurtado, M. García-Quismondo, L. Valencia- Cabrera, M.J. Pérez-Jiménez, A. Riscos-Núñez, A p–lingua based simulator for spiking neural p systems, in: M. Gheorghe, G. Pa ˘un, G. Rozenberg, A. Salomaa, S. Verlan (Eds.), Membrane Computing, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012, pp. 257–281. [28] D. Gusfield, Algorithms on stings, trees, and sequences: Computer science and computational biology, Acm Sigact News 28 (4) (1997) 41–60. [29] Coare website, last accessed December 2020. URL:https://asti.dost.gov. ph/coare/wiki/Main/. [30] J. Dong, M. Stachowicz, G. Zhang, M. Cavaliere, H. Rong, P. Paul, Automatic design of spiking neural p systems based on genetic algorithms, Int. J. Unconv. Comput. 16 (2–3) (2021) 201–216. [31] K.O. Stanley, R. Miikkulainen, Evolving neural networks through augmenting topologies, Evolutionary computation 10 (2) (2002) 99–127. [32] L.L. Custode, H. Mo, G. Iacca, Neuroevolution of spiking neural p systems, in: International Conference on the Applications of Evolutionary Computation (Part of EvoStar), Springer, 2022, pp. 435–451. [33] L.L.Custode,H.Mo,A.Ferigo,G.Iacca,Evolutionaryoptimizationofspikingneuralp systems for remaining useful life prediction, Algorithms 15 (3) (2022) 98. [34] T. Song, L. Pan, T. Wu, P. Zheng, M.D. Wong, A. Rodríguez-Patón, Spiking neural p systems with learning functions, IEEE Trans. Nanobioscience 18 (2) (2019) 176–190. [35] T. Wu, L. Pan, Q. Yu, K.C. Tan, Numerical spiking neural p systems, IEEE Trans. Neural Networks Learn. Syst. 32 (6) (2020) 2443–2457. [36] M.A. Martínez-del-Amor, D. Orellana-Martín, I. Pérez-Hurtado, F.G.C. Cabarle, H.N. Adorna, Simulation of spiking neural p systems with sparse matrix-vector operations, Processes 9 (4). doi:10.3390/pr9040690. URL:https://www. mdpi.com/2227-9717/9/4/690. [37] M.Á. Martínez-del-Amor, I. Pérez-Hurtado, D. Orellana-Martín, M.J. Pérez- Jiménez, Adaptative parallel simulators for bioinspired computing models, Future Generation Computer Systems 107 (2020) 469–484, https://doi.org/ 10.1016/j.future.2020.02.012, URL:http://www.sciencedirect.com/science/ article/pii/S0167739X19308817. [38] R. Ceterchi, A.I. Tomescu, Spiking neural p systems–a natural model for sorting networks, in: Proceedings of the Sixth Brainstorming Week on Membrane Computing, 4–8 February 2008, Fénix Editora, Sevilla, Spain, 2008, pp. 93–105. [39] D. Abrahams, A. Gurtovoy, C++ Template Metaprogramming: Concepts, Tools, and Techniques from Boost and Beyond (C++ in Depth Series), Addison-Wesley Professional, Boston, Massachusetts, United States, 2004. [40] L. Valencia-Cabrera, I. Pérez-Hurtado, M.Á. Martínez-del-Amor, Simulation challenges in membrane computing, J. Membrane Computing 2 (2020) 1–11, https://doi.org/10.1007/s41965-020-00056-w. [41] I. Pérez-Hurtado, D. Orellana-Martín, M.A. Martínez-del-Amor, L. Valencia- Cabrera, A. Riscos-Núñez, A new p-lingua toolkit for agile development in membrane computing, Inf. Sci. 587 (2022) 1–22, https://doi.org/10.1016/j. ins.2021.12.003. Rogelio V. Gungon received his B.Sc. in computer science from the University of the Philippines Diliman, in Quezon City, Philippines (working with the Algorithms and Complexity laboratory of the Department of Computer Science). His research interests include GPU computing and spiking neural P systems. Katreen Kyle M. Hernandez received her B.Sc. in computer science from the University of the Philippines Diliman, in Quezon City, Philippines (working with the Algorithms and Complexity laboratory of the Department of Computer Science). Her research interests include GPU computing and spiking neural P systems. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 160 Francis George C. Cabarle received his Ph.D. in Computer Science in 2015 from the University of the Philippines Diliman, where he is currently an associate professor. He was a visiting researcher with the Research Group of Natural Computing from 2014 to 2015, then headed by Prof. Mario de J. Pérez Jiménez, at the University of Seville, in Spain. From 2016 to 2019, FGCC performed postdoctoral work with Prof. Xiangxiang Zeng et al. at Xiamen University, in China. His research interests include the intersections of bioinspired computing, computability, and parallel computing. Ren Tristan A. de la Cruz is a Ph.D. student in computer science from the Algorithms and Complexity Laboratory, Department of Computer Science, University of the Philippines Diliman, Quezon City, Philippines. He received his M.Sc. degree in computer science from the same university. His current research interests include membrane computing and parallel computing. Henry N. Adorna is a Professor of (Theoretical) Computer Science since 2007 and Scientist at UP Diliman. He spent two years (from 2000 to 2002) at Lehrstuhl fuer Informatics I of RWTH Aachen under the DAAD PhD Sandwich Program to do his dissertation on abstract communication complexity of uniform computing models under Prof. Dr. Juraj Hromkovic. HNA went home to obtain his Ph.D. in mathematics from the UP Diliman in 2002. HNA works primarily on P systems since 2009. Other research interests of HNA include algorithmics for NP-hard problems, particularly combinatorial hard problems in computational biology, discrete mathematics, natural and unconventional computing models. Miguel Á. Martínez-del-Amor obtained his M.Sc. in Computer Engineering at the University of Murcia (2008), and his Ph.D. in Computer Science and Artificial Intelligence at the University of Seville (2013). He was a research associate at Fraunhofer IIS (Germany) from 2014 to 2017. He is an Assistant Professor at the University of Seville since 2017, and an NVIDIA DLI ambassador since 2018. His research interests are mainly focused on GPU computing and its applications to Artificial Intelligence (Machine Learning), Natural Computing (Membrane Computing) and Image Coding (JPEG standards). David Orellana-Martín received his B.Sc. in 2014, his M. Sc. in 2016 and his PhD in 2019 at the University of Seville. He received the best PhD thesis of 2019 Award from the International Membrane Computing Society. He was an ERCIM fellow at NTNU (Norway) from 2019 to 2021, and nowadays he is a postdoctoral fellow at the University of Seville. His main research interests are Computational Complexity Theory, Unconventional Computing (specially Membrane Computing), as well as Machine Learning and High Performance Computing. He is member of the Website Committee of the International Membrane Computing Society. Ignacio Pérez-Hurtado received his B.Sc. in 2003, and afterwards he worked for an IT company. From 2007 to 2013 he worked at the RGNC (University of Seville) receiving his Ph.D. in 2010. From 2014 to 2017 he was a postdoctoral fellow at Pablo Olavide University (Spain), researching in Social Robotics. He is currently an Assistant Professor at the University of Seville. His main research interests are related to P system simulators, programming languages and algorithms for social navigation in robotics. R.V. Gungon, Katreen Kyle M. Hernandez, Francis George C. Cabarle et al. Neurocomputing 503 (2022) 140–161 161