# %% Load Data Set and denoise with moving avarage
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
import numpy as np
import tikzplotlib
import time
import matplotlib.pyplot as plt

from matplotlib import rc

# Aktiviert LaTeX in Matplotlib
rc('text', usetex=True)
rc('font', **{'family': 'serif', 'serif': ['Computer Modern Roman'], 'size': 10})

# 1. Load the dataset
data = pd.read_csv('Data310823_merge.csv')
#data = pd.read_csv('simulated_machine_data_extralong.csv')


# %% ####### PLOT Only orginal power values 

import matplotlib.ticker as ticker


Powerinkw = data['Pges'] / 1000 
plt.figure(figsize=(10, 6))
plt.plot(Powerinkw, color='black')
plt.xlabel('Time in s', fontsize=16)
plt.ylabel('Power in kW ',fontsize=16)

plt.grid(True, which='major', linestyle='-', linewidth=0.5, color='grey', alpha=0.5)  # Hauptgitterlinien

plt.gca().xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16)  # Schriftgröße für Haupt-Ticks
ax.tick_params(axis='both', which='minor', labelsize=10)  # Schriftgröße für Neben-Ticks

#-->plt.savefig('emag_power_orginal.png')
#plt.savefig('emag_power_orginal.svg', format='svg')
plt.show()
# %%
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<


# Minimum and Maximum values
min_value = data['Pges'].min()
max_value = data['Pges'].max()

# Overall Mean power consumption
mean_value = data['Pges'].mean()

# Detection of outliers using the IQR method
Q1 = data['Pges'].quantile(0.25)
Q3 = data['Pges'].quantile(0.75)
IQR = Q3 - Q1
lower_bound = Q1 - 1.5 * IQR
upper_bound = Q3 + 1.5 * IQR
outliers = data[(data['Pges'] < lower_bound) | (data['Pges'] > upper_bound)]

non_outliers = data[(data['Pges'] >= lower_bound) & (data['Pges'] <= upper_bound)]

# Print the results
print(f"Minimum Value: {min_value}")
print(f"Maximum Value: {max_value}")
print(f"Mean Value: {mean_value}")
# print(f"Number of Outliers: {outliers}")

# 1. Denoise using Moving Average (Simple Low Pass Filter)
window_size = 50 ## You can adjust this based on your preference
denoised_signal = data['Pges'].rolling(window=window_size).mean()

# 2. Compute Moving Average (this is similar to the denoised signal but can be computed with a different window size if desired)
moving_avg = non_outliers['Pges'].rolling(window=window_size).mean()   


#Visualization
plt.figure(figsize=(10, 6))
plt.xlabel('Time in s', fontsize=16)
plt.ylabel('Power in kW', fontsize=16)
plt.plot(data['Pges']/1000, label='Original Signal',color='gray',linestyle='dashed', alpha=0.3)
plt.plot(denoised_signal/1000, label='Denoised Signal', color='gray')
plt.plot(moving_avg/1000, label='Moving Average', color='black')
#plt.title('Original, Denoised, and Moving Average Signals')
# Place the legend below the plot
plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=3, fontsize=16)
plt.tight_layout()
plt.grid(True, which='major', linestyle='-', linewidth=0.5, color='grey', alpha=0.5)  # Hauptgitterlinien
# Ticker-Formatter für die X-Achse ohne Tausendertrennzeichen
plt.gca().xaxis.set_major_formatter(ticker.StrMethodFormatter('{x:.0f}'))
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16)  # Schriftgröße für Haupt-Ticks
ax.tick_params(axis='both', which='minor', labelsize=10)  # Schriftgröße für Neben-Ticks

#plt.savefig('emag_power_average.png')
#plt.savefig('emag_power_average.svg', format='svg')
plt.show()


# %%  
# ################>>>>>>>MODEL DEPLOYMENT with BATCH K MEANS Cluster for Power state idtification

import time 

start_time_batch = time.time()

# 1. Load the dataset

moving_avg_cleaned = moving_avg.dropna().values.reshape(-1, 1)
scaler = StandardScaler()
moving_avg_scaled = scaler.fit_transform(moving_avg_cleaned)
# 2. Apply KMeans clustering
kmeans = KMeans(n_clusters=3, random_state=42)
clusters = kmeans.fit_predict(moving_avg_scaled)

# 3. Filter clusters based on center values and size
centers = scaler.inverse_transform(kmeans.cluster_centers_)

valid_clusters = []
for idx, center in enumerate(centers):
    if center[0] > 0 and list(clusters).count(idx) > 100:
        valid_clusters.append((idx, center[0]))

# 4. For each valid cluster, determine the boundaries
boundaries = []
for idx, center in valid_clusters:
    lower_boundary = center - 0.10 * center
    upper_boundary = center + 0.10 * center
    boundaries.append((lower_boundary, upper_boundary))

