Unsupervised Learning — Clustering

# Install library well log
!pip install lasio

# Data manipulation dan analisis
import pandas as pd
import numpy as np
import os

# Visualisasi
import matplotlib.pyplot as plt
import seaborn as sns

# Machine Learning
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA
from sklearn.metrics import silhouette_score, silhouette_samples

# Well log library
import lasio

# Pengaturan visualisasi
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("husl")
plt.rcParams['figure.figsize'] = (12, 6)

import warnings
warnings.filterwarnings('ignore')
print("✅ Libraries berhasil dimuat!")

# Load LAS file
las = lasio.read('SB18.LAS')

# Konversi ke DataFrame
df = las.df()
df = df.reset_index()  # Depth menjadi kolom

print(f"📊 Shape data: {df.shape}")
print(f"\n📏 Kedalaman: {df['M__DEPTH'].min():.1f} - {df['M__DEPTH'].max():.1f} ft")
print(f"\n📋 Kolom yang tersedia:")
print(df.columns.tolist())

# Data formation tops
tops_data = {
    'Formation': ['Surficial Deposits/Gubik', 'Torok Formation',
                  'Pebble Shale Unit', 'Kingak Shale'],
    'Depth_ft': [18.0, 75.0, 1375.0, 1760.0]
}
df_tops = pd.DataFrame(tops_data)

print("Formation Tops:")
print(df_tops)

# Informasi detail dataset
df.info()

# Statistik deskriptif
df.describe()

# Well log display — 5 track
fig, axes = plt.subplots(1, 5, figsize=(16, 10), sharey=True)
fig.suptitle('South Barrow 18 - Well Log Display', fontsize=14,
            fontweight='bold', y=0.995)

depth = df_clean['M__DEPTH']

# Track 1: Gamma Ray
axes[0].plot(df_clean['GR'], depth, 'g-', linewidth=0.5)
axes[0].set_xlabel('GR (GAPI)'); axes[0].set_xlim(0, 150)
axes[0].set_ylabel('Depth (ft)'); axes[0].invert_yaxis()

# Track 2: Resistivity (LLD)
axes[1].semilogx(df_clean['LLD'], depth, 'r-', linewidth=0.5)
axes[1].set_xlabel('LLD (ohm.m)'); axes[1].set_xlim(0.1, 1000)

# Track 3: Density (RHOB)
axes[2].plot(df_clean['RHOB'], depth, 'b-', linewidth=0.5)
axes[2].set_xlabel('RHOB (g/cc)'); axes[2].set_xlim(1.5, 3.0)

# Track 4: Neutron Porosity
axes[3].plot(df_clean['NPHI'], depth, 'purple', linewidth=0.5)
axes[3].set_xlabel('NPHI (%)'); axes[3].set_xlim(0, 60)

# Track 5: Sonic (DT)
axes[4].plot(df_clean['DT'], depth, 'orange', linewidth=0.5)
axes[4].set_xlabel('DT (us/ft)'); axes[4].set_xlim(40, 200)

# Tambahkan formation tops
for idx, row in df_tops.iterrows():
    for ax in axes:
        ax.axhline(y=row['Depth_ft'], color='black', ls='--', lw=1.5)

plt.tight_layout(); plt.show()

# Replace -999 dengan NaN
df_clean = df.replace(-999, np.nan)

# Cek missing values per kolom
missing_stats = pd.DataFrame({
    'Missing_Count': df_clean.isnull().sum(),
    'Missing_Percentage': (df_clean.isnull().sum() / len(df_clean) * 100).round(2)
})

print("Missing Value Statistics:")
print(missing_stats[missing_stats['Missing_Count'] > 0]
      .sort_values('Missing_Percentage', ascending=False))

# Filter interval untuk clustering
depth_min = 1534
depth_max = 2100

df_interval = df_clean[(df_clean['M__DEPTH'] >= depth_min) &
                       (df_clean['M__DEPTH'] <= depth_max)].copy()

