#ifndef minn_constr_H
#define minn_constr_H

/* General tools and concepts for minimization with linear constraints. */
/* Last edited on 2023-03-26 09:12:16 by stolfi */ 

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

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

#include <minn.h>

/* 
  THE SEARCH ELLIPSOID
  
  The functions in this interface search the minimum of a function in a
  {d}-dimensional /search ellipsoid/ {\RF} contained in {\RR^n}, for
  some {d} and {n}. The search ellipsoid is centered at the origin of
  {\RR^n}, and is defined by an /axial radius vector/ {arad[0..n-1]} and a
  /constraint matrix/ {C} with {p} rows and {n} columns.

  THE BASE ELLIPSOID
  
  The radius vector defines an {n}-dimensional /base ellipsoid/ {\RE},
  also centered at the origin, whose main axes are aligned with the
  coordinate axes of {\RR^n}. This ellipsoid has radius {arad[i]} along
  coordinate axis {i}. Namely, it 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. The radii {arad[i]} must not be negative.
  
  THE CONSTRAINT SUBSPACE
  
  The rows of the constraint matrix {C} must be an orthonormal basis.
  The /constraint subspace/ {\RB} is the set of vectors {v} of {\RR^n}
  such that {v*C'} is the zero {p}-vector. That is, each row {c[j]} of
  {C} represents the linear constraint {v*c[j]' = 0}.  The constraints 
  must be such that they imply {v[j] == 0} for every axis {j}
  such that {arad[j]} is zero.
  
  The search ellipsoid {\CF} is then the intersection of the linear
  subspace {\CB} with the basic ellipsoid {\CE}.
  
  The dimension {d} of the space {\CB} and of the ellipsoid {\CF}
  is therefore {n-p}. */

void minn_constr_compute_search_ellipsoid
  ( int32_t n,
    double arad[],
    int32_t p, 
    double C[],
    int32_t d,
    double U[], 
    double urad[]
  );
  /* Assumes that {C} is a {p} by {n} constraint matrix, {U} is a {d} by {n} matrix,
    where {d}, both stored by rows; and {urad} is a {d}-vector.  Requires
    {p+d == n}.
    
    Stores into the rows of {U} a /search basis/, a set of {d}
    orthonormal that are the main directions of the search ellipsoid
    {\CF}, and into {urad[0..d-1]} the correspondng major semiaxes of
    that ellipsoid.  The {d} rows of {U} will be an orthonormal basis for
    the space {\RB} of the vectors that satify the constraints {v*C' =
    (0,..0)}. */

void minn_constr_normalize_constraints
  ( int32_t n,
    double arad[],
    int32_t q,
    double A[],
    bool_t verbose, 
    int32_t *p_P,
    double **C_P
  );
  /* Converts a set of possibly redundant and/or incomplete constraints
    {A} into a set of orthonormal constraints {C} as expected by
    {minn_constr_compute_search_ellipsoid}.
    
    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 {p} orthonormal vectors is packed as a newly
    allocated {p} by {n} matrix {C}, stored by rows; and the number {p}
    and the matrix {C} are returned in {p_P} and {*C_P}. 
    
    If {verbose} is true, prints diagnostics to {stderr}. */

#endif