# Print the valid cluster centers and their boundaries
for (idx, center), (lower, upper) in zip(valid_clusters, boundaries):
    print(f"Cluster {idx}: Center at value {center:.2f} with boundaries [{lower:.2f}, {upper:.2f}]")

import csv   
# Create the data table
clustercenters = [["Powerstate", "Value"]]

# Populate the table with cluster center values
for idx, (cluster_idx, center) in enumerate(valid_clusters):
    clustercenters.append([f"Cluster {idx+1}", "{:.2f}".format(center)])


end_time_batch = time.time()


# # Save the data table to a CSV file
# with open('clustercenters.csv', 'w', newline='') as f:
#     writer = csv.writer(f)
#     writer.writerows(clustercenters) 


print('>>>>>>>Calculate Clusters Successfull<<<<<<')



print('>>>>>>>Plot Clusters <<<<<<')
# Define a color palette for clusters
cluster_colors = ['lightgray', 'gray', 'darkgray','blue','yellow']

# Plot the original power values
fig=plt.figure(figsize=(10, 6))
plt.plot(moving_avg, label='Moving Average Signal ', color='black',alpha=0.6)

# Plot the cluster centers and their boundaries with specified colors
for idx, ((cluster_idx, center), (lower, upper)) in enumerate(zip(valid_clusters, boundaries)):
    color = cluster_colors[idx]
    plt.axhline(center, color=color, alpha=0.7, label=f'Cluster {cluster_idx+1} Center')
    plt.axhline(lower, color=color, linestyle='--', alpha=0.5, label=f'Cluster {cluster_idx+1} Lower Boundary')
    plt.axhline(upper, color=color, linestyle='--', alpha=0.5, label=f'Cluster {cluster_idx+1} Upper Boundary')


def tikzplotlib_fix_ncols(obj):
    """
    workaround for matplotlib 3.6 renamed legend's _ncol to _ncols, which breaks tikzplotlib
    """
    if hasattr(obj, "_ncols"):
        obj._ncol = obj._ncols
    for child in obj.get_children():
        tikzplotlib_fix_ncols(child)


# Customize the plot
#plt.title('Power Values with KMeans Clustering Centers and Boundaries')
plt.xlabel('Data Samples',fontsize=16)
plt.ylabel('Power in Watt',fontsize=16)
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16,)  # Schriftgröße für Haupt-Ticks
# Place the legend below the plot
plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=3,fontsize=16)
plt.tight_layout()
tikzplotlib_fix_ncols(fig)

#plt.savefig('emag_batchkmeans_short.png')

#plt.savefig('emag_batchkmeans_short.svg',format='svg')
#plt.savefig('emag_batchkmeans_long.svg',format='svg')
plt.show()
# plt.show(block=False)  # Display the figure without blocking the script
# time.sleep(5)          # Wait for 5 seconds
# plt.close()            # Close the figure
# %%
from sklearn.metrics import silhouette_score
from sklearn.metrics import davies_bouldin_score

###### ##################################Bewertung Güte und Geschwindigkeit
#Calculate the Silhouette Score
silhouette_avg = silhouette_score(moving_avg_scaled, clusters)
print(f"Silhouette Score: {silhouette_avg}")

# Calculate the Davis-Bouldin Score
davies_bouldin_avg = davies_bouldin_score(moving_avg_scaled, clusters)
print(f"Davis-Bouldin Score: {davies_bouldin_avg}")

duration_batch_kmeans = end_time_batch - start_time_batch
print(f"Duration Batch K-Means {duration_batch_kmeans}")

# Anzahl der Samples bestimmen
num_samples = data.shape[0]

print(f"Anzahl der Samples: {num_samples}")

print('>>>>>>>FINISH Batch K-Means  Clustering Calculation Successfull<<<<<<')

# %%
# ################ PLOT BOX PLOT


# Organisiere die Datenpunkte nach Clustern
clustered_data = [[] for _ in range(len(centers))]
for datapoint, cluster_label in zip(moving_avg_cleaned, clusters):
    clustered_data[cluster_label].append(datapoint)
# Konvertiere jedes Cluster-Element in ein 1D-Array, falls notwendig
for idx, cluster in enumerate(clustered_data):
    cluster_array = np.array(cluster).flatten()  # Flattening des Arrays, falls es mehrdimensional ist
    clustered_data[idx] = cluster_array

# Stellen Sie sicher, dass alle Cluster jetzt 1D sind
for idx, cluster in enumerate(clustered_data):
    if np.array(cluster).ndim != 1:
        raise ValueError(f"Cluster an Index {idx} hat eine Dimension {np.array(cluster).ndim}, erwartet wird 1.")

# Erstelle den Boxplot
plt.figure(figsize=(10, 6))
bp = plt.boxplot(clustered_data, showfliers=False, patch_artist=False)


# Beschriftungen und Titel
plt.xlabel('Cluster')
plt.ylabel('Power Value')
plt.title('Boxplot of Power Values by Cluster')

# Setze individuelle Clusterbezeichnungen
cluster_labels = ['Standby', 'Disabled', 'Working']
plt.xticks(ticks=range(1, len(centers) + 1), labels=cluster_labels)
# Zeichne die Clusterzentren
for i, center in enumerate(centers):
    plt.plot([i+1], [center], 'ko', markersize=10)

# Setze die Y-Achsen-Skalierung in 250er Schritten
plt.yticks(np.arange(0, max([item for sublist in clustered_data for item in sublist]), 250))
# Zeige ein Hilfsgitter in leichtem Grau für die Y-Achse an
plt.grid(axis='y', color='0.95')

plt.legend()

# Patchen des tikzplotlib, um das Problem mit der Legenden-Attribut zu umgehen
def patch_legend(legend):
    if hasattr(legend, '_ncols'):
        legend._ncol = legend._ncols
    else:
        legend._ncols = legend._ncol
    for handle in legend.legendHandles:
        handle.set_sizes([64.0])  # Sie können die Größe nach Bedarf anpassen

patch_legend(plt.legend())


plt.show()
# Berechnen der statistischen Werte für jeden Cluster
# Berechnen der statistischen Werte für jeden Cluster und Speichern in CSV
with open('cluster_value_measure.csv', 'w', newline='') as file:
    writer = csv.writer(file)
    writer.writerow(['Cluster','Center','Min', 'Q1', 'Median', 'Q3', 'Max'])

    for i, cluster_data in enumerate(clustered_data):
        cluster_data_np = np.array(cluster_data)
        
        center = centers[i]
        # bound_cluster_up = boundaries[i][1]
        # bound_cluster_low = boundaries[i][0]
        Q1 = np.percentile(cluster_data_np, 25)
        Q2 = np.median(cluster_data_np)
        Q3 = np.percentile(cluster_data_np, 75)
        IQR = Q3 - Q1
        min_val = np.min(cluster_data_np)
        max_val = np.max(cluster_data_np)
        # lower_bound = Q1 - 1.5 * IQR
        # upper_bound = Q3 + 1.5 * IQR
        
        # Schreiben der berechneten Werte in die CSV-Datei
        writer.writerow([f'{cluster_labels[i]}', "{:.2f}".format(center[0]),"{:.2f}".format(min_val), "{:.2f}".format(Q1), 
                         "{:.2f}".format(Q2), "{:.2f}".format(Q3), "{:.2f}".format(max_val)])
                        #  ,"{:.2f}".format(lower_bound), "{:.2f}".format(upper_bound)])
        # Ausgabe der berechneten Werte
        print(f"Cluster {cluster_labels[i]}:")
        print(f"Center: {center}")
        # print(f"ClusterBoundUp: {bound_cluster_up}")
        # print(f"ClusterBoundLow: {bound_cluster_low}")        
        print(f"Minimum: {min_val}")
        print(f"Q1 (unteres Quartil): {Q1}")
        print(f"Median (Q2): {Q2}")
        print(f"Q3 (oberes Quartil): {Q3}")
        print(f"Maximum: {max_val}")
        # print(f"Ausreißer unterhalb: {lower_bound}")
        # print(f"Ausreißer oberhalb: {upper_bound}")
        print("------")

# Lade die CSV-Daten
df = pd.read_csv('cluster_value_measure.csv')

# Transponiere die Daten
df_transposed = df.T

# Speichere die transponierten Daten als neue CSV-Datei
df_transposed.to_csv('cluster_value_measure_transposed.csv', header=False)



##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Gütemaße
#from sklearn import metrics
from sklearn.metrics import silhouette_score
from sklearn.metrics import davies_bouldin_score


# Calculate the Silhouette Score
silhouette_avg = silhouette_score(moving_avg_scaled, clusters)
print(f"Silhouette Score: {silhouette_avg}")

# Calculate the Davis-Bouldin Score
davies_bouldin_avg = davies_bouldin_score(moving_avg_scaled, clusters)
print(f"Davis-Bouldin Score: {davies_bouldin_avg}")


print('>>>>>>>FINISH Batch Clustering Calculation Successfull<<<<<<')




##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %% 
# ONLINE CLUSTERING MINI BATCH K MEANS -sklearn - short data set
import pandas as pd
from sklearn.cluster import MiniBatchKMeans
from sklearn.preprocessing import StandardScaler
import numpy as np
import matplotlib.pyplot as plt
import time 
# Load the dataset

data_online = pd.read_csv('Data310823_merge.csv')

