2025-AKO/v2_unsupervised.py

# Csar Fdez, UdL, 2025
#  Unsupervised classification. Uses tslearn
#  https://tslearn.readthedocs.io/en/stable/index.html

# Be careful with v0_unsupervised  and all versions for supervised.
# because dataTrain is not stacke before create_sequences,  so, 
#  the sets are not aligned in time

# Same as v1 but changing data

import pandas as pd
import matplotlib.pyplot as plt
import datetime
import numpy as np
import os.path
from optparse import OptionParser
import copy
import pickle
from tslearn.clustering import TimeSeriesKMeans
from tslearn.neighbors import KNeighborsTimeSeries 
from collections import Counter

parser = OptionParser()
parser.add_option("-t", "--train", dest="train", help="Trains the models (false)", default=False, action="store_true")
parser.add_option("-p", "--plot", dest="plot", help="Plot Data (false)", default=False, action="store_true")
parser.add_option("-n", "--timesteps", dest="timesteps", help="TIME STEPS ", default=12)

(options, args) = parser.parse_args()


# data files arrays. Index:
# 0.  No failure
# 1.  Blocked evaporator
# 2.   Full Blocked condenser
# 3   Fan condenser not working
# 4.  Open door


NumberOfFailures=4  
datafiles={}
listofFacilitySetpoint=['5-26','5-18','5-22','30','32','34']

for j in listofFacilitySetpoint:
    datafiles[j]=[]
    for i in range(NumberOfFailures+1):
        datafiles[j].append([])

# Freezer, SP=-26
datafiles['5-26'][0]=['2025-01-25_5_','2025-01-26_5_','2025-01-29_5_','2025-01-30_5_','2025-01-31_5_','2025-02-01_5_','2025-02-02_5_','2025-02-03_5_','2025-02-04_5_']
datafiles['5-26'][1]=['2024-12-11_5_', '2024-12-12_5_','2024-12-13_5_','2024-12-14_5_','2024-12-15_5_']
datafiles['5-26'][2]=['2024-12-21_5_','2024-12-22_5_','2024-12-23_5_','2024-12-24_5_','2024-12-25_5_','2024-12-26_5_']
datafiles['5-26'][3]=['2024-12-28_5_','2024-12-29_5_','2024-12-30_5_','2024-12-31_5_','2025-01-01_5_'] 
datafiles['5-26'][4]=['2025-02-13_5_','2025-02-14_5_','2025-02-15_5_','2025-02-16_5_','2025-02-17_5_','2025-02-18_5_','2025-02-19_5_'] 

# Freezer, SP=-18
datafiles['5-18'][0]=['2025-01-21_5_','2025-01-22_5_','2025-01-23_5_'] # 
datafiles['5-18'][1]=['2025-02-17_5_','2025-02-18_5_','2025-02-19_5_'] # 
datafiles['5-18'][2]=['2025-03-10_5_','2025-03-11_5_','2025-03-12_5_'] # 
datafiles['5-18'][3]=['2025-01-04_5_','2025-01-05_5_','2025-01-06_5_','2025-01-07_5_'] # 
datafiles['5-18'][4]=['2025-04-30_5_','2025-05-01_5_','2025-05-02_5_','2025-05-03_5_','2025-05-04_5_','2025-05-05_5_'] #   

# Freezer, SP=-22
datafiles['5-22'][0]=['2025-03-13_5_','2025-03-14_5_','2025-03-15_5_','2025-03-16_5_']  
datafiles['5-22'][1]=['2025-03-21_5_','2025-03-22_5_','2025-03-23_5_','2025-03-24_5_','2025-03-25_5_']  #  
datafiles['5-22'][2]=['2025-03-26_5_','2025-03-27_5_','2025-03-28_5_']  
datafiles['5-22'][3]=['2025-03-31_5_','2025-04-01_5_','2025-04-02_5_','2025-04-03_5_']  
datafiles['5-22'][4]=['2025-03-17_5_']  

