FCMpy: a python module for constructing and analyzing fuzzy cognitive maps

FCMpy is an open-source Python module for building and analyzing Fuzzy Cognitive Maps (FCMs). The module provides tools for end-to-end projects involving FCMs. It is able to derive fuzzy causal weights from qualitative data or simulating the system behavior. Additionally, it includes machine learning algorithms (e.g., Nonlinear Hebbian Learning, Active Hebbian Learning, Genetic Algorithms, and Deterministic Learning) to adjust the FCM causal weight matrix and to solve classification problems. Finally, users can easily implement scenario analysis by simulating hypothetical interventions (i.e., analyzing what-if scenarios). FCMpy is the first open-source module that contains all the functionalities necessary for FCM oriented projects. This work aims to enable researchers from different areas, such as psychology, cognitive science, or engineering, to easily and efficiently develop and test their FCM models without the need for extensive programming knowledge.

An FCM represents a system as a directed signed graph where components are represented as nodes and the causal relationships between these components are represented by weighted directed edges. The dynamics of the system are examined by simulating its behavior over discrete simulation steps. In general, FCMs can be constructed based on the inputs of domain experts (i.e., expert-based FCMs), data collected about the system (e.g., data-driven approaches) or the combination of the two (i.e., hybrid approaches) (Mkhitaryan et al., 2020). fuzzy membership functions based on the expert ratings, (3) combine the membership functions resulting from the second step with an aggregation operation, and (4) defuzzify the aggregated membership functions (Mkhitaryan et al., 2020;Mago et al., 2012Mago et al., , 2013. In this section, we first describe methods for reading data from different file formats and then describe the methods for constructing expert-based FCMs based on qualitative data following the four steps described above.

Data handling
The available open source solutions for expert-based FCMs do not provide utilities for working with different file types thus limiting their usability. Data on FCMs include the edges (represented as pairs of source/target) and the associated linguistic ratings of the survey participants. The ExpertFcm class provides a read_data method for reading data from .csv, .xlsx, and .json files (see the code snippet below).
The corresponding files should satisfy certain requirements that are described in detail in the PyPI documentation. Before using the read_data method we first need to define the linguistic terms (explained in detail in the next section) used in the data. We can do that by using the linguistic_terms method. The read_data requires the file path as an argument. The additional arguments that depend on the file extension (e.g., csv, json, xlsx) should be specified as keyword arguments. For the .xlsx and .json files, when the optional check_consistency argument is set to True then the algorithm checks whether the experts rated the causal impact of the edges (source-target pairs) consistently in terms of the valence of the causal impact (positive or negative causality). If such inconsistencies are identified, the method outputs a separate .xlsx file that documents such inconsistencies. The read_data method returns an ordered dictionary where the keys are the experts' IDs (or the names of the excel sheets in the case of an excel file or the row index in case of a csv file) and the values are pandas dataframes with the expert inputs.
It is often useful to check the extent to which the participants agree on their opinions with respect to the causal relationships between the edges. This is often done by calculating the information entropy (Kosko, 1996) expressed as: where p i is the proportion of the answers (per linguistic term) about the causal relationship. The value of entropy is always greater than or equal to zero. If all the experts give the same answer about a particular edge (e.g., from C1 to C3), the entropy score for that connection will be 0. However, if experts disagree about the importance of a connection (e.g., in the edge between C1 and C2), the entropy value would increase. The entropy scores can be calculated with the entropy method (see the code snippet below).

FOUR STEPS FOR OBTAINING CAUSAL WEIGHTS
To convert the qualitative ratings of the domain experts to numerical weights via fuzzy logic, we must (1) define the fuzzy membership functions, (2) apply a fuzzy implication rule, (3) combine the membership functions, and (4) defuzzify the aggregated membership functions to derive the numerical causal weights. Table 2 gives a concise listing of the methods for building expert based FCMs.
Step 1: define fuzzy membership functions Fuzzy membership functions are used to map the linguistic terms to a specified numerical interval (i.e., universe of discourse). In FCMs, the universe of discourse is specified in the [-1, 1] interval where the negative causality is possible or in the [0, 1] interval otherwise. The universe of discourse can be specified with the universe setter (see the code snippet below).
>>> import numpy as np >>> fcm.universe = np.arange(-1, 1.001,.001) To generate the fuzzy membership functions we need to decide on the geometric shape that would best represent the linguistic terms. In many applications, a triangular membership function is used (Zadeh, 1971). The triangular membership function specifies the lower and the upper bounds of the triangle (i.e., where the meaning of the given linguistic term is represented the least) and the center of the triangle (i.e., where the meaning of the given linguistic term is fully expressed).
The linguistic_terms method sets the linguistic terms and the associated parameters for the triangular membership function (see the code snippet below). The keys in the above dictionary represent the linguistic terms and the values are lists that contain the parameters for the triangular membership function (i.e., the lower bound, the center and the upper bound) (see Fig. 1). After specifying the universe of discourse and the linguistic terms with their respective parameters one can use use the automf method to generate the membership functions (see the code snippet below).
>>> fcm.fuzzy_membership = fcm.automf(method = 'trimf') In addition to the triangular membership functions, the automf method also implements gaussian membership functions ('gaussmf') and trapezoidal membership functions ('trapmf') (based on sci-kit fuzzy module in python). Step 2: apply the fuzzy implication rule To determine the level of endorsement of the linguistic terms for a given pair of concepts, one must first identify the level of endorsement of the given terms by the participants. This is done by calculating the proportion of the answers to each linguistic term for a given edge. Consider a case where 50% of the participants (e.g., domain experts) rated the causal impact of an antecedent on the consequent as Positive High, 33% rated it as Positive Very High and the 16% rated it as Positive Medium. Subsequently, a fuzzy implication rule is used to activate the corresponding membership functions. Two such rules are often used, namely Mamdani's minimum and Larsen's product implication rule (Nandi, 2012). The Mamdani minimum fuzzy implication rule is expressed as: where l A ðxÞ and l B ðyÞ denote the membership value x to the linguistic term A and the membership value y to the linguistic term B respectively. The Mamdani rule cuts the membership function at the level of endorsement (see Fig. 2A). In contrast, Larsen's implication rule re-scales the membership function based on the level of endorsement (see Fig. 2B) and is expressed as: (3)  We can use fuzzy_implication method to apply the selected implication method (see the available methods in Table 3  Step 3: aggregate fuzzy membership functions In the third step, we must aggregate the activated membership functions taken from the previous step. This is commonly done by applying the family maximum aggregation operation. Alternative methods for aggregating membership functions include the family Algebraic Sum (see Eq. (4)), the family Einstein Sum (see Eq. (5)) and the family Hamacher Sum (see Eq. (6)) (Piegat, 2001).
One can use the aggregate method to aggregate the activated membership functions (see the available aggregation method in Table 4  where x and y are the membership values of the linguistic terms involved in the problem domain after the application of the implication rule presented in the previous step. Step 4: defuzzify the aggregated membership functions The last step includes the calculation of the crisp value based on the aggregated membership functions (a.k.a. defuzzification). Among the available defuzzification methods (see the available defuzzification methods in Table 5) the most commonly used method is the centroid method (a.k.a. center of gravity) (Stach et al., 2005). We can apply the dedicated defuzz method to derive the crisp value (see Fig. 3 and the code snippet below).
>>> dfuz = fcm.defuzz(x=fcm.universe, mfx=aggregated, method='centroid') The above mentioned four steps can either be done and controlled independently as we have shown, or users can rely on a single build method that implements those steps to calculate the numerical weights for all the concept pairs in the data (see the code snippet below). The method returns a pandas dataframe with the calculated weights.

Argument
Option Description Max of maximum In this section, we showed how the structure of the FCMs can be derived based on data collected from the domain experts. In the subsequent section, we illustrate how the system behavior can be simulated on top of the defined FCM structure.

SIMULATING THE SYSTEM BEHAVIOR WITH FCMS
The dynamics of the specified FCM are examined by simulating its behavior over discrete simulation steps. In each simulation step, the concept values are updated according to a defined inference method (Papageorgiou, 2011b). The FcmSimulator module implements the following three types of inference methods (see the available options in Table 6): Kosko: Modified Kosko: Rescaled: where a ðtÞ j is the value of concept j at the simulation step t and w j;i is the causal impact of concept j on concept i. Note that a (transfer) function f(x) is applied to the result. As shown in the equations above, this function is necessary to keep values within a certain range (e.g., Table 6 Defuzzification methods.

Argument
Option Description [0,1] for sigmoid function or [-1,1] for hyperbolic tangent). In the current version, four such functions are implemented (see the available options in Table 6): Sigmoid: Hyperbolic tangent: Trivalent: ; x 2 R; binds node values to fÀ1; 0; 1g (13) where, x is the value calculated by applying the above mentioned inference methods and the is a steepness parameter for the sigmoid function. The simulation is run until either of two conditions is met: (1) some concepts of interest have a difference lower than a given threshold between two consecutive steps (default value 0.001), or (2) a user-defined maximum number of iterations is reached. If we denote by S the activation vector for a subset of concepts of interest (i.e., the outputs of the FCMs), then the first condition can be stated as: The simulate method takes the initial state vector and the FCM weight matrix (a.k.a., connection matrix) and applies one of the mentioned update functions over a number of simulation steps (see the simulation results using different inference and transfer methods in Figs. 4-6). One can specify the output concepts by supplying a list of these concepts to the respective output_concepts argument. If the output_concepts argument is not specified then all the concepts in the FCM are treated as output concepts and the simulation stops when all the concepts change by less than the threshold between two consecutive steps.

LEARNING ALGORITHMS FOR FCMS
As shown in the previous sections, FCMs are often constructed based on experts' knowledge about the system. In certain domains of applications, modelers either optimize the FCMs constructed by the experts and/or constructing FCMs entirely based on the data collected about the systems. A set of machine learning algorithms were developed to meet these tasks and have previously been applied to numerous fields such as the optimization of industrial processes (Papageorgiou, Stylios & Groumpos, 2006;Stach, Kurgan & Pedrycz, 2008;Papageorgiou, 2011a), decision making (Poczeta, Papageorgiou & Gerogiannis, 2020), and classification (Nápoles et al., 2014;Nápoles, Jastrzębska & Salgueiro, 2021). In the proposed library, we include three types of algorithms used for edge optimization, FCM generation, and classification. We used the state-of-the-art methods (Nápoles, Jastrzębska & Salgueiro, 2021;Nápoles et al., 2020) and the foundational ones that have been widely adopted (Papageorgiou, Stylios & Groumpos, 2006;Stach, 2010). In the following sections, we present examples of these algorithms. Additionally, we successfully tested some of these learning methods on a case study about nutrition which included 257 participants. We refer readers to our work for more information and for an illustration of the application of the presented library on a real-world use case (Wozniak, Mkhitaryan & Giabbanelli, 2022).

Hebbian learning
One of the weaknesses of an FCM constructed by the experts is its potential convergence to undesired regions. For example, given an intervention scenario, the model may predict only extreme values such as 0 or 1 (Lavin et al., 2018). To overcome this weakness Papageorgiou, Stylios & Groumpos (2006) proposed two learning strategies, namely the Active Hebbian Learning (AHL) and the Non-Linear Hebbian Learning (NHL) algorithms that are based on the Hebbian learning rule (Hebb, 2005). The task of the proposed algorithms is to modify the initial FCM connection matrix constructed by the expert such that the chosen nodes (called Desired Output Concepts DOCs) always converge within the desired range. Both algorithms are similar to FCM simulation, with the main difference being that concepts' values and weights are updated at each time step, whereas during a simulation, only the concepts values are changing.
In the NHL algorithm, the activation values and weights are simultaneously updated at each time step. In AHL, nodes and weights are updated asynchronously based on a sequence of activation patterns specified by the user. During each simulation time step, a new node becomes an "activated node"; only this node and its incoming edges are updated, while everything else remains unchanged. Along with optimizing existing edges, AHL creates new connections between the concepts, which may be an undesirable behavior if the modeler's intent is to tweak the weights rather than create connections that have not been endorsed by experts.
The learning process continues until two termination conditions are fulfilled. First, the fitness function (F 1 ) is calculated for each DOC as shown in Eq. (15). If value of F 1 for each DOC declines at each time step, and the DOCs values are within a desired range, the first termination condition is fulfilled. Second, it is crucial to determine whether the values of the DOCs are stable, i.e., if their values vary with each step more than a threshold e shown in Eq. (16). This threshold should be determined experimentally, and it is recommended to be set between 0.001 and 0.005 (Papageorgiou, Stylios & Groumpos, 2006). If the change is lower than the threshold, the second termination condition is fulfilled.
If the termination conditions are satisfied then the learning process may stop, otherwise, it will continue until a maximum number of steps is reached (we set the default value to 100). In order to use these methods, the user has to provide the initial weight matrix, initial concept values, and the DOCs. In addition, these variables are necessary, in most cases, algorithms converge only for a specific combination of values of the hyperparameters: learning rate (g), decay coefficient (c), and the slope of the sigmoid function.  (Papageorgiou, Stylios & Groumpos, 2004;Ren, 2012;Papakostas et al., 2011) using the algorithms above is demonstrated in the code snippet below.
>>> from fcmpy import NHL >>> import numpy as np If the learning process was successful (i.e., the algorithm converged), the run method will return the optimized weight matrix as a dataframe. Successful outputs of NHL and AHL algorithms are presented in Fig. 7.

Real-coded genetic algorithm (RCGA)
In certain domains of application, one has longitudinal data about the state variables included in the FCM and wants to find an FCM connection matrix that generates data that is close enough to the collected data (Khan, Khor & Chong, 2004;Poczeta, Yastrebov & Papageorgiou, 2015). In this regards, Stach (2010) proposed a real coded genetic algorithm for searching for an optimal FCM connection matrix. The proposed algorithm, named RCGA, builds on genetic algorithms where the search process includes the following six steps: (1) initialization, (2) evaluation, (3) selection, (4) recombination, (5) mutation, and (6) replacement.
In the initialization step, the algorithm generates a population of random solutions, that is a set of random weight matrices. Each solution is an n × n connection matrix. In the evaluation step, each candidate solution in the population is evaluated based on a fitness function shown in Eq. (17) and Table 7. We can observe how the fitness function is calculated on a simple example, shown in Fig. 8.  Error ¼ a In the selection step, candidate solutions are selected for mating (a.k.a., recombination). In each step, the algorithm randomly selects between two selection mechanisms (roulette wheel and tournament selection strategies). For the recombination step, the algorithm implements the recommended one point crossover operation with a probability of crossover specified by the user (p_recombination) (Stach, 2010). The crossover operation creates new solutions based on the solutions selected in the previous step. Next, the algorithm decides whether the new solutions produced in the previous step should undergo mutations. In the mutation step, the algorithm chooses between random and non-uniform mutation operations with a probability defined by the user (p_mutation). The replacement step is determined by the evolutionary approach specified by the user. The algorithm proposed by Stach (2010) is based on a generational approach where in each  Table 8 Variables of fitness function used in RCGA algorithm.

Argument Option Description
Normalization_type Hyperbolic tangent step the new generation of solutions replace the old generation. Alternatively, the user could choose a steady state approach (a.k.a., SSGA) where in each step only two new solutions are produced and a decision is made whether the new chromosomes should be inserted back into the population. The current implementation of the SSGA uses a replacement strategy based on the concept of useful diversity (described in depth in Lozano, Herrera & Cano (2008)). To use the RCGA module, one needs to initialize the RCGA class by specifying the longitudinal data about the system, the population size, and the genetic approach to use (i.e., generational or steady state). The additional parameters that can be modified by the user are presented in Table 8. Other parameters that can be modified by the user can be found in the documentation of the package available on PyPI.
The output of the learning process is the weight matrix with the highest fitness value throughout the search process. An example of generating FCM by the RCGA using historical data on water tank case study (Papageorgiou, Stylios & Groumpos, 2004) presented in the previous section is demonstrated in the code snippet below. We give the user an option to choose: population_size which is a number of weights matrices generated at each time step and threshold that is a minimum fitness value of at least one weight matrix in a generation, for the algorithm to succeed. To ensure the user does not get stuck in an infinite loop, we define n_iterations after which algorithm will terminate if the max fitness function of the n iterations −1 generation was less than a threshold.  The learned FCM connection matrix can be validated by calculating the in-sample and out-sample errors by using the dedicated ISE and OSE modules (see the code snippet below) (Stach, 2010

Classification algorithms
It is attractive for the researchers to choose FCMs for classification tasks, over other popular tools such as neural networks. This is because FCMs are easily explainable, which is a great advantage over black box models and, in many cases, equally accurate (Nápoles, Jastrzębska & Salgueiro, 2021;Nápoles et al., 2020). We give user a choice of two methods: Evolving Long-term Cognitive Networks (ELTCN) (Nápoles, Jastrzębska & Salgueiro, 2021) and deterministic learning (LTCN-MP) 1 (Nápoles et al., 2020). LTCN-MP and ELTCN use the same topology (a fully connected FCM containing features nodes and class nodes) but the former produces numerical outputs (suitable for regression) while the latter produces nominal outputs (suitable for classification). In LTCN-MP and ELTCN algorithms (i) input variables are located in the inner layer and output variables in the outer layer, (ii) weights connecting the inputs are computed in an unsupervised way by solving a least squared problem, and (iii) weights connecting inputs with outputs are computed using the Moore-Penrose pseudo-inverse. Overall, we can say that LTCN-MP and ELTCN use the same topology but the former produces numerical outputs while the latter produces nominal outputs (decision classes). Additionally, the weights in the ELTCN model can change from one iteration to another.
An example of a model's structure with three features and three classes is shown in Fig. 9. The user has to provide the path to the directory where the data file (.arff format) is located 2 . It is necessary that values of the features are normalized in the range between 0 and 1. Multiple data sets can be utilized and the results of each of them is saved in a dictionary, using the filename as a key. After running the learning process, the output consist of k weight matrices, where k is a number of validation folds (default 5) and the weight matrix of the connections between classes and feature nodes. We also provide users with automatically generated histograms showing values of the loss function and weights 1 LSTCNs are a variant of FCMS where weights are not expected to be in the [-1,1] interval or have a causal meaning. We use regularization in order to keep them within that range.
for each fold.
In our examples we are using the Iris data set, a popular data set containing various measurements of iris flowers, such as sepal or petal length (Fisher, 1936 LSTCN-MP focuses on discovering which and how features of the data set are important for the classification task as well as finding the weights connecting the inputs and outputs. Next, the LSTCN-MP algorithm outputs a 1-D array with values in the [-1,1] range. The absolute values represent how important the features are for the classification task. In order to use the algorithm, the user has to provide the path of data sets as list under a key sources and then use that dictionary as an input to LSTCN-MP algorithm.   The function returns a list of dictionaries (one dictionary for each input data set) with keys representing hyperparameters used in the learning process, weight matrix, training error and importance of the features, which is shown in Fig. 10. Note that in Iris dataset, feature 0 is the most important feature whereas feature 1 is the least significant one in the decision-making.
In this section, we presented various implementations of learning algorithms that can be used to train FCMs based on expert inputs or quantitative data available about the system and solve classification problems. In the next section, we illustrate how FCMs can be used to analyze scenarios.

SCENARIO ANALYSIS WITH FCMS
Scenario analysis in an FCM framework is often implemented by either changing the baseline values of the concepts (single shot interventions) or by introducing the proposed scenario as a new factor in the defined FCM and specifying the causal impact the proposed intervention has on the target concepts (continuous interventions). The single shot interventions mimic interventions that stop when a desired change in the specific target variables are achieved. In the continuous case, the intervention becomes part of the system and continuously impacts the target variables (Giabbanelli & Crutzen, 2014). The Intervention module provides the respective methods for analyzing different intervention cases. The module is instantiated by passing a simulator object to the constructor.

>>> from fcmpy import FcmSimulator, FcmIntervention >>> inter = FcmIntervention(FcmSimulator)
Before specifying intervention cases and running simulations for each scenario, we need to create the baseline for the comparison (i.e., run a simulation with baseline initial conditions and take the final state vector). To do this one needs to call initialize method. As in the FcmSimulator presented in the previous section, one can specify the output concepts by supplying a list of these concepts to the respective output_concepts argument. If the output_concepts argument is not specified then all the concepts in the FCM are treated as output concepts and the simulation stops when all the concepts change by less than the threshold between two consecutive steps.  We can inspect the results of the initial simulation run (i.e., 'baseline') in the test_results field as follows: We can use the add_intervention method to specify the intervention cases. To specify a single shot intervention we must specify the name of the intervention and supply new initial states for the concept values as a dictionary (see the code snippet below). For continuous intervention cases we must specify the name of the intervention, the concepts the intervention targets and the impact the intervention has on these concepts. In some cases we might be interested in checking scenarios where the intervention fails to be delivered to its fullest. For such cases we can specify the effectiveness of a given intervention case by setting the (optional) effectiveness argument to a number in the [0,1] interval (see the code snippet below). The effectiveness will decrease the expected causal strength of the intervention: for example, if an intervention is expected to reduce stress by 0.5 but is only 20% effective, then its actual reduction will be 0.1. >>> inter.add_intervention('intervention_1', type='continuous', impact={'C1':-.3, 'C2':.5}, effectiveness=1) >>> inter.add_intervention('intervention_2', type='continuous', impact={'C4':-.5}, effectiveness=1) >>> inter.add_intervention('intervention_3', type='continuous', impact={'C5':-1}, effectiveness=1) In the example above, we specify three intervention cases. The first intervention targets concepts (nodes) C1 and C2. It negatively impacts concept C1 (-0.3) while positively impacting the concept C2 (0.5). We consider a case where the intervention has maximum effectiveness. The other two interventions follow the same logic but impact other nodes.
After specifying the proposed interventions, we can use the test_intervention method to test the effect of each case. The method requires the name of the intervention to be tested. Users also have the possibility of changing the number of iterations for the simulation; its default value is the same as specified in the initialization (see the code snippet below).

EXTENSIBILITY
Given the fact that new algorithms are continuously developed, the extensibility of the package is of paramount importance. Each module in the package uses a defined interface which ensures the scalability and cohesion of the package and their future extensions. Furthermore, we separated the creation of the objects from their use, to ensure that the future extensions do not cause major changes to the user code. Let us examine how the package can be extended by considering a case where we want to add a new feature that allows us to read.txt file format. First, in the reader.py file we would need to define a new class called TXT which implements the ReadData interface. The interface has an abstract method called read.
Subsequently, we would need to add this object to the file called methodsStore which stores all the classes that support different file formats (e.g., CSV, XLSX, JSON). The ReadStore class has a private attribute called __methods which is a dictionary that stores the argument we want the user to call in the read_data method as a key and the value is the Class we created for this particular file extension. The get method of the ReaderStore class takes the file extension as an argument and returns the corresponding object from the __methods dictionary. This ReaderStore.get method is called inside the fcmpy.ExpertFcm.read_data method, during which the appropriate file reader is supplied based on the user's specification.

CONCLUSION
We hope that by providing researchers with such a tool we will promote using FCMs as one of the potential methods for different engineering tasks. We showed that FCMPy is a valuable tool for creating transparent and explainable behavioral models based on experts' answers. We provided detailed examples of how to create and inspect the model using different membership functions. Since one of the main purposes of FCMs is to monitor how the values of the concepts change throughout time, we provided the users with all the necessary options such as different inference methods and transfer functions. As a lot of research in the FCM field focuses on machine learning, we added several algorithms for weights optimization and data-driven model generation and made them easy to use. Finally, our library allows researchers to effortlessly examine how different interventions influence their model.
The FCMpy package provides a complete set of functions necessary to conduct projects involving FCMs. We created a tool that is open-source, easy to use, and provides the necessary functionality. The design and implementation of the tool results from a collaboration with multiple experts from the field of FCMs. We believe that this tool will facilitate research and encourage new students and scientists to involve FCMs in their projects. We included both well-known algorithms as well as recently developed ones. We