% must have LOADPATH = [ LOADPATH ':~/comparasolve-1.3' ':' ] ; 

function verify_acs(matfile, kvalue)

  global a c s k C Csave D;

   if (nargin < 2)
      disp('must have LOADPATH = [ LOADPATH ":~/comparasolve-1.3" ":" ]'); 
      disp('use: verify_acs(''../comparunroll/Jsum06.mat'',4)');
      error("verify_acs: try something like verify_acs('Jsum03.mat', 2)\n");
   end%if

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

  % show 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");


  % store the original data, compute max ridge D
  Csave=C;
  D=C;
  mc=max(C);for i=1:size(D,1); D(i,:)=D(i,:)-mc; end;
  [I J]=find(D==0);
  D=zeros(256,256);
  for i=1:size(I,1); D(I(i),J(i))=C(I(i),J(i)); end;

  % now normalize weights by row/column
  Csave=C;
  Csave=kruithoff(Csave);
  D=kruithoff(D)*10000;
  %D(2     :256-64,256-64:256-3) = D(2     :256-64,256-64:256-3) * 64;% top edge
  D(2     :256-64,256-15 :256   ) = D(2    :256-64,256-15:256  ) * 1000;% top edge

%%%%%%%% FROM HERE ON a,c,k,s ARE log(a*a), c*c, k*k, s*s
% i.e. get model parameters by setting:
% > a=log(a*a); c=c*c; k=k*k; s=s*s;

  % starting from convergent ones
  [a c]=trial_ac;
  a=(e^(a));
  c=sqrt(c);
  k=sqrt(kvalue);
  s=sqrt(2);
  x = [a c s]';

  % solve for ridge parameters
  C=D;
  bfplotmat(C, [0], log(a*a), (c)^2,(k)^2, s*s); 
  x=bfgs_hess('comparasolve_f_k', 'comparasolve_d2f_k', x);
  a=x(1); c=x(2);  s=x(3); % k fixed
  bfplotmat(C, [0], log(a*a), (c)^2,(k)^2, s*s); 

  % solve original data parameters
  C=Csave+D; %PLUS RIDGE!!
  bfplotmat(C, [0], log(a*a), (c)^2,(k)^2, s*s);
  x=bfgs_hess('comparasolve_f_k', 'comparasolve_d2f_k', x);
  a=x(1); c=x(2);  s=x(3); % k fixed
  bfplotmat(Csave, [0], log(a*a), (c)^2,(k)^2, s*s);

%%%%%%%% FROM HERE ON a,c,k,s ARE REALLY a,c,k,s
% i.e. got model parameters by setting:
a=log(a*a); c=c*c; k=k*k; s=s*s;

  % save result as ppm (plot of f)
  result = zeros(256);
      
   % plot it
   ka = k^a;
   pf = (1:256)/256; % this is x, m, or f
   num = s^(1/c)*pf.^(1/c) * ka;
   den = pf.^(1/c) * (ka - 1) + s^(1/c);
   pg = (num./den).^c;    % this is y, n, or g

   for i=1:256
     if((pg(i)>2/256) && (pg(i)<=1) && (pf(i)>1/256) && (pf(i)<=1))
       result(pf(i)*256,pg(i)*256) = 255;
     end%if
     if(pg(i)>=1)
       result(pf(i)*256,256) = 255;
     end%if
   end%for

  result= 255*result./max(max(result));
  imshow(result);
  suffix = '_acs.pnm';
  ppmname = [basematfile suffix];
  fprintf(1, 'writing out %s\n', ppmname); fflush(1);
  %%%saveimage(ppmname, result, "ppm");
 
  % save a,c,s parameters
  suffix = '_acs.txt';
  acsname = [basematfile suffix];
  fout = fopen(acsname, 'w');
  % report results in 'saturation' model parameters... 
  fprintf(fout, '% 20.20g \t# "a"\n',a);
  fprintf(fout, '% 20.20g \t# "c"\n',c);
  fprintf(fout, '% 20.20g \t# "k"\n',k);
  fprintf(fout, '% 20.20g \t# "s"\n',s);
  fclose(fout);

  % save finv values 
  suffix = '_acs_finv.txt';
  acsname = [basematfile suffix];
  fout = fopen(acsname, 'w');
  pf = (1:256)/256; % this is x, m, or f
  % from saturation model
  finv = exp((-log(-(-1.0+(pf./s).^(1/c))).*c+log(pf./s))./c./a);
  for i=1:max(size(finv))
    fprintf(fout, '% 20.20g\n',finv(i));
  end
  fclose(fout);
end%function

