% takes about 5 minutes on gigahertz p3
%
% function cls(matfile, kvalue, [scraggle])
%
% run least squares solution for 256 points of
% finverse directly from comparagram saved in
% variable 'outputmat' in file 'matfile', with
% known 'k' value, 'kvalue'.
%
% uses slow method, generating 65536 equations in 
% 256 variables.
%
% see example of calling sequence in 'run_jsums.m'

function cls(matfile, kvalue)

   if (nargin < 2)
      error("cls_fast: try something like cls('Jsum03.mat', 2)\n");
   end%if

  global A C k b x n lnFinv finv  result
  % C is comparagram
  % b is rhs
  % x is unknown
  % n is the size of C

  % load comparagram
  disp("reading image"); fflush(1);
  %load("Jsum03.mat")
  load('-force',matfile)
  C=outputmat;
  basematfile = matfile(1:(max(size(matfile))-4)); % chop '.mat'


  % save decent-looking pnm version
  %Csave = C.^(1/3);
  %Csave = 255.0*(Csave./max(max(Csave)));
  %imshow(Csave);
  %ppmname = [basematfile '.pnm'];
  %fprintf(1, 'writing out %s\n', ppmname); fflush(1);
  %saveimage([basematfile '.pnm'], Csave, "ppm");


  % parameter initialization
  C=C.';
  n=size(C,1); % assumes comparagram C is square
  A=zeros(n*n,n);
  lnFinv=zeros(n,1)*NaN; % just to declare
  k=kvalue;
  b=zeros(n*n,1);
  
  disp("creating system of equations"); fflush(1);
  for i=1:n
    for j=1:n
  %    if (i != j)
         A(i+n*(j-1), i) = A(i+n*(j-1), i) + C(i,j);
         A(i+n*(j-1), j) = A(i+n*(j-1), j) - C(i,j);
         b(i+n*(j-1))    = b(i+n*(j-1))    + C(i,j);
  %    end
    end
    fprintf("%d\r",i); % show it
    fflush(1);
  end
  b=b.*k;
  fprintf("\n");

  % now A*lnFinv=b
  %solve it
  disp("solving..."); fflush(1);
  lnFinv=A\b;
  finv = e.^(lnFinv);

  % save solution  F^-1
  finvname = [basematfile '.cls.mat'];
  fprintf(1, 'writing out %s\n', finvname); fflush(1);
  save(finvname, 'lnFinv', 'finv', 'k');

  % compute f(q) vs f(kq)
  clshow(k);

  % save f(q) vs f(kq) plot
  result= 255*result./max(max(result));
  imshow(result);
  suffix = '_leastsq.pnm';
  ppmname = [basematfile suffix];
  fprintf(1, 'writing out %s\n', ppmname); fflush(1);
  saveimage(ppmname, result, "ppm");
end%function