# Preprocess the data using the previous steps and then apply MiniBatchKMeans clustering
start_time_minibatch = time.time()
# Denoise using Moving Average
window_size = 50
moving_avg_online = data_online['Pges'].rolling(window=window_size).mean()

# Preprocess the moving average data
moving_avg_cleaned_online = moving_avg_online.dropna().values.reshape(-1, 1)
scaler_online = StandardScaler()
moving_avg_scaled_online = scaler_online.fit_transform(moving_avg_cleaned_online)

# Initialize MiniBatchKMeans and fit
mbkmeans_online = MiniBatchKMeans(n_clusters=3, random_state=42, batch_size=100)
for i in range(0, len(moving_avg_scaled_online), 100):
    batch = moving_avg_scaled_online[i:i+100]
    mbkmeans_online.partial_fit(batch)

# Predict clusters for all data
predicted_clusters_online = mbkmeans_online.predict(moving_avg_scaled_online)

# Get the cluster centers
cluster_centers_online = scaler_online.inverse_transform(mbkmeans_online.cluster_centers_)
print(f"Online KMeans Mini Batch Cluster Center: {cluster_centers_online}")

end_time_minibatch = time.time()

# Plot the results
# Grautöne für Cluster
cluster_colors = ['lightgray', 'gray', 'darkgray']

fig=plt.figure(figsize=(10, 6))
plt.plot(moving_avg_online, label='Moving Average', color='black',alpha=0.6)
for idx, center in enumerate(cluster_centers_online):
    plt.axhline(center[0], label=f'Cluster {idx+1} Center', color=cluster_colors[idx], linestyle='--')


#plt.title('Moving Average with MiniBatchKMeans Cluster Centers')
plt.xlabel('Data Samples', fontsize=16)
plt.ylabel('Power in Watt', fontsize=16)
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16)  # Schriftgröße für Haupt-Ticks

plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=2,fontsize=16)

plt.tight_layout()
plt.savefig('emag_minibatchkmeans_short.svg', format='svg')
plt.show()


from sklearn.metrics import silhouette_score
from sklearn.metrics import davies_bouldin_score

# Calculate the Silhouette Score
silhouette_avg_online = silhouette_score(moving_avg_scaled_online, predicted_clusters_online)
print(f"Silhouette Score for Online Mini Batch KMeans: {silhouette_avg_online}")

# Calculate the Davis-Bouldin Score
davies_bouldin_index_online = davies_bouldin_score(moving_avg_scaled_online, predicted_clusters_online)
print(f"Davies-Bouldin Index for Online Mini Batch KMeans: {davies_bouldin_index_online}")


duration_online_mini_batch_kmeans = end_time_minibatch - start_time_minibatch
print(f"Duration Batch K-Means {duration_online_mini_batch_kmeans}")

# Anzahl der Samples bestimmen
num_samples_online = data_online.shape[0]

print(f"Anzahl der Samples: {num_samples_online}")


print('>>>>>>>FINISH Online Clustering  MINI BATCH K MEANS Calculation Successfull<<<<<<')
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %% 
# ONLINE CLUSTERING MINI BATCH K MEANS -sklearn - long data set
import pandas as pd
from sklearn.cluster import MiniBatchKMeans
from sklearn.preprocessing import StandardScaler
import numpy as np
import matplotlib.pyplot as plt
import time 
# Load the dataset

#data_online = pd.read_csv('Data310823_merge.csv')
data_online = pd.read_csv('simulated_machine_data_extralong.csv')

# Preprocess the data using the previous steps and then apply MiniBatchKMeans clustering
start_time_minibatch = time.time()
# Denoise using Moving Average
window_size = 50
moving_avg_online = data_online['Pges'].rolling(window=window_size).mean()

# Preprocess the moving average data
moving_avg_cleaned_online = moving_avg_online.dropna().values.reshape(-1, 1)
scaler_online = StandardScaler()
moving_avg_scaled_online = scaler_online.fit_transform(moving_avg_cleaned_online)

# Initialize MiniBatchKMeans and fit
mbkmeans_online = MiniBatchKMeans(n_clusters=3, random_state=42, batch_size=100)
for i in range(0, len(moving_avg_scaled_online), 100):
    batch = moving_avg_scaled_online[i:i+100]
    mbkmeans_online.partial_fit(batch)

# Predict clusters for all data
predicted_clusters_online = mbkmeans_online.predict(moving_avg_scaled_online)

# Get the cluster centers
cluster_centers_online = scaler_online.inverse_transform(mbkmeans_online.cluster_centers_)
print(f"Online KMeans Mini Batch Cluster Center: {cluster_centers_online}")

end_time_minibatch = time.time()

# Plot the results

# Grautöne für Cluster
cluster_colors = ['lightgray', 'gray', 'darkgray']

