import numpy as np
import nibabel as nib
import os
import sys

# === Filnavn ===
filename = "sub-01_task-checkerboard_run-1_bold.nii"

# === Sjekk om filen finnes ===
if not os.path.exists(filename):
    print(f"❌ File not found: {filename}")
    sys.exit(1)

# === Last inn NIfTI-fil ===
nii = nib.load(filename)
data = nii.get_fdata()  # Shape should be (x, y, z, t)

print("🔍 NIfTI shape:", data.shape)

# === Sjekk at det er 4D ===
if data.ndim != 4:
    print("❌ Error: File must have 4D shape (x, y, z, t).")
    sys.exit(1)

x, y, z, t = data.shape
print(f"✅ 4D shape confirmed — Spatial: ({x}, {y}, {z}), Time steps: {t}")

# === Sjekk antall tidsskritt ===
if t < 10:
    print("⚠️ Warning: Very few time steps")

# === Sjekk verdier ===
if np.isnan(data).any():
    print("❌ Contains NaN values")
    sys.exit(1)

if np.isinf(data).any():
    print("❌ Contains infinite values")
    sys.exit(1)

# === Statistikk ===
min_val = np.min(data)
max_val = np.max(data)
mean_val = np.mean(data)

print(f"🧪 Min: {min_val:.2f}, Max: {max_val:.2f}, Mean: {mean_val:.2f}")
print("✔️ NIfTI file passed all checks.")

import numpy as np
import nibabel as nib
import sys
import os

# === TERMINALARGUMENT ===
if len(sys.argv) < 2:
    print("❌ Usage: python3 read_nii_to_ψxt.py ")
    sys.exit(1)

nii_file = sys.argv[1]
if not os.path.exists(nii_file):
    print(f"❌ File not found: {nii_file}")
    sys.exit(1)

# === LAST IN NIfTI-FIL ===
nii = nib.load(nii_file)
data = nii.get_fdata()  # (x, y, z, t)

if data.ndim != 4:
    print("❌ Error: NIfTI file does not have 4D shape (x,y,z,t).")
    sys.exit(1)

print("🔄 Original shape (x,y,z,t):", data.shape)

# === KONVERTER TIL Ψ(x,t): FLATTEN spatial dim → shape = (t, voxels)
ψ_xt = data.reshape(-1, data.shape[3]).T
print("✅ Converted shape (t, voxels):", ψ_xt.shape)

# === LAGRE SOM .npy ===
run_id = os.path.basename(nii_file).split(".")[0].replace(".nii", "").replace(".gz", "")
output_file = f"Ψ_xt_{run_id}.npy"
np.save(output_file, ψ_xt)

print(f"✅ Saved Ψ(x,t) to: {output_file}")

import numpy as np
import sys
import os

# === Terminal Argument ===
if len(sys.argv) < 2:
    print("❌ Usage: python3 check_ψxt_file.py <ψxt-file.npy>")
    sys.exit(1)

npy_file = sys.argv[1]
if not os.path.exists(npy_file):
    print(f"❌ File not found: {npy_file}")
    sys.exit(1)

# === Load Ψ(x,t) ===
ψ = np.load(npy_file)
print("✅ Loaded:", npy_file)
print("🔍 Shape of Ψ(x,t):", ψ.shape)

# === Check dimensionality ===
if ψ.ndim != 2:
    print("❌ Error: Ψ(x,t) should be 2D (time × space).")
    sys.exit(1)

T, N = ψ.shape
print(f"🧠 Time steps (T): {T}, Voxels: {N}")

# === Basic statistics ===
print(f"🧪 Min: {np.min(ψ):.3f}, Max: {np.max(ψ):.3f}, Mean: {np.mean(ψ):.3f}")

# === Sanity checks ===
if T < 10:
    print("⚠️ Warning: Very few time steps.")

if np.isnan(ψ).any():
    print("❌ Contains NaN values.")
    sys.exit(1)

if np.isinf(ψ).any():
    print("❌ Contains infinite values.")
    sys.exit(1)

print("✔️ Ψ(x,t) file structure OK.")

import numpy as np
import json

