pukpr/lte_results_manifold.py

## gistfile1.txt
#!/usr/bin/env python3
import sys
import json
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path

def read_lpap(subdir):
    """
    Read lpap list from lt.exe.p in the given subdirectory.
    Returns a list of tuples (period_days, amplitude, phase_rad).
    """
    json_path = Path(subdir) / "lt.exe.p"
    if not json_path.is_file():
        print(f"Warning: {json_path} not found, returning empty list.")
        return []

    with open(json_path, 'r') as f:
        data = json.load(f)

    triplets = data.get("lpap", [])
    components = []
    for t in triplets:
        if len(t) >= 3:
            try:
                period = float(t[0])      # first field: period in days
                amp = float(t[1])         # second field: amplitude
                phase = float(t[2])       # third field: phase in radians
                components.append((period, amp, phase))
            except ValueError:
                # skip if conversion fails
                continue
    return components

def compute_composite(components, t_days):
    """
    Given a list of (period, amplitude, phase) and a time array (days),
    return the sum of all sine waves: Σ A_i * sin(2π * t / P_i + φ_i)
    """
    composite = np.zeros_like(t_days, dtype=float)
    for period, amp, phase in components:
        if period > 0:
            composite += amp * np.sin(2.0 * np.pi * t_days / period + phase)
    return composite

def main():
    if len(sys.argv) != 3:
        print("Usage: python plot_composite.py <subdir1> <subdir2>")
        sys.exit(1)

    subdir1, subdir2 = sys.argv[1], sys.argv[2]

    # Read component data from both directories
    comps1 = read_lpap(subdir1)
    comps2 = read_lpap(subdir2)

    if not comps1 and not comps2:
        print("No data found in either file.")
        return

    # Time axis: from 1962‑01‑01 to 1972‑01‑01 (10 years)
    start_year = 1962.0
    duration_years = 10.0
    days_per_year = 365.2422
    num_points = 20000

    t_days = np.linspace(0.0, duration_years * days_per_year, num_points)
    t_years = start_year + t_days / days_per_year

    # Compute composite signals
    if comps1:
        signal1 = compute_composite(comps1, t_days)
    else:
        signal1 = np.zeros_like(t_days)

    if comps2:
        signal2 = compute_composite(comps2, t_days)
    else:
        signal2 = np.zeros_like(t_days)

    # Compute correlation coefficient (Pearson)
    if np.any(signal1) and np.any(signal2):
        corr = np.corrcoef(signal1, signal2)[0, 1]
    else:
        corr = float('nan')

    # Plot
    fig, ax = plt.subplots(figsize=(12, 6))

    if comps1:
        ax.plot(t_years, signal1, linestyle='-', color='blue', label=f'Composite from {subdir1}', linewidth=1.5)
    if comps2:
        ax.plot(t_years, signal2, linestyle='--', color='red', label=f'Composite from {subdir2}', linewidth=1.5)

    ax.set_xlabel('Year')
    ax.set_ylabel('Amplitude')
    ax.set_title(f'Composite sine waves (10‑year span starting 1962) - Correlation: {corr:.4f}')
    ax.legend()
    ax.grid(True, linestyle='--', alpha=0.6)

    plt.tight_layout()
    plt.show()

if __name__ == "__main__":
    main()


## lte_results_manifold.py
#!/usr/bin/env python3
"""
Plot column 1 vs column 4 from lte_results.csv files in two subdirs.
Each series is normalized to unit variance.
Pearson's r over the overlapping time segment is computed and shown in the title.
"""
import sys
import os
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

def main():
    if len(sys.argv) < 3:
        print("Usage: python lte_results_manifold.py <dir1> <dir2> [output_png]")
        sys.exit(1)

    output_png = sys.argv[3] if len(sys.argv) > 3 else None

    COLUMN=3
    TITLE='Manifold'