fig=plt.figure(figsize=(10, 6))
plt.plot(moving_avg_online, label='Moving Average', color='black',alpha=0.6)
for idx, center in enumerate(cluster_centers_online):
    plt.axhline(center[0], label=f'Cluster {idx+1} Center', color=cluster_colors[idx], linestyle='--')


#plt.title('Moving Average with MiniBatchKMeans Cluster Centers')
plt.xlabel('Data Samples', fontsize=16)
plt.ylabel('Power in Watt', fontsize=16)
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16)  # Schriftgröße für Haupt-Ticks

plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=2,fontsize=16)

plt.tight_layout()
plt.savefig('emag_minibatchkmeans_long.svg', format='svg')
plt.show()


from sklearn.metrics import silhouette_score
from sklearn.metrics import davies_bouldin_score

# Calculate the Silhouette Score
silhouette_avg_online = silhouette_score(moving_avg_scaled_online, predicted_clusters_online)
print(f"Silhouette Score for Online Mini Batch KMeans: {silhouette_avg_online}")

# Calculate the Davis-Bouldin Score
davies_bouldin_index_online = davies_bouldin_score(moving_avg_scaled_online, predicted_clusters_online)
print(f"Davies-Bouldin Index for Online Mini Batch KMeans: {davies_bouldin_index_online}")


duration_online_mini_batch_kmeans = end_time_minibatch - start_time_minibatch
print(f"Duration Batch K-Means {duration_online_mini_batch_kmeans}")

# Anzahl der Samples bestimmen
num_samples_online = data_online.shape[0]

print(f"Anzahl der Samples: {num_samples_online}")


print('>>>>>>>FINISH Online Clustering  MINI BATCH K MEANS Calculation Successfull<<<<<<')
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<




# %% 
# # ONLINE CLUSTERING with  River CLUSTREAM   
from river import cluster
from river import stream
from river import utils
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import time 
# Load the dataset
#data_online = pd.read_csv('Data310823_merge.csv')
data_online = pd.read_csv('simulated_machine_data_extralong.csv')


start_time_online_Clustream  = time.time()

# Denoise using Moving Average
window_size = 50
moving_avg_online = data_online['Pges'].rolling(window=window_size).mean().dropna()

# Initialize CluStream
clustream = cluster.CluStream()


# Online CluStream update
for i, val in enumerate(moving_avg_online):
    clustream = clustream.learn_one({'value': val})

# Predict clusters for all data (we simulate this process as CluStream does not provide direct predict method)
# Instead, we can calculate the distance of each point to the cluster centroids and assign it to the nearest centroid.
def euclidean_distance(point, centroid):
    return utils.math.minkowski_distance(point, centroid, 2)

def predict_clusters(stream, model):
    centroids = model.centers
    labels = []
    for x in stream:
        # Berechne den euklidischen Abstand zum nächsten Zentrum
        nearest_centroid, min_distance = None, float('inf')
        for c_id, c_info in centroids.items():
            dist = euclidean_distance(x, c_info)
            if dist < min_distance:
                nearest_centroid, min_distance = c_id, dist
        labels.append(nearest_centroid)
    return labels

# Simulate streaming data
streaming_data = [{'value': val} for val in moving_avg_online]

# Use the predict_clusters function to assign clusters
predicted_clusters = predict_clusters(streaming_data, clustream)



# Convert to NumPy array for metric calculations
values = moving_avg_online.values.reshape(-1, 1)
predicted_clusters = np.array(predicted_clusters)

end_time_online_Clustream = time.time()

# Plot the results


# Plot the results
plt.figure(figsize=(15, 6))
plt.plot(moving_avg_online.values, label='Moving Average', alpha=0.6)

# Plot the cluster centers
for center in clustream.centers.values():
    plt.axhline(center['value'], linestyle='--')
    
plt.title('Moving Average with CluStream Cluster Centers')
plt.xlabel('Data Samples (Ordered by Timestamp)')
plt.ylabel('Power Value in Watt')
plt.legend()
plt.show()

# %%
# Grautöne für Cluster
cluster_colors = ['lightgray', 'gray', 'darkgray']
displayed_indices = [0, 1, 2]  # Hier können Sie die Cluster angeben, die in der Legende erscheinen sollen

fig=plt.figure(figsize=(10, 6))
plt.plot(moving_avg_online, label='Moving Average', color='black',alpha=0.6)
# Plot der Clusterzentren in unterschiedlichen Grautönen (nur bis zu 3 Cluster)
for idx, (key, center) in enumerate(clustream.centers.items()):
    if idx in displayed_indices:  # nur bis zur Länge von cluster_colors
        plt.axhline(center['value'], label=f'Cluster {idx+1} Center', color=cluster_colors[idx], linestyle='--')