# Refrigerator  0  
datafiles['30'][0]=['2025-01-21_3_','2025-01-22_3_','2025-01-23_3_','2025-01-24_3_','2025-01-25_3_','2025-01-26_3_']
datafiles['30'][1]=['2024-12-11_3_','2024-12-12_3_','2024-12-13_3_','2024-12-14_3_','2024-12-15_3_']
datafiles['30'][2]=['2024-12-18_3_','2024-12-19_3_','2024-12-20_3_']
datafiles['30'][3]=['2024-12-28_3_','2024-12-29_3_','2024-12-30_3_','2024-12-31_3_','2025-01-01_3_']
#datafiles['30'][4]=['2025-02-12_3_','2025-02-13_3_','2025-02-14_3_','2025-02-15_3_','2025-02-16_3_','2025-02-17_3_','2025-02-18_3_','2025-02-19_3_']   # es solapa amb ventilador no funcionant. i els dies 20 i 21 no hi son
datafiles['30'][4]=['2025-05-08_3_','2025-05-09_3_','2025-05-10_3_','2025-05-11_3_','2025-05-12_3_']


# Refrigerator  2  
datafiles['32'][0]=['2025-03-13_3_','2025-03-14_3_','2025-03-15_3_','2025-03-16_3_']
datafiles['32'][1]=['2025-03-10_3_']
datafiles['32'][2]=['2025-03-17_3_']
datafiles['32'][3]=['2025-03-24_3_','2025-03-25_3_','2025-03-26_3_']
datafiles['32'][4]=['2025-03-27_3_','2025-03-28_3_']


# Refrigerator  4  
datafiles['34'][0]=['2025-03-31_3_','2025-04-01_3_','2025-04-02_3_','2025-04-03_3_']
datafiles['34'][1]=['2025-04-25_3_','2025-04-26_3_','2025-04-27_3_','2025-04-28_3_'] # aquestes dades no hi son
#  diuen: que amb aquest set no obtenen gel i aquest experiment es pot descartar ????
datafiles['34'][2]=['2025-04-11_3_','2025-04-12_3_','2025-04-13_3_','2025-04-14_3_']
datafiles['34'][3]=['2025-04-30_3_','2025-05-01_3_','2025-05-02_3_','2025-05-03_3_','2025-05-04_3_','2025-05-05_3_']
datafiles['34'][4]=['2025-04-23_3_','2025-04-24_3_','2025-04-25_3_']



# data files arrays. Index:
# 0.  No failure
# 1.  Blocked evaporator
# 2.   Full Blocked condenser
# 3   Fan condenser not working
# 4.  Open door


facilitySetpoint='5-18'
# Features suggested by Xavier
# Care with 'tc s3' because on datafiles[0] is always nulll
# Seems to be incoropored in new tests

#r1s5 supply air flow temperature
#r1s1 inlet evaporator temperature
#r1s4 condenser outlet

# VAriables r1s4 and pa1 apiii  may not exists in cloud controlers


features=['r1 s1','r1 s4','r1 s5','pa1 apiii']
#features=['r1 s1','r1 s4','r1 s5']
featureNames={}
featureNames['r1 s1']='$T_{evap}$'
featureNames['r1 s4']='$T_{cond}$'
featureNames['r1 s5']='$T_{air}$'
featureNames['pa1 apiii']='$P_{elec}$'

unitNames={}
unitNames['r1 s1']='$(^{o}C)$'
unitNames['r1 s4']='$(^{o}C)$'
unitNames['r1 s5']='$(^{o}C)$'
unitNames['pa1 apiii']='$(W)$'


#features=['r1 s1','r1 s2','r1 s3','r1 s4','r1 s5','r1 s6','r1 s7','r1 s8','r1 s9','r1 s10','r2 s1','r2 s2','r2 s3','r2 s4','r2 s5','r2 s6','r2 s7','r2 s8','r2 s9','pa1 apiii','tc s1','tc s2']

