# Import necessary libraries
import numpy as np

# --- SECTION: Simulate High-Dimensional Protein Embeddings ---
# In a real-world scenario, you would load pre-computed embeddings from models
# like ESM, ProtBERT, or AlphaFold representations. For this demonstration,
# we generate synthetic embeddings.

# Define the number of protein embeddings (e.g., from a dataset of proteins)
num_proteins = 1000
# Define the original dimensionality of each embedding (e.g., typical for protein language models)
original_dimension = 1024

# Generate synthetic high-dimensional protein embeddings.
# We use random data for demonstration, assuming a normal distribution
# around a mean, simulating features.
# np.random.randn creates samples from a standard normal distribution.
# We then scale and shift it to make it slightly more 'structured' than pure noise.
protein_embeddings = np.random.randn(num_proteins, original_dimension) * 5 + 10

print(f"Simulated protein embeddings shape: {protein_embeddings.shape}")
print(f"First 5 dimensions of the first embedding:\n{protein_embeddings[0, :5]}")

# --- SECTION: Data Preprocessing (Standardization) ---
# Standardizing the data is a crucial prerequisite for PCA.
# PCA is sensitive to the scale of the features, and unscaled data can lead
# to components dominated by features with larger variances.

from sklearn.preprocessing import StandardScaler

# Initialize the StandardScaler. This will compute the mean and standard
# deviation for each feature (dimension) across all embeddings.
scaler = StandardScaler()

# Fit the scaler to our protein embeddings and then transform them.
# The 'fit' step calculates means and stds.
# The 'transform' step applies this scaling.
scaled_embeddings = scaler.fit_transform(protein_embeddings)

print(f"\nShape of scaled embeddings: {scaled_embeddings.shape}")
print(f"Mean of first dimension after scaling: {np.mean(scaled_embeddings[:, 0]):.4f}")
print(f"Standard deviation of first dimension after scaling: {np.std(scaled_embeddings[:, 0]):.4f}")

# Save the scaled embeddings for the next step
np.save('scaled_protein_embeddings.npy', scaled_embeddings)
print("Scaled embeddings saved as 'scaled_protein_embeddings.npy'.")

# Import necessary libraries
import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# --- SECTION: Load Scaled Embeddings ---
# Load the scaled embeddings generated in the previous step.
try:
    scaled_embeddings = np.load('scaled_protein_embeddings.npy')
    print(f"Loaded scaled embeddings with shape: {scaled_embeddings.shape}")
except FileNotFoundError:
    print("Error: 'scaled_protein_embeddings.npy' not found. Please run the previous code block first.")
    exit()

# --- SECTION: Determine Optimal Number of Components with Explained Variance Ratio ---
# We apply PCA initially with all possible components to understand the variance explained.
# This allows us to make an informed decision on the target dimensionality.

# Initialize PCA with None to retain all components (up to min(n_samples, n_features))
full_pca = PCA(n_components=None)

# Fit PCA to the scaled protein embeddings
full_pca.fit(scaled_embeddings)

# Calculate the cumulative explained variance ratio
cum_explained_variance = np.cumsum(full_pca.explained_variance_ratio_)

# Visualize the explained variance to identify the 'elbow point'
plt.figure(figsize=(10, 6))
plt.plot(range(1, len(cum_explained_variance) + 1), cum_explained_variance, marker='o', linestyle='--')
plt.xlabel('Number of Components')
plt.ylabel('Cumulative Explained Variance Ratio')
plt.title('PCA: Explained Variance vs. Number of Components')
plt.grid(True)

# Identify a target number of components, e.g., to retain 95% of variance
target_variance_ratio = 0.95
num_components_95 = np.where(cum_explained_variance >= target_variance_ratio)[0][0] + 1
plt.axvline(x=num_components_95, color='r', linestyle=':', label=f'{target_variance_ratio*100}% Variance at {num_components_95} components')
plt.legend()
plt.show()

print(f"\nTo retain {target_variance_ratio*100}% of variance, we need {num_components_95} components.")

