The source code of this ipython notebook, as well as other examples, can be found in the documentation of CoTeDe (http://cotede.castelao.net)
%matplotlib inline
import numpy as np
import skfuzzy as fuzz
import matplotlib.pyplot as plt
Let's assume a set of measurements $x$ that could be temperature measurements from an ARGO or salinity from a CTD. The goal here is to evaluate the quality of each measurement ($x_i$) using Fuzzy Logic. We would like to evaluate each measurement $x_i$ not only by the value itself, but from different perspectives, like how similar is that measurement to the climatology? We'll call each one of this perspectives as a feature, and will construct a multi-dimensional evaluation system.
The features we'll use here are:
We'll assume that the measurement $x_i$ resulted in the following values:
spike = 1.0
clim = 5.4
grad = 0.9
The first step to do is to define the membership functions on 3 levels of uncertainty: low, medium and high, for each feature.
from cotede.utils import load_cfg
cfg = load_cfg('fuzzylogic')['TEMP']
cfg.keys()
[u'gradient', u'woa_normbias', u'fuzzylogic', u'spike']
An example of membership function parameters. The fuzzy set 'low' for feature 'spike' is:
cfg['fuzzylogic']['features']['spike']
{u'high': [2, 6], u'low': [0.07, 0.2], u'medium': [0.07, 0.2, 2, 6], u'weight': 1}
data = {}
data['x_spike'] = np.linspace(0, 7, 200)
data['x_clim'] = np.linspace(0, 10, 200)
data['x_grad'] = np.linspace(0, 6, 200)
x_QC = np.linspace(0.0, 1.0, 200)
# Generate fuzzy membership functions
data['spike_lo'] = fuzz.zmf(data['x_spike'],
cfg['fuzzylogic']['features']['spike']['low'][0],
cfg['fuzzylogic']['features']['spike']['low'][1])
data['spike_md'] = fuzz.pimf(data['x_spike'],
cfg['fuzzylogic']['features']['spike']['medium'][0],
cfg['fuzzylogic']['features']['spike']['medium'][1],
cfg['fuzzylogic']['features']['spike']['medium'][2],
cfg['fuzzylogic']['features']['spike']['medium'][3]
)
data['spike_hi'] = fuzz.smf(data['x_spike'],
cfg['fuzzylogic']['features']['spike']['high'][0],
cfg['fuzzylogic']['features']['spike']['high'][1])
data['clim_lo'] = fuzz.zmf(data['x_clim'],
cfg['fuzzylogic']['features']['woa_normbias']['low'][0],
cfg['fuzzylogic']['features']['woa_normbias']['low'][1])
data['clim_md'] = fuzz.pimf(data['x_clim'],
cfg['fuzzylogic']['features']['woa_normbias']['medium'][0],
cfg['fuzzylogic']['features']['woa_normbias']['medium'][1],
cfg['fuzzylogic']['features']['woa_normbias']['medium'][2],
cfg['fuzzylogic']['features']['woa_normbias']['medium'][3]
)
data['clim_hi'] = fuzz.smf(data['x_clim'],
cfg['fuzzylogic']['features']['woa_normbias']['high'][0],
cfg['fuzzylogic']['features']['woa_normbias']['high'][1])
data['grad_lo'] = fuzz.zmf(data['x_grad'],
cfg['fuzzylogic']['features']['gradient']['low'][0],
cfg['fuzzylogic']['features']['gradient']['low'][1])
data['grad_md'] = fuzz.pimf(data['x_grad'],
cfg['fuzzylogic']['features']['gradient']['medium'][0],
cfg['fuzzylogic']['features']['gradient']['medium'][1],
cfg['fuzzylogic']['features']['gradient']['medium'][2],
cfg['fuzzylogic']['features']['gradient']['medium'][3]
)
data['grad_hi'] = fuzz.smf(data['x_grad'],
cfg['fuzzylogic']['features']['gradient']['high'][0],
cfg['fuzzylogic']['features']['gradient']['high'][1])
QC_lo = fuzz.trimf(x_QC, [0.0, 0.225, 0.45])
QC_md = fuzz.trimf(x_QC, [0.275, 0.5, 0.725])
QC_hi = fuzz.trimf(x_QC, [0.55, 0.775, 1.0])
QC_hi = fuzz.trapmf(x_QC, [0.55, 0.775, 1.0, 1e30])
Let's visualize the membership functions, and the input features being evaluated.
# Visualize these universes and membership functions
fig, (ax0, ax1, ax2, ax3) = plt.subplots(nrows=4, figsize=(8, 9))
ax0.plot(data['x_spike'], data['spike_lo'], 'g', linewidth=1.5, label='Low')
ax0.plot(data['x_spike'], data['spike_md'], 'y', linewidth=1.5, label='Medium')
ax0.plot(data['x_spike'], data['spike_hi'], 'r', linewidth=1.5, label='High')
ax0.axvline(spike, color='k', linewidth=4, alpha=0.7)
ax0.set_title('Spike')
ax0.legend()
ax1.plot(data['x_clim'], data['clim_lo'], 'g', linewidth=1.5, label='Low')
ax1.plot(data['x_clim'], data['clim_md'], 'y', linewidth=1.5, label='Medium')
ax1.plot(data['x_clim'], data['clim_hi'], 'r', linewidth=1.5, label='High')
ax1.axvline(clim, color='k', linewidth=4, alpha=0.7)
ax1.set_title('Climatology')
ax1.legend()
ax2.plot(data['x_grad'], data['grad_lo'], 'g', linewidth=1.5, label='Low')
ax2.plot(data['x_grad'], data['grad_md'], 'y', linewidth=1.5, label='Medium')
ax2.plot(data['x_grad'], data['grad_hi'], 'r', linewidth=1.5, label='High')
ax2.axvline(grad, color='k', linewidth=4, alpha=0.7)
ax2.set_title('Gradient')
ax2.legend()
ax3.plot(x_QC, QC_lo, 'g', linewidth=1.5, label='Low')
ax3.plot(x_QC, QC_md, 'y', linewidth=1.5, label='Medium')
ax3.plot(x_QC, QC_hi, 'r', linewidth=1.5, label='High')
ax3.set_title('QC output')
ax3.legend()
# Turn off top/right axes
for ax in (ax0, ax1, ax2):
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
plt.tight_layout()
spike_level_lo = fuzz.interp_membership(data['x_spike'], data['spike_lo'], spike)
spike_level_md = fuzz.interp_membership(data['x_spike'], data['spike_md'], spike)
spike_level_hi = fuzz.interp_membership(data['x_spike'], data['spike_hi'], spike)
print("Spike: %s" % spike)
print("Low uncert.: %.2f, Medium uncert.: %.2f, High uncert.: %.2f" %
(spike_level_lo, spike_level_md, spike_level_hi))
clim_level_lo = fuzz.interp_membership(data['x_clim'], data['clim_lo'], clim)
clim_level_md = fuzz.interp_membership(data['x_clim'], data['clim_md'], clim)
clim_level_hi = fuzz.interp_membership(data['x_clim'], data['clim_hi'], clim)
print("")
print("Climatology bias: %s" % clim)
print("Low uncert.: %.2f, Medium uncert.: %.2f, High uncert.: %.2f" %
(clim_level_lo, clim_level_md, clim_level_hi))
grad_level_lo = fuzz.interp_membership(data['x_grad'], data['grad_lo'], grad)
grad_level_md = fuzz.interp_membership(data['x_grad'], data['grad_md'], grad)
grad_level_hi = fuzz.interp_membership(data['x_grad'], data['grad_hi'], grad)
print("")
print("Gradient: %s" % grad)
print("Low uncert.: %.2f, Medium uncert.: %.2f, High uncert.: %.2f" %
(grad_level_lo, grad_level_md, grad_level_hi))
Spike: 1.0 Low uncert.: 0.00, Medium uncert.: 1.00, High uncert.: 0.00 Climatology bias: 5.4 Low uncert.: 0.00, Medium uncert.: 0.68, High uncert.: 0.32 Gradient: 0.9 Low uncert.: 0.68, Medium uncert.: 0.32, High uncert.: 0.00
Now we take our rules and apply them.
active_rule1 = np.mean((spike_level_lo, clim_level_lo, grad_level_lo), axis=0)
print("active rule low: %s, %s, %s -> %s" %
(spike_level_lo, clim_level_lo, grad_level_lo, active_rule1))
# Now we apply this by clipping the top off the corresponding output
# membership function with `np.fmin`
QC_activation_lo = np.fmin(active_rule1, QC_lo)
#print("QC_activation_lo: %s" % QC_activation_lo)
active rule low: 0.0, 0.0, 0.679768187672 -> 0.226589395891
active_rule2 = np.mean((spike_level_md, clim_level_md, grad_level_md), axis=0)
print("active rule medium: %s, %s, %s -> %s" %
(spike_level_md, clim_level_md, grad_level_md, active_rule2))
#QC_activation_md = np.fmin(clim_level_md, QC_md)
QC_activation_md = np.fmin(active_rule2, QC_md)
#print("QC_activation_md: %s" % QC_activation_md)
active rule medium: 1.0, 0.678745486225, 0.320231812328 -> 0.666325766184
# The OR operator means we take the maximum of these.
active_rule3 = np.fmax(grad_level_hi, np.fmax(spike_level_hi, clim_level_hi))
print("active rule high: %s, %s, %s -> %s" %
(spike_level_hi, clim_level_hi, grad_level_hi, active_rule3))
QC_activation_hi = np.fmin(active_rule3, QC_hi)
#print("QC_activation_hi: %s" % QC_activation_hi)
active rule high: 0.0, 0.321254513775, 0.0 -> 0.321254513775
# Visualize this
fig, ax0 = plt.subplots(figsize=(8, 3))
QC0 = np.zeros_like(x_QC)
ax0.fill_between(x_QC, QC0, QC_activation_lo, facecolor='g', alpha=0.7, label='low')
ax0.plot(x_QC, QC_lo, 'g', linewidth=0.5, linestyle='--', )
ax0.fill_between(x_QC, QC0, QC_activation_md, facecolor='y', alpha=0.7, label='medium')
ax0.plot(x_QC, QC_md, 'y', linewidth=0.5, linestyle='--')
ax0.fill_between(x_QC, QC0, QC_activation_hi, facecolor='r', alpha=0.7, label='high')
ax0.plot(x_QC, QC_hi, 'r', linewidth=0.5, linestyle='--')
ax0.set_title('Output membership activity')
# Turn off top/right axes
for ax in (ax0,):
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
plt.tight_layout()
# Aggregate all three output membership functions together
aggregated = np.fmax(QC_activation_lo,
np.fmax(QC_activation_md, QC_activation_hi))
# Calculate defuzzified result
#QC = fuzz.defuzz(x_QC, aggregated, 'centroid')
QC = fuzz.defuzz(x_QC, aggregated, 'bisector')
QC_activation = fuzz.interp_membership(x_QC, aggregated, QC) # for plot
print("The QC was evaluated as: %.2f" % QC)
The QC was evaluated as: 0.53
# Visualize this
fig, ax0 = plt.subplots(figsize=(8, 3))
ax0.plot(x_QC, QC_lo, 'b', linewidth=0.5, linestyle='--', )
ax0.plot(x_QC, QC_md, 'g', linewidth=0.5, linestyle='--')
ax0.plot(x_QC, QC_hi, 'r', linewidth=0.5, linestyle='--')
ax0.fill_between(x_QC, QC0, aggregated, facecolor='Orange', alpha=0.7)
ax0.plot([QC, QC], [0, QC_activation], 'k', linewidth=4, alpha=0.9)
ax0.set_title('Aggregated membership and result (line)')
# Turn off top/right axes
for ax in (ax0,):
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.get_xaxis().tick_bottom()
ax.get_yaxis().tick_left()
plt.tight_layout()