% takes about 30 seconds to run
%
% function cls_fast(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'.
% 
% if 'scraggle' is 1, use a scragglier formulation 
% for the 256 equations generated (use alternate 
% formulation). this may be omitted
% 
% see example of calling sequence in 'run_jsums.m'
%
function cls_fast(matfile, kvalue, scraggle)

  if (nargin < 2) 
    disp('cls_fast: try cls_fast(comparagram, k value of comparagram)')
    disp('cls_fast: comparagram .mat file must generate variable "outputmat"')
    error('cls_fast: try something like cls_fast("Jsum03.mat", 2)\n')
  end%if

  if (nargin < 3); scraggle=0; end%if
  
  global A C k b x n lnFinv finv  result
  
  % load comparagram
  disp("reading image"); fflush(1);
  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.'; % image to matrix
  k=kvalue;   % <======  SET 'k' HERE!!!
  n=size(C,1);
  A=zeros(n,n);
  lnFinv=zeros(n,1);
  b=zeros(n,1);
  
  disp("creating system of equations"); fflush(1);
  for i=1:n
    for j=1:n
      % grouping the 65536 equations together in contiguous groups of 256
      if (~scraggle)
        A(i, i) = A(i, i) + C(i,j);
        A(i, j) = A(i, j) - C(i,j);
        b(i)    = b(i)    + C(i,j);
      else
      % grouping the 65536 equations together in interlaced groups of 256
        % this looks very scraggly (same eqns, different form)
        A(j, i) = A(j, i) + C(i,j);
        A(j, j) = A(j, j) - C(i,j);
        b(j)    = b(j)    + C(i,j);
      end%if(scraggle)
    end
    fprintf("%d\r",i); % show it
    fflush(1);
  end
  b=b.*k;
  fprintf("\n");
  
  % now A*lnFinv=b, so solve it
  disp("solving..."); fflush(1);
  lnFinv=A\b;
  finv = e.^(lnFinv);

  % save solution  F^-1
  if (~scraggle)
    suffix = '.cls_fast.mat';
  else
    suffix = '.cls_scraggle.mat';
  end%if
  finvname = [basematfile suffix];
  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);
  if (~scraggle)
    suffix = '_fast.pnm';
  else
    suffix = '_scraggle.pnm';
  end%if
  ppmname = [basematfile suffix];
  fprintf(1, 'writing out %s\n', ppmname); fflush(1);
  saveimage(ppmname, result, "ppm");

end%function