# --- SECTION: Apply PCA with the Chosen Number of Components ---
# Now, we re-initialize PCA with our chosen number of components for the final transformation.
# For demonstration, let's target reducing to 128 dimensions, which is a common practice
# for efficient vector search while retaining significant information. 
# In a real scenario, you'd use `num_components_95` or another determined value.

target_n_components = 128 # Example target dimensionality for search efficiency

# Ensure target_n_components does not exceed available dimensions
if target_n_components > scaled_embeddings.shape[1]:
    target_n_components = scaled_embeddings.shape[1]
    print(f"Warning: Target components adjusted to maximum possible: {target_n_components}")

# Initialize PCA with the chosen number of components
pca = PCA(n_components=target_n_components)

# Fit PCA to the scaled embeddings and transform them
# The 'fit' step learns the principal components (eigenvectors).
# The 'transform' step projects the original data onto these new components.
reduced_embeddings = pca.fit_transform(scaled_embeddings)

print(f"\nOriginal embedding shape: {scaled_embeddings.shape}")
print(f"Reduced embedding shape: {reduced_embeddings.shape}")
print(f"Variance explained by {target_n_components} components: {pca.explained_variance_ratio_.sum():.4f}")

# Save the reduced embeddings for downstream tasks
np.save('reduced_protein_embeddings.npy', reduced_embeddings)
print("Reduced embeddings saved as 'reduced_protein_embeddings.npy'.")

# Import necessary libraries
import numpy as np
from sklearn.metrics.pairwise import cosine_similarity
import time

# --- SECTION: Load Reduced and Original Embeddings ---
# Load the reduced embeddings for efficient search and the original scaled ones for comparison.
try:
    reduced_embeddings = np.load('reduced_protein_embeddings.npy')
    scaled_embeddings = np.load('scaled_protein_embeddings.npy') # Need original for comparison
    print(f"Loaded reduced embeddings with shape: {reduced_embeddings.shape}")
    print(f"Loaded original scaled embeddings with shape: {scaled_embeddings.shape}")
except FileNotFoundError:
    print("Error: Embedding files not found. Please run previous code blocks first.")
    exit()

# --- SECTION: Simulate Vector Search with Cosine Similarity ---
# We'll compare the time taken for a simple similarity search using
# both original (scaled) and reduced embeddings.

# Define a query embedding (e.g., the first embedding from our dataset)
query_index = 0
query_original = scaled_embeddings[query_index].reshape(1, -1)
query_reduced = reduced_embeddings[query_index].reshape(1, -1)

print(f"\nPerforming similarity search for query embedding at index {query_index}...")

# Measure time for search using original (scaled) embeddings
start_time_original = time.time()
similarities_original = cosine_similarity(query_original, scaled_embeddings)
end_time_original = time.time()
time_original = end_time_original - start_time_original

# Measure time for search using reduced embeddings
start_time_reduced = time.time()
similarities_reduced = cosine_similarity(query_reduced, reduced_embeddings)
end_time_reduced = time.time()
time_reduced = end_time_reduced - start_time_reduced

print(f"Time taken for search with ORIGINAL embeddings: {time_original:.6f} seconds")
print(f"Time taken for search with REDUCED embeddings: {time_reduced:.6f} seconds")

# Verify that the most similar protein is still the query itself (index 0)
# For original
top_k_original_indices = np.argsort(similarities_original[0])[::-1][:5] # Top 5
print(f"Top 5 similar (Original): {top_k_original_indices}")
print(f"Similarity scores (Original): {[similarities_original[0][i] for i in top_k_original_indices]}")

# For reduced
top_k_reduced_indices = np.argsort(similarities_reduced[0])[::-1][:5] # Top 5
print(f"Top 5 similar (Reduced): {top_k_reduced_indices}")
print(f"Similarity scores (Reduced): {[similarities_reduced[0][i] for i in top_k_reduced_indices]}")

# --- SECTION: Discussing Downstream Impact and Storage ---
print("\n--- Storage and Efficiency Impact ---")
original_size_mb = scaled_embeddings.nbytes / (1024 * 1024)
reduced_size_mb = reduced_embeddings.nbytes / (1024 * 1024)
print(f"Original embeddings storage: {original_size_mb:.2f} MB")
print(f"Reduced embeddings storage: {reduced_size_mb:.2f} MB")
print(f"Storage reduction factor: {original_size_mb / reduced_size_mb:.2f}x")