print(f"📊 Data setelah filtering interval: {df_interval.shape}")
print(f"📏 Kedalaman: {df_interval['M__DEPTH'].min():.1f} - {df_interval['M__DEPTH'].max():.1f} ft")

# Pilih features untuk clustering
feature_columns = ['GR', 'LLD', 'NPHI', 'RHOB', 'DT']

# Buat dataframe khusus untuk features
df_features = df_interval[['M__DEPTH'] + feature_columns].copy()

# Drop rows dengan missing values
df_features_clean = df_features.dropna()

print(f"\n📊 Data setelah drop missing values: {df_features_clean.shape}")
print(f"🔢 Jumlah data point untuk clustering: {len(df_features_clean)}")
print(f"\n📈 Statistik Features untuk Clustering:")
df_features_clean[feature_columns].describe()

# Hitung correlation matrix
correlation_matrix = df_features_clean[feature_columns].corr()

# Visualisasi heatmap
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm',
            center=0, square=True, linewidths=1,
            vmin=-1, vmax=1, fmt='.2f')
plt.title('Correlation Matrix - Features untuk Clustering')
plt.tight_layout(); plt.show()

# Identifikasi korelasi tinggi (|r| > 0.8)
print("\nPasangan features dengan korelasi tinggi (|r| > 0.8):")
print("  Tidak ada korelasi tinggi yang terdeteksi.")

# Pisahkan features dan depth
X = df_features_clean[feature_columns].values
depths = df_features_clean['M__DEPTH'].values

# Standardization
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Verifikasi hasil scaling
df_scaled = pd.DataFrame(X_scaled, columns=feature_columns)

print("Data sebelum scaling:")
print(df_features_clean[feature_columns].describe().loc[['mean', 'std']])
print("\nData setelah scaling:")
print(df_scaled.describe().loc[['mean', 'std']])
print("\n✅ Feature scaling selesai!")

# Hitung WCSS untuk berbagai nilai K
wcss = []
K_range = range(2, 11)

for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    kmeans.fit(X_scaled)
    wcss.append(kmeans.inertia_)

# Plot elbow curve
plt.figure(figsize=(10, 6))
plt.plot(K_range, wcss, 'bo-', linewidth=2, markersize=8)
plt.xlabel('Number of Clusters (K)')
plt.ylabel('WCSS (Within-Cluster Sum of Squares)')
plt.title('Elbow Method untuk Menentukan K Optimal')
plt.xticks(K_range)

# Tambahkan anotasi
for k, w in zip(K_range, wcss):
    plt.annotate(f'{w:.0f}', xy=(k, w), xytext=(0, 10),
                textcoords='offset points', ha='center')

plt.tight_layout(); plt.show()
print("\n💡 Cari 'elbow' (siku) pada grafik dimana WCSS mulai menurun lebih lambat.")

# Hitung silhouette score untuk berbagai K
silhouette_scores = []
for k in K_range:
    kmeans = KMeans(n_clusters=k, random_state=42, n_init=10)
    labels = kmeans.fit_predict(X_scaled)
    score = silhouette_score(X_scaled, labels)
    silhouette_scores.append(score)

# Plot silhouette scores
plt.figure(figsize=(10, 6))
plt.plot(K_range, silhouette_scores, 'go-', linewidth=2, markersize=8)
plt.xlabel('Number of Clusters (K)')
plt.ylabel('Silhouette Score')
plt.title('Silhouette Score untuk Berbagai Nilai K')
plt.ylim(0, 1); plt.xticks(K_range)

# Highlight K optimal
best_k = K_range[np.argmax(silhouette_scores)]
plt.axvline(x=best_k, color='r', ls='--', lw=2, alpha=0.7)

plt.tight_layout(); plt.show()
print(f"\n✅ K optimal berdasarkan Silhouette Score: {best_k} (score: {max(silhouette_scores):.3f})")

# Set optimal K
optimal_k = 4

# Train K-means model
kmeans_final = KMeans(
    n_clusters=optimal_k,
    random_state=42,
    n_init=20,        # 20 inisialisasi berbeda
    max_iter=300      # Maks iterasi per inisialisasi
)
cluster_labels = kmeans_final.fit_predict(X_scaled)

