function [U, rk] =  lsq_fit(nt,nx,nf,gen_data_point,verbose)

%nt - Number of data points to generate.
%nx - Number of independent variables (argument coordinates per data point). 
%nf - Number of dependent variables (functions samples per data point). 
%lsq_gen_data_point_t *gen_data_point,

%U - Fitted linear transformation matrix. 
%?bool_t verbose

%Parameter {gen_data_point} is a function with parameters (i, nx, nf);
%It is a client-implemented procedure that yields one /data
%\point/ for the purpose of least squares fitting. The data point
%consists of a set of /argument coordinates/ (independent variables)
%{Xi(1..nx}, a corresponding set of /function samples/ (dependent
%variables) {Fi(1..nf)}, and a scalar /reliability weight/ {Wi}.

%The parameter {i} is the sequential index of the data point,
%starting from 0. The procedure should return into {Xi(1..nx)} the
%argument coordinates of the data point, and into {Fi(1..nf)} the
%corresponding function samples. The procedure should also return
%the (non-negative) weight {Wi} of the data point.

%The function {lsq_fit} finds the linear function {s} from {R^nx} to {R^nf} 
%that best approximates some function {f} from {R^nx} to {R^nf}, over a 
%given set of data points, by weighted least squares.

%The procedure calls {gen_data_point(i,nx,Xi,nf,Fi,&Wi)}, with {i}
%varying from 0 to {nt-1}, to generate the argument coordinates
%{Xi(0..nx-1)}, the function samples {Fi(0..nf-1)}, and the reliability 
%weight {Wi} for each data point {i}. The weighted least squares method 
%assumes that the function samples {Fi[0..nf-1]} are measurements of some 
%``true'' linear function {f} from {R^nx} to {R^nf}, whose actual (unknown) 
%value at {Xi[0..nx-1]} is {Ti[0..nf-1]}, perturbed by the addition of {nf} 
%measurement errors {Ni[0..nf-1]}, with independent Gaussian distributions 
%with zero mean and the same variance {1/Wi}.

%The procedure stores into argument {U} the coefficient matrix {U} of
%the linear function {s} that maximizes the likelihood of the data points. 
%The matrix {U} should have {nx*nf} elements, interpreted as {nx} rows and 
%{nf} columns. The approximation is defined as the vector-matrix product 
%{s(x) = X U}. Namely, column {i} of {U} is the coefficient vector of the 
%linear combination of {x[0..nx-1]} that best approximates coordinate {i} 
%of the function {f}, in the sense of minimizing the sum of squared errors 
%{Wi*|Xi U(*,i) - Fi(i)|^2} over all data points data points {Xi,Fi,Wi}.

%Note that only the relative values of the weights are relevant. A
%data point with zero weight is ignored. Moreover, for any integer {N}, a 
%data point that has weight {Wi} counts the same as {N} data points with
%same arguments and function samples, but with weight {Wi/N}.

%!!! To implement: if {Wi} is set to {+INF} by {gen_data_point}, it means
%that the function {s} should interpolate {Fi} at {Xi}. !!!

%!!! Eliminate {nt} parameter and make {gen_data_point} return TRUE on EOF. !!!

%The procedure returns the rank {r} of the least-squares system. If {r} is 
%less than {nx}, it means that the {Xi} vectors generated by {gen_data_point}
%belonged to a proper linear subspace of {R^{nx}}. In that case, the
%best-fit map is not unique.

%This procedure uses storage proportional to {nx*(nx+nf)}, independently of
%the number of data points {nt}. Therefore, it may be preferrable to 
%{lsq_fit_arrays} when the data points are too numerous to fit in memory at 
%once.

%Build the Lsq fitting system {A U = B} for {f[0..nf-1]} in terms of 
%{x[0..nx-1]}: 

[A,B] = lsq_compute_matrix_and_rhs(nt, nx, nf, gen_data_point,verbose);

[U, rk] = lsq_solve_system(nx, nf, A, B, verbose);