#features=['r2 s2', 'tc s1','r1 s10','r1 s6','r2 s8']

NumFeatures=len(features)

df_list=[]
for i in range(NumberOfFailures+1):
    df_list.append([])

for i in range(NumberOfFailures+1):
    dftemp=[]
    for f in datafiles[facilitySetpoint][i]:
        print("                 ", f)
        #df1 = pd.read_csv('./data/'+f+'.csv', parse_dates=['datetime'], dayfirst=True, index_col='datetime')
        df1 = pd.read_csv('./data/'+f+'.csv')
        dftemp.append(df1)
    df_list[i]=pd.concat(dftemp)


# subsampled to 5'  =  30 * 10"
# We consider smaples every 5' because in production, we will only have data at this frequency
subsamplingrate=30

dataframe=[]
for i in range(NumberOfFailures+1):
    dataframe.append([])

for i in range(NumberOfFailures+1):
    datalength=df_list[i].shape[0]
    dataframe[i]=df_list[i].iloc[range(0,datalength,subsamplingrate)][features]
    dataframe[i].reset_index(inplace=True,drop=True)
    dataframe[i].dropna(inplace=True)


# Train data is first 2/3 of data
# Test data is: last 1/3 of data 
dataTrain=[]
dataTest=[]
NumberOfSamplesForTest=0

for i in range(NumberOfFailures+1):
    dataTrain.append(dataframe[i].values[0:int(dataframe[i].shape[0]*2/3),:])
    if NumberOfSamplesForTest==0:  # Take all
        dataTest.append(dataframe[i].values[int(dataframe[i].shape[0]*2/3):,:])
    else:
        dataTest.append(dataframe[i].values[int(dataframe[i].shape[0]*2/3):int(dataframe[i].shape[0]*2/3)+NumberOfSamplesForTest,:])

# Calculate means and stdev
a=dataTrain[0]
for i in range(1,NumberOfFailures+1):
    a=np.vstack((a,dataTrain[i]))

means=a.mean(axis=0) 
stdevs=a.std(axis=0)
def normalize2(train,test):
    return( (train-means)/stdevs, (test-means)/stdevs )

dataTrainNorm=[]
dataTestNorm=[]
for i in range(NumberOfFailures+1):
    dataTrainNorm.append([])
    dataTestNorm.append([])

for i in range(NumberOfFailures+1):
    (dataTrainNorm[i],dataTestNorm[i])=normalize2(dataTrain[i],dataTest[i])

def plotData():    
    fig, axes = plt.subplots(
        nrows=NumberOfFailures+1, ncols=NumFeatures, figsize=(15, 20), dpi=80, facecolor="w", edgecolor="k",sharex=True
    )
    for i in range(NumberOfFailures+1):
        for  j in range(NumFeatures):
            axes[i][j].plot(np.concatenate((dataTrainNorm[i][:,j],dataTestNorm[i][:,j])),label="Fail "+str(i)+",  feature 0")
    #axes[1].legend()
    #axes[0].set_ylabel(features[0])
    #axes[1].set_ylabel(features[1])
    plt.show()

if options.plot:
    plotData()
    exit()

def create_sequences(values, time_steps):
    output = []
    for i in range(len(values) - time_steps ):
        output.append(values[i : (i + time_steps)])
    return np.stack(output)

def listToString(l):
    r=''
    for i in l:
        r+=str(i)
    return(r.replace(' ',''))

timesteps=int(options.timesteps)

X=dataTrainNorm[0]
for i in range(1,NumberOfFailures+1):
    X=np.vstack((X,dataTrainNorm[i]))

xtrain=create_sequences(X,timesteps)

km = TimeSeriesKMeans(n_clusters=NumberOfFailures+1, metric="dtw", random_state=0)
modelpath="model_v2_unsupervised_"+facilitySetpoint+"_"+str(timesteps)+listToString(features)+".pk"
if options.train:
    km.fit(xtrain)
    km.to_pickle(modelpath)
