#! /usr/bin/python3 
# Last edited on 2026-10-23 20:37:40 by stolfi

import sys, re, os
from sys import stderr as err, stdin as inn, stdout as out
from math import sqrt, sin, cos, log, exp, floor, pi, inf, isfinite
import numpy

def main(np,nc,fmt):
  # Reads from {stdin} a list of {np} points, one per line.
  # Each linemust have the same number {nc} of float coordinates.
  # 
  # Outputs to stdout {nc+1} lines each with {nc} coordinates.
  # The first line is the average {pavg} of the data points.
  # The next {nc} lines are the eigenvectors of the covariance
  # matrix (the semiaxes of the Gaussian ellipsoid that fits 
  # the displacements from that average).  All coordinate are 
  # written with the format {fmt}, which should be good for a
  # single float value.
  
  D = numpy.zeros((np,nc,))
  kp = 0
  for (iline,line) in enumerate(inn.readlines()):
    nline = iline + 1
    if re.match(r" *([#]|$)", line): continue
    assert kp < np,  f"stdin:{nline} ** too many data points"
    p = line.split()
    if len(p) != nc: assert False, f"stdin:{nline}: ** wrong line size [[{line}]]"
    p = [ float(pj) for pj in p ]
    D[kp,:] = p
    kp += 1
  assert kp == np, f"stdin:{nline} ** too few data points {kp} {np}"
  
  pavg = get_point_cloud_mean(D)
  Evec, eval = get_point_cloud_axes(D, pavg)
  assert numpy.shape(Evec) == (nc,nc,)
  assert numpy.shape(eval) == (nc,)
  write_vec(pavg,nc,fmt)
  for ke in range(nc):
    evalk = eval[ke]
    assert evalk >= 0
    sigmak = sqrt(evalk)
    eveck = Evec[ke,:]
    esum = 0;
    for jc in range(nc): esum += eveck[jc]
    if esum < 0: sigmak = -sigmak
    eveck = sigmak*eveck
    write_vec(eveck,nc,fmt)
  return
  # ----------------------------------------------------------------------
  
def write_vec(v,nc,fmt):
  for jc in range(nc):
    if jc > 0: out.write(" ")
    out.write(fmt % v[jc])
  out.write("\n")
  # ----------------------------------------------------------------------

def get_point_cloud_mean(D):
  # Expects {D} to be a matrix of point data. It must have shape
  # {np,nc}, where {np} is the number of data points and {nc} is the
  # number of data coordinates.
  #
  # Returns the average {pavg} of all data points.
  
  np, nc = numpy.shape(D)
  
  err.write(f"computing mean of data points ...\n")
  pavg = numpy.sum(D, 0)/numpy.size(D, 0)
  err.write("max(abs(pavg)) = %.6f\n" % numpy.max(numpy.abs(pavg)))
  assert numpy.shape(pavg) == (nc,), "bug: pavg shape"
  return pavg
  # ----------------------------------------------------------------------

def get_point_cloud_axes(D, pctr):
  # Expects {D} to be a matrix of point data. It must have shape
  # {np,nc}, where {np} is the number of data points and {nc} is the
  # number of data coordinates.
  #
  # Computes the covariance matrix {V} of those data points, with shape {(nc,nc)},
  # relative to the given center {pctr}. Then computes the eigenvectors of
  # {V} as the COLUMNS of a matrix {evec} with shape {(nc,nc)}, and the
  # vector of the respective eigenvalues {eval}. 
  #
  # Returns {evec, eval}.
  
  np, nc = numpy.shape(D)
  
  # Compute {D} whose rows are the displacements from barycenter:
  err.write(f"subtracting given center from data points ...\n")
  D = numpy.array([ D[kp,:] - pctr for kp in range(np) ])
  err.write("max(abs(D)) = %.6f\n" % numpy.max(numpy.abs(D)))
  
  # Compute the covariance matrix {V=D'*D}:
  V = (D.transpose() @ D) / float(np)
  err.write("max(V) = %.6f\n" % numpy.max(V))

  # Compute the eigenvalues {evec} and eigenvectors {eval}:
  err.write(f"computing eigenvalues of covariance matrix ...\n")
  eig = numpy.linalg.eigh(V)

  eval = eig.eigenvalues
  err.write(f"eigenvalues:\n{eval}\n");
  assert numpy.shape(eval) == (nc,), "bug: eigenvalues shape"

  evec = eig.eigenvectors
  err.write(f"eigenvectors:\n{evec}\n")
  assert numpy.shape(evec) == (nc,nc,), "bug: eigenvectors shape"
  return evec, eval
  # ----------------------------------------------------------------------

iarg = 1
np = int(sys.argv[iarg]); iarg += 1
nc = int(sys.argv[iarg]); iarg += 1
fmt = sys.argv[iarg]; iarg += 1
main(np,nc,fmt)