#ifndef rmxn_ellipsoid_H
#define rmxn_ellipsoid_H

/* General tools and concepts for geometry of ellipsoids. */
/* Last edited on 2024-12-05 10:28:31 by stolfi */ 

#include <stdio.h>
#include <stdint.h>

#include <bool.h>
#include <r2.h>
#include <i2.h>

/* 
  AXIS-ALIGNED ELLIPSOID
  
  Let {arad[0..n-1]} be an /axial radius vector/, a list of non-negative
  numbers. In what follows we denote by {\RE(arad)} the {n}-dimensional
  ellipsoid in {\RR^n} whose center is the origin, whose main axes are
  aligned with the coordinate axes of {\RR^n}, and whose radius along
  coordinate axis {i} is {arad[i]}.
  
  Namely, {\RE(arad)}  consists of all vectors {v} of {\RR^n}
  such that {v[i]} is zero where {arad[i]} is zero, and the sum of the
  squares of {v[i]/arad[i]}, over all {i} where {arad[i]} is positive, is at
  most 1.
  
  GENERIC ELLIPSOID
  
  Let {U} be a {d} by {n} matrix with orthonormal rows, and
  {urad[0..d-1]} be a list of positive numbers. We denote by
  {\RF(U,urad)} the {d}-dimensional ellipsoid in {\RR^n} whose center is
  the origin, whose main axes are parallel to the rows of {U}, and whose
  radius along row {k} of {U} is {urad[k]}.
  
  That is, {\RF(U,urad)} consists of all vectors {v} in the row span of {U}
  such that {v*U'*D*U*v' <= 1}, where {D} is the {d} by {d} diagonal matrix
  such that {D[k,k] = 1/(urad[k])^2}.
  
  In particular, if {d==n} and {U} is the identity matrix, then
  {\RF(U,urad)==\RE(urad)}. */

void rmxn_ellipsoid_pierce
  ( uint32_t n,
    double arad[],
    uint32_t d, 
    double S[],
    double U[], 
    double urad[]
  );
  /* Assumes that {arad} is an {n}-vector and {S} is a {d} by {n} matrix,
    stored by rows, with orthonormal rows.
    
    Computes matrix {U} and the vector {urad} such that {\RF(U,urad)} is
    the {d}-dimensional ellipsoid that results from the intersection of
    {\RE(arad)} the linear subspace {<S>} of {\RR^n} generated by the rows of
    {S}.  The argument {U} must be a {d} by {n} matrix, stored by rows,
    and {urad} must be a {d}-vector.
    
    That is, {\RF(U,urad)} will be the set of all vectors {u*S} such that 
    {u\in\RE(urad)}.
    
    The radii {arad[0..n-1]} must be non-negative.  For every 
    radius {arad[i]} that is zero, column {i} of {S} must be zero;
    that is the subspace {<S>} must be perpendicular to coordinate axis {i}. */
    
void rmxn_ellipsoid_cut
  ( uint32_t n,
    double arad[],
    uint32_t m, 
    double C[],
    uint32_t d,
    double U[], 
    double urad[]
  );
  /* Assumes that {arad} is an {n}-vector with non-negative elements, 
    and {C} is an {m} by {n} matrix with orthonormal rows.
    
    Computes {U} and {urad} such that {\RF(U,urad)} is the
    {d}-dimensional ellipsoid that results from the intersection of
    {\RE(arad)} and the linear subspace {<C*>} that is the orthogonal
    complement of the span {<C>} of the rows of {C}; where {d=n-m}. The
    argument {U} must be a {d} by {n} matrix, stored by rows, and {urad}
    must be a {d}-vector.
    
    That is, {\RF(U,urad)} will be the set of every vector {v}
    in {\RE(arad)} such that {v*C' = (0,..0)}. 
    
    The radii {arad[0..n-1]} must be non-negative.  For every 
    radius {arad[i]} that is zero, the coordinate axis {i} must be in 
    the supspace {<C>}. */

void rmxn_ellipsoid_normalize_constraints
  ( uint32_t n,
    double arad[],
    uint32_t q,
    double A[],
    bool_t verbose, 
    uint32_t *m_P,
    double **C_P
  );
  /* Converts a set of possibly redundant and/or incomplete linear constraints
    {A} into a set of orthonormal constraints {C} as expected by
    {rmxn_ellipsoid_constrain}.
    
    The parameter {A} must be a {q} by {n} matrix, stored by rows, which
    is interpreted as a set of {q} linear constraints {v*A' = (0,..0)}.
    The procedure implicitly appends to the rows of {A} the equations
    {v[j] == 0} for every axis {j} where {arad[j] == 0}.
    
    The procedure then applies Gram-Schmidt ortonormalization to the
    rows of {A}, eliminating any rows that are appear to be redundant.
    The resulting set of {m} orthonormal vectors is packed as a newly
    allocated {m} by {n} matrix {C}, stored by rows; and the number {m}
    and the matrix {C} are returned in {*m_P} and {*C_P}. 
    
    If {verbose} is true, prints diagnostics to {stderr}. */

#endif