else:
    km.from_pickle(modelpath)
    #km.fit_predict(xtrain)

colorline=['violet','lightcoral','cyan','lime','grey']
colordot=['darkviolet','red','blue','green','black']


featuresToPlot=features
indexesToPlot=[]
for i in featuresToPlot:
    indexesToPlot.append(features.index(i))


def plot(data,ranges):
    km.fit_predict(data)
# Expand data to plot with the timesteps samples of last sample
    datatoplot=data[:,0,:]
    datatoplot=np.vstack((datatoplot,data[ranges[len(ranges)-1][1],:,:]))

    labels=[]    #   Labels are assigned randomly by classifer
    for i in range(len(ranges)):
        b=Counter(km.labels_[ranges[i][0]:ranges[i][1]])
        labels.append(b.most_common(1)[0][0])

    print("\n\n\n LABELS: ",labels,"\n\n")

    NumFeaturesToPlot=len(indexesToPlot)
    plt.rcParams.update({'font.size': 16})
    fig, axes = plt.subplots(
        nrows=NumFeaturesToPlot, ncols=1, figsize=(15, 10), dpi=80, facecolor="w", edgecolor="k",sharex=True
    )
    for i in range(NumFeaturesToPlot):
        init=0
        end=ranges[0][1]
        labelsplotted=[]
        for j in range(len(ranges)):
            if j==(len(ranges)-1): # Plot the last timesteps
                classtype=labels.index(labels[j])
                if classtype in labelsplotted:
                    axes[i].plot(range(init,end+timesteps),datatoplot[ranges[j][0]:ranges[j][1]+timesteps,indexesToPlot[i]]*stdevs[i]+means[i], color=colorline[classtype],linewidth=1)
                else:
                    axes[i].plot(range(init,end+timesteps),datatoplot[ranges[j][0]:ranges[j][1]+timesteps,indexesToPlot[i]]*stdevs[i]+means[i], label="Class: "+str(classtype), color=colorline[classtype],linewidth=1)
                    labelsplotted.append(classtype)
            else:
                classtype=labels.index(labels[j])
                if classtype in labelsplotted:
                    axes[i].plot(range(init,end),datatoplot[ranges[j][0]:ranges[j][1],indexesToPlot[i]]*stdevs[i]+means[i], color=colorline[classtype],linewidth=1)
                else:
                    axes[i].plot(range(init,end),datatoplot[ranges[j][0]:ranges[j][1],indexesToPlot[i]]*stdevs[i]+means[i], label="Class: "+str(classtype), color=colorline[classtype],linewidth=1)
                    labelsplotted.append(classtype)
            init=end
            if j<(len(ranges)-1):
                end+=(ranges[j+1][1]-ranges[j+1][0])
        x=[]
        y=[]
        for j in range(len(ranges)):
            x.append([])
            y.append([])

        for j in range(len(ranges)):
            for k in range(ranges[j][0],ranges[j][1]):
                try:  # Idont know why sometimes fails index !!!!
                    x[labels.index(km.labels_[k])].append(k+timesteps)
                    y[labels.index(km.labels_[k])].append(datatoplot[k+timesteps,indexesToPlot[i]]*stdevs[i]+means[i])
                except:
                    x[0].append(k+timesteps)
                    y[0].append(datatoplot[k+timesteps,indexesToPlot[i]]*stdevs[i]+means[i])

        labelsplotted=[]
        for j in range(len(ranges)):
            classtype=labels.index(labels[j])
            if classtype in labelsplotted:
                axes[i].plot(x[j],y[j] ,color=colordot[labels.index(labels[j])],marker='.',linewidth=0)
            else:
                axes[i].plot(x[j],y[j] ,color=colordot[labels.index(labels[j])],marker='.',linewidth=0,label="Class type "+str(j) )
                labelsplotted.append(classtype)

        if i==(NumFeatures-1):
            axes[i].legend(loc='right')
        s=''
        s+=featureNames[features[indexesToPlot[i]]]
        s+=' '+unitNames[features[indexesToPlot[i]]]
        axes[i].set_ylabel(s)
        axes[i].grid()

    axes[NumFeaturesToPlot-1].set_xlabel("Sample number")
    plt.show()
        
   