#    COLUMN=2
#    TITLE='Index'

    # Collect one CSV per directory (take the first found if multiple)
    csv_files = []
    for d in sys.argv[1:3]:
        found = glob.glob(os.path.join(d, "lte_results.csv"))
        if found:
            csv_files.append(found[0])  # pick first

    if len(csv_files) != 2:
        print("Need exactly one lte_results.csv in each of the two directories.")
        sys.exit(1)

    # Read both files
    df1 = pd.read_csv(csv_files[0], header=None)
    df2 = pd.read_csv(csv_files[1], header=None)
    x1, y1 = df1[0].values, df1[COLUMN].values
    x2, y2 = df2[0].values, df2[COLUMN].values

    # Normalize each y to unit variance (one‑liner)
    y1 = (y1 - y1.mean()) / y1.std()
    y2 = (y2 - y2.mean()) / y2.std()

    # ---- ONE‑LINER for Pearson's r over the shared time segment ----
    # Defines common x as the union of points inside the overlap, interpolates both,
    # computes correlation, and stores it in 'r'
    x_min, x_max = max(x1.min(), x2.min()), min(x1.max(), x2.max())
    common_x = np.unique(np.concatenate([x1[(x1 >= x_min) & (x1 <= x_max)], x2[(x2 >= x_min) & (x2 <= x_max)]]))
    r = np.corrcoef(np.interp(common_x, x1, y1), np.interp(common_x, x2, y2))[0, 1] if len(common_x) > 1 else np.nan
    # -----------------------------------------------------------------


    # Plot both series
    fig, ax = plt.subplots(figsize=(10, 6))

    # Add vertical gray dotted grid lines every 3.81 years
    spacing = 3.81
    grid_years = np.arange(1880, 2021, spacing)
    for g in grid_years:
        ax.axvline(g, color='green', linestyle=':', lw=0.8)

    ax.plot(x1, y1, label=sys.argv[1],  color='red', linestyle='-', linewidth=1)
    ax.plot(x2, y2, label=sys.argv[2], color='blue', linestyle='--', linewidth=1)


    # Place r in the title (with 3 decimals)
    title = f"{TITLE} — Pearson r = {r:.3f}" if not np.isnan(r) else "No overlapping time segment"
    ax.set_title(title)
    ax.set_xlabel("Year")
    ax.set_ylabel(TITLE)
    ax.legend()
    # ax.grid(True, linestyle='--', alpha=0.7)
    plt.tight_layout()

    if output_png:
        plt.savefig(output_png, dpi=100, bbox_inches='tight')
        print(f"Saved: {output_png}")
    else:
        plt.show()

