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

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.

>>> from fcmpy import ExpertFcm
>>> fcm = ExpertFcm()
>>> fcm.linguistic_terms = {
                ’−vh’: [−1, −1, −0.75],
               ’−h’: [−1, −0.75, −0.50],
               ’−m’: [−0.75, −0.5, −0.25],
               ’−l’: [−0.5, −0.25, 0],
               ’−vl’: [−0.25, 0, 0],
               ’na’: [−0.001, 0, 0.001],
               ’+vl’: [0, 0, 0.25],
               ’+l’: [0, 0.25, 0.50],
               ’+m’: [0.25, 0.5, 0.75],
               ’+h’: [0.5, 0.75, 1],
               ’+vh’: [0.75, 1, 1]
               }
>>> data = fcm.read_data (file_path = ’data_test.xlsx’,
               check_consistency = False, engine = ’openpyxl’)

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:

(1) $$R=-\sum _{i=1}^n{p}_{i}lo{g}_2(p_i)$$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).

>>> entropy = fcm.entropy (data)
From    To  Entropy
C1     C1  0.000000
    C2  1.459148
    C3  0.000000
    C4  0.000000
C2     C1  1.459148
    C2  0.000000
    C3  0.000000
    C4  0.000000
C3     C1  0.820802
    C2  0.000000
    C3  0.000000
    C4  0.930827
C4     C1  0.000000
    C2  0.000000
    C3  0.000000
    C4  0.000000

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).

>>> fcm.linguistic_terms = {
                ’−VH’: [−1, −1, −0.75],
               ’−H’: [−1, −0.75, −0.50],
               ’−M’: [−0.75, −0.5, −0.25],
               ’−L’: [−0.5, −0.25, 0],
               ’−VL’: [−0.25, 0, 0],
               ’No Causality’: [−0.001, 0, 0.001],
               ’+VL’: [0, 0, 0.25],
               ’+L’: [0, 0.25, 0.50],
                ’+M’: [0.25, 0.5, 0.75],
                ’+H’: [0.5, 0.75, 1],
                ’+VH’: [0.75, 1, 1]
                }

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:

(2) $${\mu }_R(x,y)=min\left({\mu }_A(x),{\mu }_B(y)\right).$$where ${\mu }_A(x)$ and ${\mu }_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) $${\mu }_R(x,y)={\mu }_A(x)\cdot {\mu }_B(y)$$

We can use fuzzy_implication method to apply the selected implication method (see the available methods in Table 3 and the code snippet below).

>>> mfs = fcm.fuzzy_membership
>>> act_pvh = fcm.fuzzy_implication(mfs[’+vh’],
                weight = 0.33, method = ’Mamdani’)
>>> act_pm = fcm.fuzzy_implication(mfs[’+m’],
                weight = 0.16, method = ’Mamdani’)
>>> act_ph = fcm.fuzzy_implication(mfs[’+h’],
                weight = 0.5, method = ’Mamdani’)
>>> activatedMamdani = {’+vh’: act_pvh,
                 ’+h’: act_ph, ’+m’: act_pm}

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).

(4) $$f(x,y)=x+y-x\cdot y$$

(5) $$f(x,y)={\displaystyle \frac{(x+y)}{(1+x\cdot y)}}$$

(6) $$f(x,y)={\displaystyle \frac{(x+y-2\cdot x\cdot y)}{(1-x\cdot y)}}$$

One can use the aggregate method to aggregate the activated membership functions (see the available aggregation method in Table 4 and the code snippet below).

>>> import functools
>>> aggregated = functools.reduce(lambda x,y:
           fcm.aggregate(x=x, y=y, method=’fMax’),
        [activatedMamdani[i] for i in activatedMamdani.keys()])

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.

>>> weight_matrix = fcm.build(data=data,
          implication_method = ’Mamdani’,
aggregation_method = ’fMax’,
defuzz_method=’centroid’)
     C1   C2    C3    C4