'''
Ranges=[]
r=0
for i in range(NumberOfFailures+1):
    Ranges.append([r,r+dataTrainNorm[i].shape[0]])
    r+=dataTrainNorm[i].shape[0]
# Drop the last TIME_STEPS for plotting
Ranges[NumberOfFailures][1]=Ranges[NumberOfFailures][1]-timesteps-1

plot(xtrain,Ranges)
'''
# Try with test data
X=dataTestNorm[0]
Ranges=[[0,dataTestNorm[0].shape[0]]]
r=dataTestNorm[0].shape[0]
for i in range(1,NumberOfFailures+1):
    X=np.vstack((X,dataTestNorm[i]))
    Ranges.append([r,r+dataTestNorm[i].shape[0]])
    r+=dataTestNorm[i].shape[0]

#X=np.vstack((X,dataTestNorm[0]))  # We add a last segment of no fail data
#Ranges.append([r,r+dataTestNorm[0].shape[0]])

Ranges[len(Ranges)-1][1]=Ranges[len(Ranges)-1][1]-timesteps-1

xtest=create_sequences(X,timesteps)




def anomalyMetric(labels,ranges):
    # Takes ONLY the firsts segments of Ranges
    # FP, TP: false/true positive
    # TN, FN: true/false negative
    # Sensitivity (recall): probab failure detection if data is fail: TP/(TP+FN)
    # Precision: Rate of positive results:  TP/(TP+FP)  
    # F1-score: predictive performance measure: 2*Precision*Sensitity/(Precision+Sensitity)

    lab=[]    #   Labels are assigned randomly by classifer
    TP=[]
    FN=[]
    TPFP=[]
    for i in range(NumberOfFailures+1):
        TP.append([])
        FN.append([])
        TPFP.append([])
        b=Counter(labels[ranges[i][0]:ranges[i][1]])
        lab.append(b.most_common(1)[0][0])

    for i in range(NumberOfFailures+1):
        counttp=0
        countfn=0
        for j in range(ranges[i][0],ranges[i][1]-timesteps):
            if lab.index(labels[j])==i:
                counttp+=1
            else:
                countfn+=1
        TP[i]=counttp
        FN[i]=countfn

    for i in range(NumberOfFailures+1):
        count=0
        for ii in range(NumberOfFailures+1):
            for j in range(ranges[ii][0],ranges[ii][1]-timesteps):
                if lab.index(labels[j])==i:
                    count+=1
            TPFP[i]=count

    segmentLength=[]
    for i in range(NumberOfFailures+1):
        segmentLength.append(ranges[i][1]-timesteps-ranges[i][0])
    totalSegmentLength=0
    for i in range(NumberOfFailures+1):
        totalSegmentLength+=segmentLength[i]


    Sensitivity=0
    Precision=0
    for i in range(NumberOfFailures+1):
        Sensitivity+=TP[i]/(TP[i]+FN[i])*segmentLength[i]/totalSegmentLength
        Precision+=TP[i]/(TPFP[i])*segmentLength[i]/totalSegmentLength


    print(lab)
    print("TP: ",TP)
    print("FN: ",FN)
    print("TPFP: ",TPFP)
    print("Sensitivity: ",Sensitivity)
    print("Precision: ",Precision)
    print("F1-Score: ",2*Precision*Sensitivity/(Sensitivity+Precision))


km.fit_predict(xtest)
anomalyMetric(km.labels_,Ranges)
plot(xtest,Ranges)