#plt.title('Moving Average with MiniBatchKMeans Cluster Centers')
plt.xlabel('Data Samples', fontsize=16)
plt.ylabel('Power in Watt', fontsize=16)
ax = plt.gca()  # Greift auf das aktuelle Achsen-Objekt zu
ax.tick_params(axis='both', which='major', labelsize=16)  # Schriftgröße für Haupt-Ticks

plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15), ncol=2,fontsize=16)

plt.tight_layout()
#plt.savefig('emag_clustream_short.svg', format='svg')
plt.savefig('emag_clustream_long.svg', format='svg')
plt.show()


# Plot the cluster centers und Drucke die Zentren
for center_id, center in enumerate(clustream.centers.values()):
    print(f"Cluster {center_id}: Center at value {center['value']:.2f}")
# %%
# #print(f"Cluster 0: Center at value {center:.2f} 

# # Calculate Silhouette Score and Davies-Bouldin Index
# silhouette_avg = silhouette_score(values, predicted_clusters)
# davies_bouldin = davies_bouldin_score(values, predicted_clusters)

# print(f"Silhouette Score for CluStream: {silhouette_avg}")
# print(f"Davies-Bouldin Index for CluStream: {davies_bouldin}")

# duration_online_Clustream = end_time_online_Clustream  - start_time_online_Clustream 
# print(f"Duration ClusStream Clustering {duration_online_Clustream }")

# # Anzahl der Samples bestimmen
# num_samples_online_Clustream = data_online.shape[0]

# print(f"Anzahl der Samples: {num_samples_online_Clustream}")

print('>>>>>>>FINISH Online Clustering CluStream Calculation Successfull<<<<<<')



##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %% 
# # BESTIMME POWER STATE INTERVALL
def find_power_states_intervals(data, column_name):
    """
    Finds the start and end indices for each power state in a given mode column.
    
    Args:
    - data (DataFrame): The dataset.
    - column_name (str): The name of the mode column.
    
    Returns:
    - List[Tuple]: A list of tuples, where each tuple contains the start and end index of a power state.
    """
    intervals = []
    start = None
    current_state = None
    
    for i, state in enumerate(data[column_name]):
        if start is None:
            start = i
            current_state = state
        elif current_state != state:
            intervals.append((start, i-1, current_state))
            start = i
            current_state = state
    
    # Handle the last interval
    if start is not None:
        intervals.append((start, i, current_state))
    
    return intervals

# Find intervals for each mode


modes = ['Mode_Working', 'Mode_Operation', 'Mode_Standby', 'Mode_Disabled']
intervals_dict = {}

for mode in modes:
    intervals_dict[mode] = find_power_states_intervals(data, mode)

# Print the intervals
for mode, intervals in intervals_dict.items():
    print(f"{mode} intervals: {intervals}")


##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %%  
# # BESTIMME HOLING DURATIONS
# Compute holding duration min max for each power state
def compute_durations(intervals):
    """
    Compute the duration for each interval.
    
    Args:
    - intervals (List[Tuple]): List of intervals.
    
    Returns:
    - List[Tuple]: A list of tuples, where each tuple contains the start, end index, power state value, and duration.
    """
    durations = [(start, end, value, end-start+1) for start, end, value in intervals]
    return durations

def find_min_max_durations(durations):
    """
    Finds the minimum and maximum durations for each power state.
    
    Args:
    - durations (List[Tuple]): List of tuples with durations.
    
    Returns:
    - Dict: Dictionary with power state values as keys and (min_duration, max_duration) as values.
    """
    unique_states = set([duration[2] for duration in durations])
    min_max_durations = {}
    for state in unique_states:
        state_durations = [duration[3] for duration in durations if duration[2] == state]
        min_max_durations[state] = (min(state_durations), max(state_durations))
    
    return min_max_durations

# Compute durations and find min and max durations for each mode
durations_dict = {}
min_max_durations_dict = {}

for mode in modes:
    durations = compute_durations(intervals_dict[mode])
    durations_dict[mode] = durations
    min_max_durations_dict[mode] = find_min_max_durations(durations)

# Print the min and max durations for each power state in each mode
for mode, min_max_durations in min_max_durations_dict.items():
    print(f"{mode} min-max durations:")
    for state, (min_duration, max_duration) in min_max_durations.items():
        print(f"  Power state {state}: Min duration = {min_duration}, Max duration = {max_duration}")

print('>>>>>>>Calculate Holding Duration for each power State Successfull<<<<<<')

##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %% 
#  # BESTIMME REACTION DURATION 