C1  0.000000 0.702032 0.000000 0.000000
C2  0.607453 0.000000 0.000000 0.000000
C3  0.555914 0.000000 0.000000 0.172993
C4  0.000000 0.000000 0.000000 0.000000

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.

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.

(15) $$F_1=\sqrt{|DO{C}_j^{(k)}-{\displaystyle \frac{DO{C}_j^{min}-DO{C}_j^{max}}{2}{|}^2}}$$

(16) $$F_2=|DO{C}_j^{(k+1)}-DO{C}_j^{(k)}|<e$$

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 ( $\eta$), decay coefficient ( $\gamma$), and the slope of the sigmoid function. The sample values used in several case studies are slope [0.9, 1.01], decay for NHL [0.99, 1.0], decay for AHL [0.01, 0.1] and learning rate [0.001, 0.1]. The optimization of an FCM from a water tank case study (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
# initial values of weight matrix
>>> w_init_WT = np.asarray([[0,−0.4,−0.25,0,0.3],
              [0.36,0,0,0,0],
              [0.45,0,0,0,0],
              [−0.9,0,0,0,0],
              [0,0.6,0,0.3,0]])
>>> w_init_WT = pd.DataFrame(w_init_WT,
                columns=[’C1’, ’C2’, ’C3’, ’C4’, ’C5’],
               index = [’C1’, ’C2’, ’C3’, ’C4’, ’C5’])
# initial values of the concepts
>>> init_states_WT = {’C1’: 0.40, ’C2’: 0.7077,
           ’C3’: 0.612, ’C4’: 0.717, ’C5’: 0.30}
# DOCs
>>> doc_values_WT = {’C1’:[0.68,0.74], ’C5’:[0.74,0.8]}
# NHL
>>> nhl = NHL(state_vector=init_states_WT,
       weight_matrix=w_init_WT, doc_values=doc_values_WT)
>>> res_nhl = nhl.run(learning_rate = 0.01, l=.98, iterations=100)
The NHL learning process converged at step 63 with the
learning rate eta = 0.01 and decay = 1!
   C1     C2     C3     C4    C5
C1   0.000000 −0.200310 −0.023806  0.000000  0.472687
C2   0.539068   0.000000   0.000000  0.000000  0.000000
C3   0.571531   0.000000   0.000000  0.000000  0.000000
C4 −0.832174   0.000000   0.000000  0.000000  0.000000
C5   0.000000   0.710523   0.000000  0.496934  0.000000

The AHL.run method has an additional auto_learn argument; if set to True then the algorithm automatically updates the hyperparameters during the learning process.

# AHL
>>> activation_pattern_WT = {0:[’C1’], 1:[’C2’, ’C3’],
                2: [’C5’], 3: [’C4’]}
>>> ahl = AHL(state_vector=init_states_WT, weight_matrix=w_init_WT,
         activation_pattern=activation_pattern_WT,
         doc_values=doc_values_WT)
>>> res_ahl = ahl.run(decay=0.03, learning_rate = 0.01, l=1,
              iterations=100, transfer= ’sigmoid’,
              thresh = 0.002, auto_learn=False,
              b1=0.003, lbd1=0.1, b2=0.005, lbd2=1)
The AHL learning process converged at step 19 with
the learning rate eta = 0.01 and decay = 0.03!
    C1    C2    C3     C4    C5
C1   0.000000 −0.128532 −0.060395  0.071200  0.218170
C2   0.245859  0.000000   0.068981  0.076592  0.074289
C3   0.288257  0.069457   0.000000  0.070342  0.068190
C4 −0.386349  0.073807   0.067187  0.000000  0.073991
C5   0.070113  0.368913   0.069145  0.223312  0.000000

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.

(17) $$\begin{array}{rl}& Error=\alpha \sum _{t=1}^{T-1}\sum _{j=1}^{N-1}|A_j(t)-{\hat{A}}_j(t){|}^p\\ & Fitness={\displaystyle \frac{1}{a\ast Error+1}}\end{array}$$

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 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\mathrm{\_}iterations-1$ −1 generation was less than a threshold.

# Generate Longitudinal Data
>>> sim = FcmSimulator()
>>> data_WT = sim.simulate(initial_state=init_states_WT,
             weight_matrix=w_init_WT, transfer=’sigmoid’,
            inference=’mKosko’, thresh=0.001,
            iterations=50, l=1)
# Select two time points
>>> data_WT = data_WT.iloc[:3]
# Generational Approach
>>> rcga = RCGA(data=data_WT, population_size=100,
             ga_type=’generational’)
>>> rcga.run(n_iterations=30000, threshold=0.99)
# Steady State Approach
>>> rcga = RCGA(data=data_WT, population_size=100,
             ga_type=’ssga’)
>>> rcga.run(n_iterations=30000, threshold=0.99)
The RCGA solution and the associated fitness score can be accessed in the rcga.solution and rcga.fitness fields.

>>> rcga.fitness
0.9790443073348711
>>> rcga.solution
   C1     C2    C3     C4     C5
C1   0.851948   0.550687  0.045434   0.648337  0.077349
C2 −0.870796 −0.944261 −0.684952 −0.829811   0.660863
C3   0.807396 −0.075391  0.223011   0.299976 −0.442272
C4 −0.479566   0.724956  0.048086   0.352487   0.120083
C5 −0.030473   0.068480  0.330551   0.003072 −0.347919

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).

>>> from fcmpy import ISE
>>> from fcmpy import OSE
>>> val_ise = ISE()
>>> val_ose = OSE()
>>> error_ise = val_ise.validate(initial_state=init_states_WT,
            weight_matrix=rcga.solution, data=data_WT,
            transfer=’sigmoid’, inference=’mKosko’, l=1)
>>> error_ose, std = val_ose.validate(weight_matrix=rcga.solution,
              data=data_WT, low=0, high=1,
              k_validation=100, transfer=’sigmoid’,
              inference=’mKosko’, l=1)

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)¹ (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² . 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 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).

>>> from fcmpy.ml.classification.eltcn import run
>>> path = ’data’
>>> results = run(path)

# average of the feature weight matrices
>>> print(results[’irisnorm.arff’][’avgW’])

# weight matrix connecting feature nodes with class nodes
>>> print(results[’irisnorm.arff’][’classW’])

[[ 0.4382333, −0.09812646, 0.79643613, 0.97163904],
[−0.2016794, 0.2925336, −0.21699643, 0.12452463],
[−0.0038293, −0.08575615, 0.47770625, 0.52030253],
[ 0.24917993, 0.21717176, 0.30113676, 0.34600133]]

[[−0.9135318, 0.15331155, 0.24873467],
[−1., 0.42972276, 0.17217618],
[−0.5610481, −0.58390266, 0.79575956],
[−0.60481995, −0.6718271, 0.72851056]]

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.

>>> sources = [‘datasets/iris.arff’,‘datasets/pima.arff’]
>>> params = {‘sources’:sources}

Various hyperparameters can be used. Their default values and description can be found in Nápoles et al. (2020), our library documentation, or Table 9.

>>> import fcmpy.ml.classification.FCM_MP as mp
>>> import matplotlib.pylab as plt
>>> sources = [‘iris.arff’]
>>> params = {‘sources’:sources}
>>> out = mp.run(**params)
# feature importance for classification purposes
>>> fig, ax = plt.subplots()
>>> ax.bar(range(len(out[0][’importance’].flatten())),
height=out[0][’importance’].flatten())
# connections between features and class nodes
>> print(out[‘weights’])
[[ 0.46838854, −0.03855411, 1. ],
[−0.04176139, 0.46838854, −0.48478635],
[ 0.16611878, −0.07434725, 0.46838854]]

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.