added mutual info plot

## Correlations

### Cross-correlation

The [Pearson product-moment correlation][2] was calculated for each pair of variables:

![Cross-correlation](img/7-pearson-correlation.png)

### Mutual-information

Normalized [mutual information][3] was calculated for each pair of variables. A value of 0 indicated no relationship. A value of 1 indicated perfect predictive ability between the two variables.

![Mutual information (MI)](img/7-mutual-information.png)

### Fan speed and power consumption

*Hypothesis*: Fan power, `PowFan[A | B]` depends on ambient temperature `TempAmbient`, relative humidity `PerHumidity`, and fan speed setting `PerFreqFan[A | B]`.

[1]: http://scikit-learn.org/stable/modules/clustering.html#dbscan

[2]: https://en.wikipedia.org/wiki/Pearson_correlation_coefficient?oldformat=true

[3]: https://en.wikipedia.org/wiki/Mutual_information

 No newline at end of file

%% Cell type:code id: tags:

``` python

%matplotlib notebook

%matplotlib inline

%reload_ext autoreload

%autoreload 2

import datetime

from os import path, environ

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

from sklearn.preprocessing import MinMaxScaler

from sklearn.preprocessing import MinMaxScaler, StandardScaler

from sklearn.pipeline import Pipeline

from sklearn.linear_model import LinearRegression

from sklearn.model_selection import train_test_split

from sklearn.cluster import DBSCAN

from sklearn.metrics import mutual_info_score

import ipyvolume as ipv

from thermo import CONSTANTS

# source file, see docs/5-dataset.md for info on field names

Chiller1File = path.join(environ['DATADIR'], 'EngineeringScienceBuilding', 'Chillers.csv')

plot_path = path.join('..', 'docs', 'img')

```

%% Cell type:code id: tags:

``` python

# Read pre-processed data

df = pd.read_csv(Chiller1File, index_col='Time', parse_dates=True, dtype=float)

```

%% Cell type:markdown id: tags:

# Correlations

## Cross-correlation

Pearson product-moment correlation coefficients between variables. Measure of linear relationship between variables.

$$

r_{X,Y} = \frac{\textrm{cov}(X,Y)}{\sigma_X \sigma_Y}

$$

%% Cell type:code id: tags:

``` python

corr = np.corrcoef(df, rowvar=False)

fig, ax = plt.subplots(figsize=(10, 10))

axImg = ax.matshow(corr, vmin=-1, vmax=1)

fig.colorbar(axImg)

plt.xticks(np.arange(len(df.columns)), df.columns, rotation='vertical')

plt.yticks(np.arange(len(df.columns)), df.columns);

plt.savefig(path.join(plot_path, '7-pearson-correlation.png'))

```

%% Cell type:markdown id: tags:

### Mutual information

Mutual information is a measure of how much information about the distribution variable $X$ is contained in the distribution of variable $Y$.

$$

\mathrm{MI}(X,Y) = \sum_{i}^{\mid X \mid} \sum_{j}^{\mid Y \mid} \frac{\mid X_i \cap Y_j \mid}{N} \log{\frac{N \; \mid X_i \cap Y_j \mid}{\mid X_i \mid \; \mid Y_j \mid}}

$$

Where $N$ is the total number of samples, $i$ and $j$ are class labels for samples (in this case, histogram bins). $X_i$ is the number of samples of $X$ with a label $i$. The MI score is normalized to lie between 0 and 1.

%% Cell type:code id: tags:

``` python

bins = 32

mi = np.empty((len(df.columns), len(df.columns)), dtype=float)

# Standardizing values for 0 mean and unit variance

scaled = StandardScaler().fit_transform(df)

# Calculating MI matrix for upper triangular half explicitly,

# the lower half is symmetric.

for i, col1 in enumerate(df.columns):     # row

    for j, col2 in enumerate(df.columns): # column

        if j < i:

            continue

        else:

            hist,_ , _ = np.histogram2d(scaled[:,i], scaled[:,j], bins=bins, density=False)

            N = hist.sum()

            # Normalize MI, taken from source of sklearn's

            # normalized_mutual_info_score to work with a

            # contingency table.

            p_i, p_j = np.sum(hist, axis=0) / N, np.sum(hist, axis=1) / N

            non_zero_i, non_zero_j = p_i > 0, p_j > 0

            ent_i = -np.sum(p_i[non_zero_i] * np.log(p_i[non_zero_i]))

            ent_j = -np.sum(p_j[non_zero_j] * np.log(p_j[non_zero_j]))

            norm = 0.5 * (ent_i + ent_j)

            mi[i,j] = mutual_info_score(None, None, contingency=hist) / norm

            mi[j,i] = mi[i,j]

fig, ax = plt.subplots(figsize=(10, 10))

axImg = ax.matshow(mi)

fig.colorbar(axImg)

plt.xticks(np.arange(len(df.columns)), df.columns, rotation='vertical')

plt.yticks(np.arange(len(df.columns)), df.columns);

plt.savefig(path.join(plot_path, '7-mutual-information.png'))

```