# Annahme, dass WZM GT immer mindestens im Standby betrieb (ist wenn Flexibilität gezogen)
# im Vergleich zum Referenzbetrieb. # Daher daher ist die Reaktionszeit die 0.5 sekunden
# minimum und maximal 90 Sekunden wenn gerade ein Bauteil gestartet 
# wurde  90 Sekunden weil Durchschnittlich Bauteilbearbeitungszeit = 90 sekunden ist. 
# In 15 Minuten Intervallen, können 10 Teile Produziert werden 
# 15 Minuten * 60 Sekunden / 10 Teile = 90 Sekunden
# Annhame: Es ist bekannt, das  in 15 Minuten maximal 10 Teile gefertigt werden können aus B & D Understanding 
import math
import latexify


cycle_time_quarter = 15*60 # sekunden
parts_per_quarter = 10

reaction_duration_max = cycle_time_quarter/parts_per_quarter

reaction_duration_min = 0.5 #sekunden


print(f" Reaction Duration min: {reaction_duration_min}, Reaction Duration max: {reaction_duration_max}")

print('>>>>>>>Calculate Reaction Duration Successfull<<<<<<')

#%%
@latexify.function
def calculate_reaction_duration(cycle_time_quarter, parts_per_quarter):
  reaction_duration_max = cycle_time_quarter / parts_per_quarter
  reaction_duration_min = 0.5
  return reaction_duration_max, reaction_duration_min

# Now, when you print the function, it should display the LaTeX representation.
print(calculate_reaction_duration(cycle_time_quarter,parts_per_quarter))
calculate_reaction_duration


# %% Calculate gradients
# CALCUALTE POWER GRADIENTS 
#data = pd.read_csv('Data310823_merge.csv')
def current_mode(row):
    """Determine the current mode based on the dataset's structure."""
    for col in ['Mode_Working', 'Mode_Operation', 'Mode_Standby', 'Mode_Disabled']:
        if row[col] == 1:
            return col
    return None

# Add a new column to the dataframe to represent the current mode for each row
data['Mode'] = data.apply(current_mode, axis=1)

def compute_gradients(data, from_state, to_state):
    """
    Compute the gradients for transitions between specified power states.
    
    Args:
    - data (DataFrame): The dataset.
    - from_state (str): The initial power state.
    - to_state (str): The next power state.
    
    Returns:
    - List[float]: A list of gradient values for the specified transitions.
    """
    gradients = []
    in_transition = False
    for i in range(1, len(data)):
        if data['Mode'].iloc[i-1] == from_state and data['Mode'].iloc[i] == to_state:
            if not in_transition:  # start of a new transition
                start_power = data['Pges'].iloc[i-1]
                start_time = i-1
                in_transition = True
            end_power = data['Pges'].iloc[i]
            end_time = i
        elif in_transition:  # end of the current transition
            gradient = (end_power - start_power) / (end_time - start_time)  # in W/second if time is in seconds
            gradients.append(gradient)
            in_transition = False
    return gradients

# Calculate gradients for specified transitions
transitions = [("Mode_Disabled", "Mode_Standby"),
               ("Mode_Standby", "Mode_Working"),
               ("Mode_Working", "Mode_Standby"),
               ("Mode_Working", "Mode_Disabled"),
               ("Mode_Standby", "Mode_Disabled")]

for from_state, to_state in transitions:
    gradients = compute_gradients(data, from_state, to_state)
    average_gradient = sum(gradients) / len(gradients) if gradients else 0
    print(f"Average Gradient for transition from {from_state} to {to_state}: {average_gradient:.2f} W/second")
    
    
    
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %% B estimmung Energie pro Teil bzw. Zyklus
# Vorgehen: 10 Teile pro 15 minuten -> Bestimmung Energie pro 15 minuten Zyklus -> Leistung*15 Minuten durch 10 Teile 


#manual Calculation with raw data
import pandas as pd

# Load the dataset (replace 'your_dataset.csv' with the actual file name)
data = pd.read_csv('Data310823_merge.csv')

# Make sure the 'Timestamp' column is in the datetime format
data['Timestamp'] = pd.to_datetime(data['Timestamp'])

# Define the start and end of the cycle (replace with the actual timestamps)
cycle_start = '2023-08-31 05:29:00'
cycle_end = '2023-08-31 05:44:00'

# Convert the cycle start and end to datetime
cycle_start = pd.to_datetime(cycle_start)
cycle_end = pd.to_datetime(cycle_end)

# Filter the dataset to only include the data within the cycle
cycle_data = data[(data['Timestamp'] >= cycle_start) & (data['Timestamp'] <= cycle_end)].copy()

cycle_data.loc[:, 'Time_diff'] = cycle_data['Timestamp'].diff().dt.total_seconds().fillna(0)
cycle_data.loc[:, 'Energy'] = cycle_data['Pges'] * cycle_data['Time_diff']  # Power in kW and Time in seconds, Energy in kWs