if __name__ == "__main__":
    main()
	#!/usr/bin/env python3
	import sys
	import json
	import numpy as np
	import matplotlib.pyplot as plt
	from pathlib import Path

	def read_lpap(subdir):
	"""
	Read lpap list from lt.exe.p in the given subdirectory.
	Returns a list of tuples (period_days, amplitude, phase_rad).
	"""
	json_path = Path(subdir) / "lt.exe.p"
	if not json_path.is_file():
	print(f"Warning: {json_path} not found, returning empty list.")
	return []

	with open(json_path, 'r') as f:
	data = json.load(f)

	triplets = data.get("lpap", [])
	components = []
	for t in triplets:
	if len(t) >= 3:
	try:
	period = float(t[0]) # first field: period in days
	amp = float(t[1]) # second field: amplitude
	phase = float(t[2]) # third field: phase in radians
	components.append((period, amp, phase))
	except ValueError:
	# skip if conversion fails
	continue
	return components

	def compute_composite(components, t_days):
	"""
	Given a list of (period, amplitude, phase) and a time array (days),
	return the sum of all sine waves: Σ A_i * sin(2π * t / P_i + φ_i)
	"""
	composite = np.zeros_like(t_days, dtype=float)
	for period, amp, phase in components:
	if period > 0:
	composite += amp * np.sin(2.0 * np.pi * t_days / period + phase)
	return composite

	def main():
	if len(sys.argv) != 3:
	print("Usage: python plot_composite.py <subdir1> <subdir2>")
	sys.exit(1)

	subdir1, subdir2 = sys.argv[1], sys.argv[2]

	# Read component data from both directories
	comps1 = read_lpap(subdir1)
	comps2 = read_lpap(subdir2)

	if not comps1 and not comps2:
	print("No data found in either file.")
	return

	# Time axis: from 1962‑01‑01 to 1972‑01‑01 (10 years)
	start_year = 1962.0
	duration_years = 10.0
	days_per_year = 365.2422
	num_points = 20000

	t_days = np.linspace(0.0, duration_years * days_per_year, num_points)
	t_years = start_year + t_days / days_per_year

	# Compute composite signals
	if comps1:
	signal1 = compute_composite(comps1, t_days)
	else:
	signal1 = np.zeros_like(t_days)

	if comps2:
	signal2 = compute_composite(comps2, t_days)
	else:
	signal2 = np.zeros_like(t_days)

	# Compute correlation coefficient (Pearson)
	if np.any(signal1) and np.any(signal2):
	corr = np.corrcoef(signal1, signal2)[0, 1]
	else:
	corr = float('nan')

	# Plot
	fig, ax = plt.subplots(figsize=(12, 6))

	if comps1:
	ax.plot(t_years, signal1, linestyle='-', color='blue', label=f'Composite from {subdir1}', linewidth=1.5)
	if comps2:
	ax.plot(t_years, signal2, linestyle='--', color='red', label=f'Composite from {subdir2}', linewidth=1.5)

	ax.set_xlabel('Year')
	ax.set_ylabel('Amplitude')
	ax.set_title(f'Composite sine waves (10‑year span starting 1962) - Correlation: {corr:.4f}')
	ax.legend()
	ax.grid(True, linestyle='--', alpha=0.6)

	plt.tight_layout()
	plt.show()

	if __name__ == "__main__":
	main()
	#!/usr/bin/env python3
	"""
	Plot column 1 vs column 4 from lte_results.csv files in two subdirs.
	Each series is normalized to unit variance.
	Pearson's r over the overlapping time segment is computed and shown in the title.
	"""
	import sys
	import os
	import glob
	import numpy as np
	import pandas as pd
	import matplotlib.pyplot as plt

	def main():
	if len(sys.argv) < 3:
	print("Usage: python lte_results_manifold.py <dir1> <dir2> [output_png]")
	sys.exit(1)

	output_png = sys.argv[3] if len(sys.argv) > 3 else None

	COLUMN=3
	TITLE='Manifold'

	# COLUMN=2
	# TITLE='Index'

	# Collect one CSV per directory (take the first found if multiple)
	csv_files = []
	for d in sys.argv[1:3]:
	found = glob.glob(os.path.join(d, "lte_results.csv"))
	if found:
	csv_files.append(found[0]) # pick first

	if len(csv_files) != 2:
	print("Need exactly one lte_results.csv in each of the two directories.")
	sys.exit(1)

	# Read both files
	df1 = pd.read_csv(csv_files[0], header=None)
	df2 = pd.read_csv(csv_files[1], header=None)
	x1, y1 = df1[0].values, df1[COLUMN].values
	x2, y2 = df2[0].values, df2[COLUMN].values

	# Normalize each y to unit variance (one‑liner)
	y1 = (y1 - y1.mean()) / y1.std()
	y2 = (y2 - y2.mean()) / y2.std()

	# ---- ONE‑LINER for Pearson's r over the shared time segment ----
	# Defines common x as the union of points inside the overlap, interpolates both,
	# computes correlation, and stores it in 'r'
	x_min, x_max = max(x1.min(), x2.min()), min(x1.max(), x2.max())
	common_x = np.unique(np.concatenate([x1[(x1 >= x_min) & (x1 <= x_max)], x2[(x2 >= x_min) & (x2 <= x_max)]]))
	r = np.corrcoef(np.interp(common_x, x1, y1), np.interp(common_x, x2, y2))[0, 1] if len(common_x) > 1 else np.nan
	# -----------------------------------------------------------------



	# Plot both series
	fig, ax = plt.subplots(figsize=(10, 6))

	# Add vertical gray dotted grid lines every 3.81 years
	spacing = 3.81
	grid_years = np.arange(1880, 2021, spacing)
	for g in grid_years:
	ax.axvline(g, color='green', linestyle=':', lw=0.8)

	ax.plot(x1, y1, label=sys.argv[1], color='red', linestyle='-', linewidth=1)
	ax.plot(x2, y2, label=sys.argv[2], color='blue', linestyle='--', linewidth=1)


	# Place r in the title (with 3 decimals)
	title = f"{TITLE} — Pearson r = {r:.3f}" if not np.isnan(r) else "No overlapping time segment"
	ax.set_title(title)
	ax.set_xlabel("Year")
	ax.set_ylabel(TITLE)
	ax.legend()
	# ax.grid(True, linestyle='--', alpha=0.7)
	plt.tight_layout()

	if output_png:
	plt.savefig(output_png, dpi=100, bbox_inches='tight')
	print(f"Saved: {output_png}")
	else:
	plt.show()

	if __name__ == "__main__":
	main()