# We can also conceptually show how reduced dimensions benefit Approximate Nearest Neighbors (ANN)
print("\nReduced dimensionality significantly benefits Approximate Nearest Neighbor (ANN) algorithms.")
print("Algorithms like Faiss, Annoy, or HNSW perform better with lower-dimensional data ")
print("due to reduced memory footprint, faster distance calculations, and improved index build times.")
print("This directly translates to quicker identification of similar proteins in large datasets.")

# Import necessary libraries
import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.manifold import TSNE

# --- SECTION: Load Reduced Embeddings ---
# Load the reduced embeddings for validation.
try:
    reduced_embeddings = np.load('reduced_protein_embeddings.npy')
    print(f"Loaded reduced embeddings with shape: {reduced_embeddings.shape}")
except FileNotFoundError:
    print("Error: 'reduced_protein_embeddings.npy' not found. Please run previous code blocks first.")
    exit()

# --- SECTION: Evaluate Biological Signal Preservation (Clustering) ---
# We perform a simple K-Means clustering to see if biologically relevant clusters
# can still be identified in the reduced space.
# In a real scenario, you'd have ground truth labels to evaluate clustering performance.

# For demonstration, let's assume we expect 5 distinct protein families/clusters
num_clusters = 5

# Initialize and fit K-Means clustering on the reduced embeddings
kmeans = KMeans(n_clusters=num_clusters, random_state=42, n_init=10) # n_init for robustness
cluster_labels = kmeans.fit_predict(reduced_embeddings)

print(f"\nK-Means clustering performed on reduced embeddings. Found {num_clusters} clusters.")
print(f"First 10 cluster labels: {cluster_labels[:10]}")

# --- SECTION: Visualize Reduced Embeddings (t-SNE for 2D/3D visualization) ---
# PCA itself can reduce to 2 or 3 components for direct visualization,
# but t-SNE (or UMAP) is often better for visualizing complex high-dimensional structures in 2D.

# Reduce to 2 components for t-SNE visualization
# Note: t-SNE is computationally intensive for very large datasets.
# For demonstration, we might sample a subset if original `num_proteins` was very large.

# Let's perform t-SNE on the reduced_embeddings for better cluster separation visualization
print("\nApplying t-SNE for 2D visualization...")
tsne = TSNE(n_components=2, random_state=42, perplexity=30, n_iter=1000, learning_rate='auto')
tsne_results = tsne.fit_transform(reduced_embeddings)

plt.figure(figsize=(10, 8))
sns.scatterplot(
    x=tsne_results[:, 0],
    y=tsne_results[:, 1],
    hue=cluster_labels, # Color points by their assigned cluster
    palette=sns.color_palette('tab10', num_clusters),
    legend='full',
    alpha=0.8
)
plt.title('t-SNE Visualization of PCA-Reduced Protein Embeddings (Colored by K-Means Clusters)')
plt.xlabel('t-SNE Component 1')
plt.ylabel('t-SNE Component 2')
plt.grid(True)
plt.show()

print("\nThis visualization helps to qualitatively assess if clusters are well-separated in the reduced space.")

# --- SECTION: Discussing Future Frontiers and Advanced Techniques ---
print("\n--- Strategic Implications and Future Directions ---")
print("•  **High-throughput Screening**: Faster similarity searches enable rapid filtering of vast protein libraries.")
print("•  **Drug Discovery**: Quickly identify potential drug targets or lead compounds based on protein interaction profiles.")
print("•  **Protein Engineering**: Design novel proteins with desired functions by exploring the reduced embedding space.")
print("•  **Combining Methods**: Explore integrating PCA with non-linear methods like UMAP or Autoencoders for more nuanced dimensionality reduction.")
print("•  **Incremental PCA**: For extremely large datasets that don't fit in memory, Incremental PCA can process data in mini-batches.")