# Sum up the total energy in the cycle
total_energy_kWs = cycle_data['Energy'].sum()  # Energy in kilowatt-seconds (kWs)
total_energy_kWh = total_energy_kWs / 3600  # Convert kWs to kilowatt-hours (kWh)
total_energy_part = total_energy_kWh / 10 # Energy per part 

print(f'Total energy consumed in the cycle: {total_energy_kWh} kWh')
print(f'Total energy consumed per Part: {total_energy_part} kWh')
print('Use this Value because auf its most acuarate, while if use the base from Average it is only 50 kWh per Part which is 20 percent wrong')



##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %% Create final EFDM für EMAG GT
print(f"P_working: {clustercenters[3][1]} kW mit T_H_min 90 sec T_H_max 1365") # 90 manuell berechnet, kleinster identifizierter wert 403
print(f"P_standby: {clustercenters[1][1]} kW mit T_H_min 2 T_H_max 1800 sec") # 1800 weil max 30m min in Standby dann automatisch in adisabled
print(f"P_disabled: {clustercenters[2][1]} kW mit T_H_min 0 T_H_max inf")


print(f" Reaction Duration min: {reaction_duration_min}, Reaction Duration max: {reaction_duration_max}")
#print(f" T_R for Regeneration Duration is zero, while switch from Working in a flex mode to reference working mode is direct possible: The regeneration duration tR denotes the time until the FLM can be activated again after deactivation. In the case of the grinding machine, this waiting period is set to zero since all relevant constraints (the remaining production time and latency in the reaction duration; the gradients and workpiece storage capacity; and the production plan) are already covered in other key figures of the grinding machine or workpiece storage data model.")
print(f"T_Reg for Regeneration Duration is {reaction_duration_max} s")
    #   while switch from Working in a flex
    #   mode to reference working mode is direct possible: The regeneration duration tR 
    #   denotes the time until the FLM can be activated again after deactivation.
    #   In the case of the grinding machine, this waiting period is set to maximal Value of T_reaction
    #   since all relevant constraints (the remaining production time and latency in the reaction duration;
    #   the gradients and workpiece storage capacity; and the production plan) 
    #   are already covered in other key figures of the grinding machine or workpiece storage data model.")





##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
##### >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<<<<<<<<<<<<<<<<<<<<<<<<<<<
# %%
#  PLOT FLEX SPACE EMAG GT
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.colors as mcolors


def plot_machine(ax, leistung, dauer, aenderungsrate, unabhaengigkeit, cmap, norm):
    heat = plt.cm.gray(norm(unabhaengigkeit))
    #heat = plt.cm.plasma(unabhaengigkeit)
    
    # Plot the top surface
    Z = aenderungsrate * np.ones_like(leistung)
    ax.plot_surface(leistung, dauer, Z, color=heat,shade=False)#, alpha=1)
    
    # Plot the bottom surface
    Z = np.zeros_like(leistung)
    ax.plot_surface(leistung, dauer, Z, color=heat,shade=False, alpha=0.5)
    
    # Plot vertical edges to create a 3D volume
    for i in [0, 1]:
        for j in [0, 1]:
            ax.plot([leistung[i, j], leistung[i, j]], 
                    [dauer[i, j], dauer[i, j]], 
                    [0, aenderungsrate], 
                    color=heat)
# Daten berechnet aus EFDM 
leistung1 = np.array([[0.2, 2.289], [0.2, 2.289]]) # Power state max von WZM GT
dauer1 = np.array([[0.5, 0.5], [0, 0]]) # Holding Duration max / min in h 
aenderungsrate1 = 0.320 # Gradienten eintragen
unabhaengigkeit1 = 0.4 #lagged Correlation Wert

#plt.rcParams['font.size'] = 12
# %%
# Erstellen des Plots für WZM
fig = plt.figure(figsize=(10, 10))
cmap = plt.cm.plasma
norm = mcolors.Normalize(vmin=0, vmax=1)

# Plot für Maschine 1
ax = fig.add_subplot(111, projection='3d')
plot_machine(ax, leistung1, dauer1, aenderungsrate1, unabhaengigkeit1, cmap, norm)
#ax.set_title("WZM EMAG GT")

# Skalierung anpassen
ax.set_xlim([0, 3]) #Power
ax.set_ylim([0, 1]) # Dauer
ax.set_zlim([0, 0.6]) # Changerate

# Labels setzen
ax.set_xlabel('Power in kW', fontsize=16)
ax.set_ylabel('Duration in h',fontsize=16)
ax.set_zlabel(r'$\nabla$P in kW/s',fontsize=16)
ax.tick_params(axis='both', labelsize=16) 
plt.tight_layout(pad=0.1)

plt.savefig("flexspace_emag.svg", format='svg')
#-->plt.savefig('flexspace_emag.png')

plt.show()



