#!/usr/bin/env python
# coding: utf-8

# # Fuzzy Logic Quality Control

# ### Example of the fuzzy approach to quality control
# 
# The source code of this ipython notebook, as well as other examples, can be found in the documentation of CoTeDe (http://cotede.castelao.net)

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')


# In[2]:


import numpy as np
import skfuzzy as fuzz
import matplotlib.pyplot as plt


# ## Introduction

# 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:
# 
# - Spike
# - Climatology
# - Gradient
# 
# We'll assume that the measurement $x_i$ resulted in the following values:

# In[3]:


spike = 1.0
clim = 5.4
grad = 0.9


# ## Membership Functions

# The first step to do is to define the membership functions on 3 levels of uncertainty: low, medium and high, for each feature. 

# In[4]:


from cotede.utils import load_cfg
cfg = load_cfg('fuzzylogic')['TEMP']
cfg.keys()


# An example of membership function parameters. The fuzzy set 'low' for feature 'spike' is:

# In[5]:


cfg['fuzzylogic']['features']['spike']


# In[6]:


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.

# In[7]:


# 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()


# ### Fuzzy Membership

# In[8]:


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


# ##  Fuzzy Rules

# Now we take our rules and apply them.

# ### Low uncertainty

# In[9]:


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)


# ### Medium uncertainty

# In[10]:


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)


# ### High uncertainty

# In[11]:


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


# In[12]:


# 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()


# ## Defuzzifying

# In[13]:


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


# In[14]:


# 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()


# In[ ]: