function [U, rk] = lsq_array_fit(nt, nx, nf, X, F, W, verbose)
%int lsq_array_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.

%Assumes that the data points are given as two matrices {X,F} and a vector
%{W} (instead of a generating procedure {gen_data_point}, as in {lsq_fit}
%from {lsq.h}). The matrices {X} and {F} should have {nt} rows and
%{nx} or {nf} columns, respectively, linearized by rows; and the vector
%{W} should have {nt} elements, interpreted as a column. Each data point
%is therefore {Xi[j] = X[i*nx + j]} for {j} in {0..nx-1}, 
%{Fi[k] = F[i*nf + k]} for {k} in {0..nf-1}, and {Wi = W[i]}. If {W} is  
%null, assumes unit weights for all arguments.

%The procedure stores into {U} the coefficient matrix of the
%fitted linear approximation {s}.
%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 U}. Namely, column {j} of {U} is the coefficient vector of the 
%linear combination of {Xi[0..nx-1]} that best approximates coordinate {j} 
%of the function {Fi}, in the sense of minimizing the sum of squared errors 
%{Wi*|X[i,*] U[*,j] - F[i,j]|^2} over all data points {i} in {0..nt-1}.

%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}.

%Returns the rank {r} of the least-squares system. If {r} is less
%than {nx}, it means that the rows of {X} span a proper linear subspace of
%{R^{nx}}. In that case, the best-fit map is not unique.

%This method uses internal working storage proportional {nx*(nx+nf)},
%like {lsq_fit}.


%nt - Number of data points to generate. 
%nx - Number of independent variables (argument coordinates per data point). 
%nf - Number of dependent variables (function samples per data point). 
%X - Argument coordinates of data points ({nt} by {nx}). 
%F - Corresponding function samples ({nt} by {nf}. 
%W - Corresponding weights ({nt} elements). 
%U - (OUT) Fitted linear transformation matrix ({nx} by {nt}). 
%bool_t verbose

A = lsq_array_compute_matrix(nt, nx, X, W);
B = lsq_array_compute_rhs(nt, nx, nf, X, F, W);
[U, rk] = lsq_solve_system(nx, nf, A, B, verbose);

