% takes about 5 minutes on gigahertz p3
%
% function Finv=cls(J, k)
%
% run least squares solution for 256 points of
% finverse directly from comparagram J and known k.
%
% uses slow method, generating 65536 equations in 
% 256 unknowns.
%
% see usage example by running 'runme.m'


function Finv=cls(J, kvalue);

  % n=size(J,1); % assumes comparagram J is square
  [m,n]=size(J);
  if m~=n
    error('cls_fast: input comparagram must be square\n');
  end%if


  if (nargin < 2)
     error("use: try something like Finv=cls(Jsum03,2)\n");
  end%if
  if (nargout < 1)
     error("use: try something like Finv=cls(Jsum03,2)\n");
  end%if

  % J is the comparagram
  % b is rhs
  % n is the size of C

  % load comparagram
  disp("reading image"); fflush(1);
  %load("Jsum03.mat")

  % parameter initialization
  J=J.';
  A=zeros(m*n,n);
  Finv=zeros(n,1); % zeros(n,1)*NaN; % just to declare as 0xDEADBEEF
  k=kvalue;
  K=log(k)/log(2);
  b=zeros(m*n,1); % zeros(n*n,1)*NaN; % as 0xDEADBEEF
  
  disp("creating system of equations"); fflush(1);
  for i=1:n
    for j=1:n
         A(i+n*(j-1), i) = A(i+n*(j-1), i) + J(i,j);
         A(i+n*(j-1), j) = A(i+n*(j-1), j) - J(i,j);
         b(i+n*(j-1))    = b(i+n*(j-1))    + K*J(i,j);
    end
    fprintf("%d\r",i); % show it
    fflush(1);
  end
  fprintf("\n");

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