# === You must have these two functions defined ===
# You can implement them yourself or import from HHM library
def collapse_pattern(ψ_xt):
    """
    CollapsePattern operator: projects each time slice Ψₜ into a feature vector.
    Placeholder implementation: use mean across time.
    """
    return np.array([np.mean(ψ_t) for ψ_t in ψ_xt])  # shape: (T,)

def modal_id_score(collapse_vector):
    """
    Modal ID operator: compares modal identity over time.
    Calculates identity and variation scores.
    """
    modal_id = np.corrcoef(collapse_vector)[0, 0] if collapse_vector.ndim == 2 else 1.0
    variation = 1.0 - np.std(collapse_vector) / (np.abs(np.mean(collapse_vector)) + 1e-8)
    mean_norm = float(np.mean(np.abs(collapse_vector)))
    return {
        "modal_id": float(modal_id),
        "variation": float(variation),
        "mean_norm": float(mean_norm)
    }

# === Load Ψ(x,t) ===
filename = "Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy"
ψ_xt = np.load(filename)

# === Apply operators ===
collapse = collapse_pattern(ψ_xt)
score = modal_id_score(collapse)

# === Evaluate result ===
result = {
    "modal_ID_score": score["modal_id"],
    "variation_score": score["variation"],
    "mean_norm_change": score["mean_norm"],
    "passed": score["modal_id"] > 0.95 and score["variation"] > 0.95
}

# === Save result to JSON ===
with open("result.json", "w") as f:
    json.dump(result, f, indent=2)

# === Print outcome ===
print("✅ MT0001 Validation Result:")
for k, v in result.items():
    print(f"   {k}: {v}")

import numpy as np
import json

# === Load Ψ(x,t) ===
filename = "Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy"
ψ_xt = np.load(filename)  # shape: (T, N)

# === Calculate pairwise L2 norm differences between Ψₜ and Ψₜ₊₁
differences = np.linalg.norm(ψ_xt[1:] - ψ_xt[:-1], axis=1)
mean_diff = float(np.mean(differences))
min_diff = float(np.min(differences))
max_diff = float(np.max(differences))

# === Define threshold
threshold = 0.1

# === Result
result = {
    "mean_frame_diff": mean_diff,
    "min_frame_diff": min_diff,
    "max_frame_diff": max_diff,
    "threshold": threshold,
    "passed": mean_diff > threshold
}

# === Save result
with open("result_mt0002.json", "w") as f:
    json.dump(result, f, indent=2)

# === Print
print("✅ MT0002 Validation Result:")
for k, v in result.items():
    print(f"  {k}: {v}")

import numpy as np
import matplotlib.pyplot as plt
import csv
from sklearn.decomposition import PCA
from sklearn.metrics import mean_squared_error

ψ_xt = np.load("Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy")
component_range = range(10, 510, 10)
errors = []

for k in component_range:
    pca = PCA(n_components=k)
    ψ_pca = pca.fit_transform(ψ_xt)
    ψ_recon = pca.inverse_transform(ψ_pca)
    mse = mean_squared_error(ψ_xt, ψ_recon)
    errors.append((k, mse))
    print(f"Components: {k}, MSE: {mse:.2f}")