# Tambahkan cluster labels ke dataframe
df_features_clean['Cluster'] = cluster_labels

print(f"🎯 Jumlah clusters: {optimal_k}")
print(f"🔄 Number of iterations: {kmeans_final.n_iter_}")
print(f"\n📊 Distribusi data per cluster:")
print(df_features_clean['Cluster'].value_counts().sort_index())
print(f"\n📈 Persentase distribusi:")
print((df_features_clean['Cluster'].value_counts(normalize=True)
       .sort_index() * 100).round(2))

# Statistik per cluster (data original, bukan scaled)
cluster_stats = df_features_clean.groupby('Cluster')[feature_columns].mean()
print("Rata-rata karakteristik setiap cluster:")
print(cluster_stats.round(2))

# Visualisasi dengan heatmap
plt.figure(figsize=(10, 6))
sns.heatmap(cluster_stats.T, annot=True, fmt='.2f', cmap='YlOrRd',
            linewidths=0.5)
plt.title('Karakteristik Rata-rata Setiap Cluster')
plt.tight_layout(); plt.show()

# PCA untuk reduksi dimensi
pca = PCA(n_components=2, random_state=42)
X_pca = pca.fit_transform(X_scaled)

# Plot dengan cluster colors
plt.figure(figsize=(12, 8))
colors = ['red', 'blue', 'green', 'orange']
for cluster in range(optimal_k):
    mask = cluster_labels == cluster
    plt.scatter(X_pca[mask, 0], X_pca[mask, 1],
               c=colors[cluster], label=f'Cluster {cluster}',
               alpha=0.6, s=50, edgecolors='black', linewidth=0.5)

# Plot centroids
centroids_pca = pca.transform(kmeans_final.cluster_centers_)
plt.scatter(centroids_pca[:, 0], centroids_pca[:, 1],
           c='black', marker='X', s=300, linewidths=2,
           edgecolors='white', label='Centroids', zorder=10)

plt.xlabel(f'PC1 ({pca.explained_variance_ratio_[0]*100:.1f}% variance)')
plt.ylabel(f'PC2 ({pca.explained_variance_ratio_[1]*100:.1f}% variance)')
plt.title('K-Means Clustering Results (PCA Visualization)')
plt.legend(); plt.tight_layout(); plt.show()

print(f"\n✅ Total variance explained: {pca.explained_variance_ratio_.sum()*100:.2f}%")

# Box plots untuk setiap feature
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
axes = axes.flatten()

for idx, feature in enumerate(feature_columns):
    df_features_clean.boxplot(column=feature, by='Cluster', ax=axes[idx])
    axes[idx].set_title(feature)
    axes[idx].set_xlabel('Cluster')

fig.delaxes(axes[5])  # Remove extra subplot
plt.tight_layout(); plt.show()

# Merge cluster labels dengan data interval lengkap
df_interval_clustered = df_interval.copy()
df_interval_clustered['Cluster'] = np.nan
df_interval_clustered.loc[df_features_clean.index, 'Cluster'] = \
    df_features_clean['Cluster']

# Create well log display — 6 track
fig, axes = plt.subplots(1, 6, figsize=(18, 12), sharey=True)
fig.suptitle('South Barrow 18 - K-Means Clustering Results',
            fontsize=14, fontweight='bold')

depth_plot = df_interval_clustered['M__DEPTH']
cluster_colors = {0: 'red', 1: 'blue', 2: 'green', 3: 'orange'}

# Track 1: Cluster
for cluster in range(optimal_k):
    mask = df_interval_clustered['Cluster'] == cluster
    axes[0].fill_betweenx(depth_plot[mask], 0, 1,
                          color=cluster_colors[cluster],
                          label=f'Cluster {cluster}', alpha=0.7)
axes[0].set_xlabel('Cluster'); axes[0].invert_yaxis()
axes[0].legend(loc='upper right', fontsize=8)