%% Cell type:markdown id: tags:

## Fan speed & fan power

*Hypothesis*: Fan power depends on ambient temperature `TempAmbient`, relative humidity `PerHumidity`, and fan speed setting `PerFreqFan[A | B]`.

*Result*: Very little relationship from current data. Primarily because fan speed is ~100%. Need more variation.

%% Cell type:code id: tags:

``` python

# Data

X = df.loc[:, ('TempAmbient', 'PerHumidity', 'PerFreqFanA')]

Y = df.loc[:, 'PowFanA']

Xtrain, Xtest, Ytrain, Ytest = train_test_split(X, Y, test_size=0.1)

```

%% Cell type:code id: tags:

``` python

# Pipeline

scaler = MinMaxScaler()        # scale to [0-1] range

regressor = LinearRegression() # regression model

reg = Pipeline([('scaler', scaler), ('regressor', regressor)])

```

%% Cell type:code id: tags:

``` python

reg.fit(Xtrain, Ytrain)

reg.score(Xtest, Ytest)

```

%% Output

    0.1107982030593837

%% Cell type:code id: tags:

``` python

f = plt.figure(figsize=(8,8))

ax = f.add_subplot(111, projection='3d')

ax.scatter3D(X.iloc[:, 0], X.iloc[:, 1], Y, c=Y)

ax.set_xlabel('Ambient Temperature (K)')

ax.set_ylabel('Relative Humidity')

ax.set_zlabel('Fan Power (W)')

plt.savefig(path.join(plot_path, '7-fan-power-vs-temp-humidity.png'))

```

%% Cell type:markdown id: tags:

# Clusters

## Temperature

%% Cell type:code id: tags:

``` python

X = df.loc[:, ('TempAmbient', 'TempWetBulb')]

X['DeltaTemp'] = df['TempCondOut'] - df['TempCondIn']

labels = DBSCAN().fit_predict(X)

```

%% Cell type:code id: tags:

``` python

ipv.clear()

cmap = plt.cm.Accent(labels)

f = ipv.scatter(X['TempAmbient'], X['DeltaTemp'], X['TempWetBulb'], color=cmap,

            size=10, marker='point_2d')

ipv.xyzlim(280, 310)

ipv.ylim(0, 10)

ipv.xyzlabel('TempAmbient', 'DeltaTemp', 'TempWetBulb')

ipv.show()

```

%% Cell type:code id: tags:

``` python

def view(fig, n, frac):

    ipv.view(360*frac)

ipv.movie(path.join(plot_path, '7-ct-temp-clusters.gif'), view, fps=12, frames=36)

```

%% Output

%% Cell type:markdown id: tags:

![](../docs/img/7-ct-temp-clusers.gif)

%% Cell type:markdown id: tags:

## Power

%% Cell type:code id: tags:

``` python

X = df.loc[:, ('PowConP', 'PowFanA', 'PowFanB')]

labels = DBSCAN().fit_predict(X)

```

%% Cell type:code id: tags:

``` python

ipv.clear()

cmap = plt.cm.Accent(labels)

f = ipv.scatter(X['PowConP'], X['PowFanA'], X['PowFanB'], color=cmap,

            size=10, marker='point_2d')

# ipv.xyzlim(280, 310)

ipv.xyzlabel('PowConP', 'PowFanA', 'PowFanB')

ipv.show()

```

%% Cell type:code id: tags:

``` python

def view(fig, n, frac):

    ipv.view(360*frac)

ipv.movie(path.join(plot_path, '../docs/img/7-ct-power-clusters.gif'), view, fps=12, frames=36)

```