# Save CSV
with open("mt0003_reconstruction_log.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["PCA_Components", "Reconstruction_MSE"])
    writer.writerows(errors)

# Plot
components = [x[0] for x in errors]
mse_values = [x[1] for x in errors]

plt.figure(figsize=(10, 6))
plt.plot(components, mse_values, marker='o')
plt.axhline(y=100, color='red', linestyle='--', label='MT0003 Threshold (MSE = 100)')
plt.axvline(x=470, color='green', linestyle='--', label='Validation Point (470 components)')
plt.title("Modal Trace Reconstruction Error vs. PCA Components")
plt.xlabel("Number of PCA Components")
plt.ylabel("Reconstruction Error (MSE)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("mt0003_complexity_curve.png")
plt.show()

import numpy as np
import json
from sklearn.metrics.pairwise import cosine_similarity

ψ_xt = np.load("Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy")
ψ_norm = ψ_xt / (np.linalg.norm(ψ_xt, axis=1, keepdims=True) + 1e-8)
similarity_matrix = cosine_similarity(ψ_norm)

T = ψ_xt.shape[0]
lag = 5
mask = np.triu(np.ones_like(similarity_matrix), k=lag)

similarity_threshold = 0.95
recurrence_hits = np.sum((similarity_matrix > similarity_threshold) * mask)
possible_pairs = np.sum(mask)
recurrence_score = recurrence_hits / possible_pairs if possible_pairs > 0 else 0.0

threshold = 0.4
passed = bool(recurrence_score > threshold)

result = {
    "recurrence_score": float(recurrence_score),
    "similarity_threshold": float(similarity_threshold),
    "recurrence_lag_min": int(lag),
    "comparison_pairs": int(possible_pairs),
    "recurring_pairs": int(recurrence_hits),
    "threshold": float(threshold),
    "passed": passed
}

with open("result_mt0004.json", "w") as f:
    json.dump(result, f, indent=2)

print("✅ MT0004 Validation Result:")
for k, v in result.items():
    print(f"  {k}: {v}")

import numpy as np
import csv
import matplotlib.pyplot as plt
from sklearn.metrics.pairwise import cosine_similarity

# === Load and normalize Ψ(x,t)
ψ_xt = np.load("Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy")
ψ_norm = ψ_xt / (np.linalg.norm(ψ_xt, axis=1, keepdims=True) + 1e-8)

# === Cosine similarity matrix
similarity_matrix = cosine_similarity(ψ_norm)
T = ψ_xt.shape[0]
lag = 5
mask = np.triu(np.ones_like(similarity_matrix), k=lag)

# === Define similarity thresholds
thresholds = [round(t, 2) for t in np.arange(0.90, 0.991, 0.01)]
results = []

for sim_thresh in thresholds:
    recurrence_hits = np.sum((similarity_matrix > sim_thresh) * mask)
    possible_pairs = np.sum(mask)
    recurrence_score = recurrence_hits / possible_pairs if possible_pairs > 0 else 0.0
    results.append((sim_thresh, recurrence_score))
    print(f"Similarity: {sim_thresh}, Recurrence Score: {recurrence_score:.4f}")

# === Save to CSV
with open("mt0004_recurrence_profile.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["Similarity_Threshold", "Recurrence_Score"])
    writer.writerows(results)

# === Plot recurrence curve
thresholds_only = [r[0] for r in results]
scores_only = [r[1] for r in results]
pass_threshold = 0.4
max_sim_valid = max([t for t, s in results if s >= pass_threshold])

plt.figure(figsize=(10, 6))
plt.plot(thresholds_only, scores_only, marker='o')
plt.axhline(y=pass_threshold, color='red', linestyle='--', label=f'MT0004 Validation Line (Score ≥ {pass_threshold})')
plt.axvline(x=max_sim_valid, color='green', linestyle='--', label=f'Validation Max Similarity = {max_sim_valid:.2f}')
plt.xlabel("Similarity Threshold")
plt.ylabel("Recurrence Score")
plt.title("MT0004: Modal Recurrence vs Similarity Threshold")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("mt0004_recurrence_curve.png")
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
import csv

filename = "Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy"
ψ_xt = np.load(filename)
T, N = ψ_xt.shape
component_range = list(range(1, 11)) + list(range(20, 210, 10))
threshold = 1.2
results = []

def evaluate(ψ, n_components):
    pca = PCA(n_components=n_components)
    ψ_reduced = pca.fit_transform(ψ)
    X = np.hstack([ψ_reduced[:-2], ψ_reduced[1:-1]])
    Y = ψ_reduced[2:]
    model = Ridge(alpha=1.0)
    model.fit(X, Y)
    Y_pred = model.predict(X)
    mse_true = mean_squared_error(Y, Y_pred)
    mse_shuffled = mean_squared_error(np.random.permutation(Y), Y_pred)
    score = mse_shuffled / mse_true if mse_true > 0 else 0.0
    return score

for n_components in component_range:
    score_real = evaluate(ψ_xt, n_components)
    noise = np.random.normal(0, 1, size=(T, N))
    score_noise = evaluate(noise, n_components)
    results.append((n_components, score_real, score_noise))
    print(f"PCA {n_components} → Ψ(x,t): {score_real:.2f}, Noise: {score_noise:.2f}")

with open("mt0005_causality_noise_comparison.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["PCA_Components", "Causality_Ψ", "Causality_Noise"])
    writer.writerows(results)

components = [r[0] for r in results]
scores_real = [r[1] for r in results]
scores_noise = [r[2] for r in results]

plt.figure(figsize=(10, 6))
plt.plot(components, scores_real, marker='o', color='blue', label="Ψ(x,t)")
plt.plot(components, scores_noise, marker='x', color='black', label="Noise")
plt.axhline(y=threshold, color='red', linestyle='--', label=f"Validation threshold = {threshold}")
plt.title("MT0005 – Modal Causality: Ψ(x,t) vs Noise")
plt.xlabel("PCA Components")
plt.ylabel("Causality Score (Shuffled / True)")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("mt0005_causality_noise_comparison.png")
plt.show()

import numpy as np
import json
from sklearn.decomposition import PCA
from math import factorial
from collections import Counter

# Permutation entropy function
def permutation_entropy(time_series, order=3, delay=1, normalize=True):
    n = len(time_series)
    if n < order * delay:
        return np.nan
    permutations = []
    for i in range(n - delay * (order - 1)):
        window = time_series[i:(i + delay * order):delay]
        ranks = tuple(np.argsort(window))
        permutations.append(ranks)
    counts = Counter(permutations)
    probs = np.array(list(counts.values()), dtype=np.float64)
    probs /= np.sum(probs)
    pe = -np.sum(probs * np.log(probs + 1e-10))
    if normalize:
        pe /= np.log(factorial(order))
    return pe

ψ_xt = np.load("Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy")[:500]
pca_components = 10
pca = PCA(n_components=pca_components)
ψ_pca = pca.fit_transform(ψ_xt)

pe_true = [permutation_entropy(ψ_pca[:, i], order=3) for i in range(pca_components)]
entropy_true_mean = float(np.nanmean(pe_true))

ψ_shuffled = np.copy(ψ_pca)
np.random.shuffle(ψ_shuffled)
pe_shuffled = [permutation_entropy(ψ_shuffled[:, i], order=3) for i in range(pca_components)]
entropy_shuffled_mean = float(np.nanmean(pe_shuffled))

entropy_score = entropy_shuffled_mean / entropy_true_mean if entropy_true_mean > 0 else 0.0
threshold = 1.05
passed = entropy_score > threshold

result = {
  "pca_components": pca_components,
  "sample_length": 500,
  "permutation_entropy_true": entropy_true_mean,
  "permutation_entropy_shuffled": entropy_shuffled_mean,
  "entropy_score": float(entropy_score),
  "threshold": threshold,
  "passed": bool(passed)
}

with open("result_t014.json", "w") as f:
  json.dump(result, f, indent=2)

print("✅ T014 – Permutation Entropy Result:")
for k, v in result.items():
  print(f"  {k}: {v}")

import numpy as np
import json

# === Load both Ψ(x,t) files ===
ψ_baseline = np.load("Ψ_xt_sub-PS50_ses-baseline_pet.npy")
ψ_blocked = np.load("Ψ_xt_sub-PS50_ses-blocked_pet.npy")

# === Check shapes ===
assert ψ_baseline.shape == ψ_blocked.shape
T, N = ψ_baseline.shape

# === Compute difference ΨΔ(x,t) ===
delta = ψ_blocked - ψ_baseline

# === Metrics ===
delta_norms = np.linalg.norm(delta, axis=1)  # per timestep
ΔΨ_mean = float(np.mean(delta_norms))
ΔΨ_total = float(np.linalg.norm(delta))
ΔΨ_score = ΔΨ_mean / (np.mean(np.linalg.norm(ψ_baseline, axis=1)) + 1e-8)
ΔΨ_min = float(np.min(delta_norms))
ΔΨ_max = float(np.max(delta_norms))

# === Result ===
result = {
    "ΔΨ_mean_per_timestep": ΔΨ_mean,
    "ΔΨ_total_norm": ΔΨ_total,
    "ΔΨ_score": ΔΨ_score,
    "ΔΨ_min_timestep": ΔΨ_min,
    "ΔΨ_max_timestep": ΔΨ_max,
    "passed": ΔΨ_score > 0.3
}

with open("result_mt0009.json", "w") as f:
    json.dump(result, f, indent=2)

print("✅ MT0009 Validation Result:")
for k, v in result.items():
    print(f"  {k}: {v}")

import numpy as np
import matplotlib.pyplot as plt
import csv
from sklearn.decomposition import PCA
from sklearn.metrics.pairwise import cosine_distances

# === Load Ψ(x,t) ===
ψ_xt = np.load("Ψ_xt_sub-01_task-checkerboard_run-1_bold.npy")
T, N = ψ_xt.shape
max_lag = 20
threshold = 0.05

# === Option 1: full Ψ(x,t)
dist_full = []
for lag in range(1, max_lag + 1):
    ψ_now = ψ_xt[:T - lag]
    ψ_lag = ψ_xt[lag:]
    cos_dist = cosine_distances(ψ_now, ψ_lag).diagonal()
    dist_full.append((lag, np.mean(cos_dist)))

# === Option 2: PCA-reduced
pca = PCA(n_components=50)
ψ_pca = pca.fit_transform(ψ_xt)
dist_pca = []
for lag in range(1, max_lag + 1):
    ψ_now = ψ_pca[:T - lag]
    ψ_lag = ψ_pca[lag:]
    cos_dist = cosine_distances(ψ_now, ψ_lag).diagonal()
    dist_pca.append((lag, np.mean(cos_dist)))

# === Save CSV
with open("mt0010_full_vs_pca.csv", "w", newline="") as f:
    writer = csv.writer(f)
    writer.writerow(["Lag", "Full_Ψ", "PCA_Ψ_50"])
    for i in range(max_lag):
        writer.writerow([dist_full[i][0], dist_full[i][1], dist_pca[i][1]])

# === Plot
lags = [d[0] for d in dist_full]
full_vals = [d[1] for d in dist_full]
pca_vals = [d[1] for d in dist_pca]

plt.figure(figsize=(10, 6))
plt.plot(lags, full_vals, marker='o', label="Full Ψ(x,t)")
plt.plot(lags, pca_vals, marker='x', label="PCA-reduced (50 comps)")
plt.axhline(y=threshold, color='red', linestyle='--', label=f"Threshold = {threshold}")
plt.title("MT0010 – Cosine Distance Across Lags (Full vs PCA)")
plt.xlabel("Time Lag Δt")
plt.ylabel("Avg Cosine Distance")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.savefig("mt0010_full_vs_pca.png")
plt.show()

import numpy as np
import json
from scipy.spatial.distance import cosine

# === Filnavn for to Ψ(x,t)-filer ===
file_a = "Ψ_xt_subj_A.npy"
file_b = "Ψ_xt_subj_B.npy"

# === Last inn begge filer ===
ψ_a = np.load(file_a)  # shape: (T, N)
ψ_b = np.load(file_b)  # shape: (T, N)

# === Sjekk at de har samme form ===
if ψ_a.shape != ψ_b.shape:
    raise ValueError(f"Shape mismatch: {ψ_a.shape} vs {ψ_b.shape}")

T = ψ_a.shape[0]
cosine_scores = []

# === Beregn gjennomsnittlig cosine-similaritet for hver Ψₜ ===
for t in range(T):
    sim = 1.0 - cosine(ψ_a[t], ψ_b[t])  # cosine similarity = 1 - distance
    cosine_scores.append(sim)

avg_similarity = float(np.mean(cosine_scores))
min_similarity = float(np.min(cosine_scores))
max_similarity = float(np.max(cosine_scores))

# === Terskel for modal resonans ===
threshold = 0.9
passed = avg_similarity > threshold

# === Resultat ===
result = {
    "T": T,
    "avg_similarity": avg_similarity,
    "min_similarity": min_similarity,
    "max_similarity": max_similarity,
    "threshold": threshold,
    "passed": passed
}

# === Lagre resultat ===
with open("result_mt0011.json", "w") as f:
    json.dump(result, f, indent=2)

print("✅ MT0011 Validation Result:")
for k, v in result.items():
    print(f"  {k}: {v}")