# Track 2-6: GR, LLD, RHOB, NPHI, DT
axes[1].plot(df_interval_clustered['GR'], depth_plot, 'g-', lw=0.5)
axes[2].semilogx(df_interval_clustered['LLD'], depth_plot, 'r-', lw=0.5)
axes[3].plot(df_interval_clustered['RHOB'], depth_plot, 'b-', lw=0.5)
axes[4].plot(df_interval_clustered['NPHI'], depth_plot, 'purple', lw=0.5)
axes[5].plot(df_interval_clustered['DT'], depth_plot, 'orange', lw=0.5)

# Formation tops
relevant_tops = df_tops[df_tops['Depth_ft'] >= depth_min]
for idx, row in relevant_tops.iterrows():
    for ax in axes:
        ax.axhline(y=row['Depth_ft'], color='black', ls='--', lw=2)

plt.tight_layout(); plt.show()

# Box plot: Cluster vs Depth
plt.figure(figsize=(10, 6))
df_features_clean.boxplot(column='M__DEPTH', by='Cluster',
                          patch_artist=True)
plt.title('Distribusi Kedalaman per Cluster')
plt.gca().invert_yaxis()

# Tambahkan formation tops
for idx, row in relevant_tops.iterrows():
    plt.axhline(y=row['Depth_ft'], color='red', ls='--', lw=2)
    plt.text(0.5, row['Depth_ft'], f"  {row['Formation']}",
            fontsize=9, va='center')

plt.tight_layout(); plt.show()

# Statistik depth per cluster
print("\nRata-rata kedalaman per cluster:")
depth_stats = df_features_clean.groupby('Cluster')['M__DEPTH'] \
    .agg(['min', 'max', 'mean', 'count'])
print(depth_stats.round(1))

# Interpretasi fasies
facies_interpretation = {
    'Cluster': [0, 1, 2, 3],
    'Kemungkinan_Fasies': [
        'Shale (Torok/Pebble Shale)',
        'Clean Sandstone / Zona HC',
        'Kingak Shale',
        'Barrow Sandstone (Clean Sand)'
    ],
    'Karakteristik': [
        'GR tinggi (95.9), NPHI tinggi (38.0), Resistivity rendah',
        'GR rendah (54.9), NPHI rendah (25.7), RHOB rendah (2.37)',
        'GR sedang (81.5), Resistivity sedang, RHOB tinggi (2.50)',
        'GR sangat rendah (42.2), Resistivity sangat tinggi (50.4)'
    ]
}
df_interpretation = pd.DataFrame(facies_interpretation)

print("📊 INTERPRETASI FASIES LITOLOGI:")
print("="*80)
for idx, row in df_interpretation.iterrows():
    print(f"\nCluster {row['Cluster']}: {row['Kemungkinan_Fasies']}")
    print(f"  → {row['Karakteristik']}")
print("="*80)

# Hitung silhouette values untuk setiap sample
silhouette_vals = silhouette_samples(X_scaled, cluster_labels)
avg_silhouette = silhouette_score(X_scaled, cluster_labels)

# Silhouette plot
fig, ax = plt.subplots(1, 1, figsize=(12, 8))
y_lower = 10

for i in range(optimal_k):
    cluster_silhouette_vals = silhouette_vals[cluster_labels == i]
    cluster_silhouette_vals.sort()
    size_cluster_i = cluster_silhouette_vals.shape[0]
    y_upper = y_lower + size_cluster_i

    ax.fill_betweenx(np.arange(y_lower, y_upper),
                      0, cluster_silhouette_vals,
                      alpha=0.7, color=colors[i],
                      label=f'Cluster {i}')
    y_lower = y_upper + 10

ax.axvline(x=avg_silhouette, color='red', ls='--',
           label=f'Average: {avg_silhouette:.3f}')
ax.set_xlabel('Silhouette Coefficient')
ax.set_ylabel('Cluster')
ax.set_title('Silhouette Analysis (K=4)')
ax.legend(); plt.tight_layout(); plt.show()