function occqp_cmd(infile,outfile,solve)
% OccQP_cmd
%
% Command-line user interface to the Occupancy program.
%
% Usage:  occqp_cmd(infile)
%   Infile must be a character variable giving the name of a Matlab script
%   M-file declaring values for the Matlab variables Atoms, Sites, Qobs,
%   Cobs, C, Coord, dA. Optionally, the script may also be declared values
%   for any of the Matlab variables x, Q, R, Z, Caggr, aIonic, fixed, full,
%   W_Q, W_V, W_D, W_C, W_I. See the documentation for further details.
%
%     NOTE: The character string in Infile must end in '.m'.
%
%   The variable aIonic is a target value for the (positive) ionic sum that
%   must be met in the ionic balance equation. This target value is usually
%   known in advance and specifying it greatly improves the results.
%
%   Infile is overwritten by an output file giving the data, all model
%   parameter values used, the solution occupancy x, and the model
%   residuals.
%
%   If no value for x is specified in Infile, then optimization is
%   attempted. If x is specified in Infile, then optimization is not
%   attempted and the output uses the given x to calculate the residuals.
%
% Optional usage:
%
%   occqp_cmd(infile,outfile)
%     Outfile must be a valid filename for a file that will be overwritten
%     with the results instead of overwriting Infile. If outfile is empty
%     or missing, then Infile will be overwritten.
%
%   occqp_cmd(infile,outfile,solve)
%     Solve is a logical value indicating whether the optimizer should be
%     run. If Solve is true, then the optimizer is used regardless of
%     whether Infile declares values for x. If Solve is false, then the
%     optimizer is not called and the residuals are based on the given x
%     if it is specified in Infile. If Solve is false and Infile does not
%     provide x, then x is taken to be evenly distributed over all cation
%     sites.
%

% Copyright 1999-2019 by Stephen E. Wright
% Department of Statistics
% Miami University, Oxford, OH.
% Version number 1.00; 05 July December 2019

% Note: "error" handling via print-and-return rather than Matlab's error
% utility. This is more graceful for users who unfamiliar with Matlab.

if nargin<1
    fprintf('Error: requires names of an input data M-file.')
    return
elseif nargin<2
    outfile = infile;
end 

% obtain the user's input
[data,errorflag] = get_data(infile);
if errorflag
    return
end

if isempty(data.x) % use uniform distribution for initial guess
    data.x=ones(data.n,data.m)/data.n;
    if nargin<3
        solve = true;
    end
elseif nargin<3
    solve = false;
end

% solve the occupancy optimization problem
if solve
    % build weights for solver
    full=reshape(deal(data.full),data.m,1);
    filled=zeros(size(full));
    filled(full)=inf*ones(sum(full),1);
    W_res=reshape(deal(data.W_res)',3*data.m,1);
    weights=[filled;W_res;deal(data.W_C);data.W_I];

    % call the solver
    [X,errstr]=... %[X,errstr,Verts]
        occup(data.Q,data.Qobs,data.Z,data.V,data.C,data.Cobs,...
        data.aIonic,data.D,data.dbar,weights,data.fixed,data.x,...
        data.Caggr);
    if ~isempty(errstr)
        fprintf('Error!\n');
        fprintf(errstr)
        return
    else
        data.x=X;
    end
end

% create output file
data=residual_calc(data);
if nargin<2||isempty(outfile)
    outfile = infile;
end
write_Mfile(data,outfile)
end % occqp_cmd


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The remaining code is divided into two main sections:
% 1. Interface: 
%    a. get_data, read_Mfile - input data, parameters, occupancy (if given)
%    b. parameter_lookup - get default parameters
%    c. residual_calc - calculate residuals from current occupancy
%    d. write_Mfile - output data, parameters, occupancy and residuals
% 2. Optimization: 
%    a. occup - construct the mathematical model
%    b. use_LHnnls, LHnnls, h1, h2 - perform numerical optimization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Interface code
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function [data,errorflag] = get_data(infile)
% place data in a structure to simplify passing among functions

% obtain the data as individual variables
[Name,Atoms,Sites,...                % required
    Qobs,Cobs,C,Coord,dA,...         % required
    Caggr,W_res,W_C,full,x,fixed,... % optional
    Q,Z,R,...                        % optional
    W_I,aIonic,...                   % optional
    V,D,dbar,...                     % calculated
    errorflag]...
    = read_Mfile(infile);
if errorflag, return, end

m=length(Sites); % number of sites
n=length(Atoms); % number of species (atoms)
if isempty(fixed), fixed=false(n,m); end
if isempty(full),  full =false(1,m); end

data=struct(...
    'm',m,'n',n,...
    'name',{Name},'atoms',{Atoms},'sites',{Sites},...
    'Qobs',{Qobs},'Cobs',{Cobs},...
    'C',{C},'Coord',{Coord},'dA',{dA},...
    'Q',{Q}, 'Z',{Z}, 'R',{R},...
    'V',{V}, 'D',{D}, 'dbar',{dbar},...
    'Caggr',{Caggr},...
    'W_res',{W_res},'W_C',{W_C},...
    'W_I',{W_I},'aIonic',{aIonic},'cIonic',{[]},...
    'full',{full},'vac',{zeros(1,m)},...
    'fixed',{fixed},'x',{x},...
    'Cx',{zeros(size(x))},'chem_calc',{[]},...
    'Cvac',{[]},'cc_vac',{[]},...
    'res',{[]},'r_calc',{[]},'r_obs',{[]});
if isempty(aIonic), data.aIonic=0; end
end % OF FUNCTION get_data


function [Name,Atoms,Sites,...           % required
        Qobs,Cobs,C,Coord,dA,...         % required
        Caggr,W_res,W_C,full,x,fixed,... % optional
        Q,Z,R,...                        % optional
        W_I,aIonic,...                   % optional
        V,D,dbar,...                     % calculated
        errorflag]...
        = read_Mfile(datafile1093784746) % nonsense name to avoid conflicts
% "Reads" the data by running the M-file specified in the input argument

% Make sure these are all defined somehow, in case they're not in the file
Name=[]; Atoms=[]; Sites=[];
Qobs=[]; Cobs=[]; C=[]; Coord=[]; dA=[];
Caggr=[]; W_res=[]; W_Q=[]; W_V=[]; W_D=[]; W_C=[]; full=[]; x=[]; fixed=[];
Q=[]; Z=[]; R=[];
V=[]; D=[]; dbar=[];
W_I=[];

% Attempt to run the data M-file
if ~ischar(datafile1093784746)
    fprintf('Error: Data file name must be a character variable');
    return
end
fprintf('Reading data file: %s\n',datafile1093784746);
curdir1093784746=pwd;
[pname1093784746,fname1093784746]=fileparts(datafile1093784746);
if ~isempty(pname1093784746), cd(pname1093784746), end
try
    eval(fname1093784746);
    errorflag = true; % assume the worst, but set it after running file
    cd(curdir1093784746)
catch
    fprintf(['Error! Data file is not a valid Matlab'...
        ' command file!\n To see the nature of the error'...
        ' use the command\n    >>  %s \n at the Matlab'...
        ' command prompt.\n\n'],datafile1093784746);
    errorflag = true;
    cd(curdir1093784746)
    return
end

% Is all the data defined?
if ~exist('Atoms','var')||isempty(Atoms)
    fprintf('Error! The data file did not set up Atoms.\n');
    return
elseif ~exist('Sites','var')||isempty(Sites)
    fprintf('Error! The data file did not set up Sites.\n');
    return
elseif ~exist('Qobs','var')||isempty(Qobs)
    fprintf('Error! The data file did not set up Qobs.\n');
    return
elseif ~exist('C','var')||isempty(C)
    fprintf('Error! The data file did not set up C.\n');
    return
elseif ~exist('Cobs','var')||isempty(Cobs)
    fprintf('Error! The data file did not set up Cobs.\n');
    return
elseif ~exist('Coord','var')||isempty(Coord)
    fprintf('Error! The data file did not set up Coord.\n');
    return
elseif ~exist('dA','var')||isempty(dA)
    fprintf('Error! The data file did not set up dA.\n');
    return
end

% Reshape as rows
Atoms=Atoms(:)';
Sites=Sites(:)';
Qobs=Qobs(:)';
Cobs=Cobs(:)';
C=C(:)';

% Check dimensions
n=length(Atoms); % number of atoms (species)
m=length(Sites); % number of sites
if length(Qobs)~=m
   fprintf('Error! Qobs incorrectly dimensioned!\n');
   return
elseif length(C)~=m
   fprintf('Error! C incorrectly dimensioned!\n');
   return
elseif length(Coord)~=m
   fprintf('Error! Coord incorrectly dimensioned!\n');
   return
elseif any(size(dA)~=[sum(Coord),2])
   fprintf('Error! dA incorrectly dimensioned!\n');
   return
end

% Check weights
if exist('weights','var') % no backward compatibility
   fprintf(['\n   Warning! ''weights'' has been replaced by'...
            ' ''full'', ''W_Q'', etc---\n   using default weights.\n']);
end

% Set up defaults for site residuals
W_res=[ones(2,m);zeros(1,m)];

if ~isempty(W_Q)
   if length(W_Q(:))~=m
      fprintf(['\n   Warning: Incorrect number of weights'...
                 ' specified for scattering---using defaults.\n']);
   elseif any(W_Q<0)
      fprintf(['\n   Warning: Negative weights not allowed'...
                 '---using default weights for scattering.\n']);
   else
      W_res(1,:)=W_Q(:)';
   end
end

if ~isempty(W_V)
   if length(W_V(:))~=m
      fprintf(['\n   Warning: Incorrect number of weights'...
                 ' specified for valence---using defaults.\n']);
   elseif any(W_V<0)
      fprintf(['\n   Warning: Negative weights not allowed'...
                 '---using default weights for valence.\n']);
   else
      W_res(2,:)=W_V(:)';
   end
end

if ~isempty(W_D)
   if length(W_D(:))~=m
      fprintf(['\n   Warning: Incorrect number of weights'...
              ' specified for bond length---using defaults.\n']);
   elseif any(W_D<0)
      fprintf(['\n   Warning: Negative weights not allowed'...
                 '---using default weights for bond length.\n']);
   else
      W_res(3,:)=W_D(:)';
   end
end

if ~isempty(W_I)
   if length(W_I)~=1
      fprintf(['\n   Warning: Incorrect number of weights'...
              ' specified for ionic balance---using defaults.\n']);
   elseif (W_I<0) %#ok<BDSCI>
      fprintf(['\n   Warning: Negative weights not allowed'...
                 '---using default weight for ionic balance.\n']);
   else
      W_I=W_I(:);
   end
else
   W_I=0;
end

if ~isempty(W_C)
   if length(W_C(:))~=n
      fprintf(['\n   Warning: Incorrect number of weights'...
               ' specified for chemistry---using defaults.\n']);
      W_C=ones(n,1);
   elseif any(W_C<0)
      fprintf(['\n   Warning: Negative weights not allowed'...
                 '---using default weights for chemistry.\n']);
      W_C=ones(n,1);
   else
      W_C=W_C(:);
   end
else
   W_C=ones(n,1);
end

if ~isempty(full)&&any(size(full)~=[1,m])
   fprintf(['Error!\n Zero-one array filling vacancy is'...
            ' incorrectly dimensioned!\n']);
   return
end
full=~(~full);

% Check occupancy
if ~isempty(x)&&any(size(x)~=[n,m])
   fprintf('Error!\n Occupancy incorrectly dimensioned!\n');
   return
end

% Check occupancy weights
if exist('xfixed','var') % backward compatibility
   fprintf(['\n   Warning: ''xfixed'' has been replaced by'...
            '''fixed''---\n   future versions will'...
            ' ignore ''weights''.\n']);
   fixed=xfixed;
end
if ~isempty(fixed)&&any(size(fixed)~=[n,m])
   fprintf(['Error!\n Zero-one array fixing occupancy is'...
            ' incorrectly dimensioned!\n']);
   return
end
fixed=~(~fixed);

% Periodic table
if ~isempty(Q)
   if length(Q)~=n
      fprintf(['Error!\n User-defined formal scattering incorrectly'...
              '  dimensioned!\n']);
      return
   end
   Q=Q(:)';
else
   Q=parameter_lookup({'ptQ'},Atoms);
   if any(isnan(Q(:)))
      fprintf(['Error!\n Formal charge not defined for the'...
                 ' following cations:\n']);
      iC=find(isnan(Q));
      for i=1:length(iC), fprintf('   %s\n',Atoms{iC(i)}); end
      return
   end
end

if ~isempty(Z)
   if length(Z)~=n
      fprintf(['Error!\n User-defined formal valence incorrectly'...
              ' dimensioned!\n']);
      return
   end
   Z=Z(:)';
else
   Z=parameter_lookup({'ptZ'},Atoms);
   if any(isnan(Z(:)))
      fprintf(['Error!\n'...
               ' Formal valence not defined for the'...
                    ' following cations:\n']);
      iC=find(isnan(Z));
      for i=1:length(iC), fprintf('   %s\n',Atoms{iC(i)}); end
      return
   end
end

% Bond-valence parameters
if ~isempty(R)
   if (size(R,2)~=n+1)
      fprintf(['Error!\n'...
              '  User-defined bond-valence parameter array'...
              ' incorrectly dimensioned!\n']);
      return
   end
   rR=zeros(sum(Coord),n);
   for i=1:sum(Coord)
      iA=find(startsWith(R(:,end),dA{i,2}),1,'first');
      if ~isempty(iA), rR(i,:)=[R{iA,1:n}];
      else
         fprintf(['Error!\n'...
                  ' User-defined bond-valence parameter'...
                   ' values omitted for anion']);
         fprintf('  %s\n',dA{i,2});
         return
      end
   end
else
   rR=parameter_lookup(dA(:,2),Atoms);
   [anions,akey]=sort(dA(:,2));
   aselect=false(size(anions));
   i0=1;
   iinf=sum(Coord);
   while i0<=iinf
      aselect(i0)=1;
      i0=find(strcmp(anions,anions(i0)),1,'last')+1;
   end
   R=[ num2cell(rR(akey(aselect),:)), anions(aselect) ];
   if any(isnan([R{:,1:n}]))
      fprintf(['Error!\n'...
               ' Bond-valence parameter not defined for'...
                 ' the following anion/cation pairs:\n']);
      [iA,iC]=find(isnan([R{1:n,:}]));
      for i=1:length(iA)
         fprintf('   %s/%s\n',R{iA(i),2},Atoms{iC(i)});
      end
      return
   end
end

% Compute V and average bond lengths from Coord, Atoms, dA

b=0.37; % universal parameter
V=zeros(m,n);
D=zeros(m,n);
dbar=zeros(1,m);
d=[dA{:,1}]';

i0=1; i1=0;
for j=1:m
  ncoord=Coord(j);
  i1=i1+ncoord;

  V(j,:)=sum(exp( (rR(i0:i1,:)-d(i0:i1)*ones(1,n)) / b));
  D(j,:)=-b*log(Z/ncoord)+sum(rR(i0:i1,:))/ncoord;
  dbar(j)=sum(d(i0:i1))/ncoord;

  i0=i1+1;
end

% Aggregate chemistry
if ~isempty(Caggr)&&(length(Caggr)~=n)
   fprintf('Error!\n Caggr incorrectly dimensioned!')
   return
end
Caggr=fix(Caggr);

% anion value for ionic balance
if exist('aIonic','var')
   if length(aIonic)>1 %#ok<NODEF>
      fprintf('Error!\n aIonic incorrectly dimensioned!')
      return
   end
else
   W_I=0;
   aIonic=[];
end

errorflag = false; % all is well
end % OF FUNCTION read_Mfile


function parm=parameter_lookup(whichP,Cations)
% Looks up parameters whichP for the indicated Cations.
% Both whichP and Cations should be 1D cell-arrays of strings.
% The returned values are arranged in a 2D array of numbers, with
% whichP labelling the rows and Cations labelling the columns.
%
% NaN's signify 'Not Available'.
%

% Table:
%  * formal charge and valence from periodic table
%  * bond-valence parameters from Brese & O'Keeffe, Acta Cryst.
%                                            (1991),B47,192-197
P=            { 'ptQ'  'ptZ'  'O'    'F'   'Cl'   'S' };
C={'Ac'          89     3    2.24    2.13  2.63   nan
   'Ag'          47     1    1.805   1.80  2.09   2.15
   'Al'          13     3    1.651   1.545 2.03   2.13
   'Am'          95     3    2.11    2.00  2.48   nan
   'As(III)'     33     3    1.789   1.70  2.16   2.26
   'As(V)'       33     5    1.767   1.62  2.14   2.26
   'Au'          79     3    1.833   1.81  2.17   2.03
   'B'            5     3    1.371   1.31  1.74   1.82
   'Ba'          56     2    2.29    2.19  2.69   2.77
   'Be'           4     2    1.381   1.28  1.76   1.83
   'Bi(III)'     83     3    2.09    1.99  2.48   2.55
   'Bi(V)'       83     5    2.06    1.97  2.44   2.55
   'Bk'          97     3    2.08    1.96  2.46   nan
   'Br'          35     7    1.81    1.72  2.19   nan
   'C'            6     4    1.39    1.32  1.76   1.82
   'Ca'          20     2    1.967   1.842 2.37   2.45
   'Cd'          48     2    1.904   1.811 2.23   2.29
   'Ce(III)'     58     3    2.151   2.036 2.52   2.62
   'Ce(IV)'      58     4    2.028   1.995 2.41   2.62
   'Cf'          98     3    2.07    1.95  2.45   nan
   'Cl'          17     7    1.632   1.55  2.00   nan
   'Cm'          96     3    2.23    2.12  2.62   nan
   'Co(II)'      27     2    1.692   1.64  2.01   2.06
   'Co(III)'     27     3    1.70    1.62  2.05   2.06
   'Cr(II)'      24     2    1.73    1.67  2.09   2.18
   'Cr(III)'     24     3    1.724   1.64  2.08   2.18
   'Cr(VI)'      24     6    1.794   1.74  2.12   2.18
   'Cs'          55     1    2.42    2.33  2.79   2.89
   'Cu(I)'       29     1    1.593   1.6   1.85   1.86
   'Cu(II)'      29     2    1.679   1.60  2.00   1.86
   'Dy'          66     3    2.036   1.922 2.41   2.47
   'Er'          68     3    2.010   1.906 2.39   2.46
   'Eu(II)'      63     2    2.147   2.04  2.53   2.53
   'Eu(III)'     63     3    2.076   1.961 2.455  2.53
   'Fe(II)'      26     2    1.734   1.65  2.06   2.16
   'Fe(III)'     26     3    1.759   1.67  2.09   2.16
   'Ga'          31     3    1.730   1.62  2.07   2.17
   'Gd'          64     3    2.065   1.95  2.445  2.53
   'Ge'          32     4    1.748   1.66  2.14   2.22
   'H'            1     1    0.95    0.92  1.28   1.35
   'Hf'          72     4    1.923   1.85  2.30   2.39
   'Hg(I)'       80     1    1.90    1.81  2.28   2.32
   'Hg(II)'      80     2    1.93    1.90  2.25   2.32
   'Ho'          67     3    2.023   1.908 2.401  2.48
   'I(V)'        53     5    2.00    1.90  2.38   nan
   'I(VII)'      53     7    1.93    1.83  2.31   nan
   'In'          49     3    1.902   1.79  2.28   2.36
   'Ir'          77     5    1.916   1.82  2.30   2.38
   'K'           19     1    2.13    1.99  2.52   2.59
   'La'          57     3    2.172   2.057 2.545  2.64
   'Li'           3     1    1.466   1.360 1.91   1.94
   'Lu'          71     3    1.971   1.876 2.361  2.43
   'Mg'          12     2    1.693   1.581 2.08   2.18
   'Mn(II)'      25     2    1.790   1.698 2.13   2.20
   'Mn(III)'     25     3    1.760   1.66  2.14   2.20
   'Mn(IV)'      25     4    1.753   1.71  2.13   2.20
   'Mn(VII)'     25     7    1.79    1.72  2.17   2.20
   'Mo'          42     6    1.907   1.81  2.28   2.35
   'N(III)'       7     3    1.361   1.37  1.75   nan
   'N(V)'         7     5    1.432   1.36  1.80   nan
   'Na'          11     1    1.80    1.677 2.15   2.28
   'Nb'          41     5    1.911   1.87  2.27   2.37
   'Nd'          60     3    2.117   2.008 2.492  2.59
   'Ni'          28     2    1.654   1.599 2.02   2.04
   'Os'          76     4    1.811   1.72  2.19   nan
   'P'           15     5    1.604   1.521 1.99   2.11
   'Pb(II)'      82     2    2.112   2.03  2.53   2.55
   'Pb(IV)'      82     4    2.042   1.94  2.43   2.55
   'Pd'          46     2    1.792   1.74  2.05   2.10
   'Pr'          59     3    2.135   2.022 2.50   2.60
   'Pt(II)'      78     2    1.768   1.68  2.05   2.08
   'Pt(IV)'      78     4    1.879   1.759 2.17   2.08
   'Pu'          94     3    2.11    2.00  2.48   nan
   'Rb'          37     1    2.26    2.16  2.65   2.70
   'Re'          75     7    1.97    1.86  2.23   2.37
   'Rh'          45     3    1.791   1.71  2.17   2.15
   'Ru'          44     4    1.834   1.74  2.21   2.16
   'S(IV)'       16     4    1.644   1.60  2.02   2.07
   'S(VI)'       16     6    1.624   1.56  2.03   2.07
   'Sb(III)'     51     3    1.973   1.90  2.35   2.45
   'Sb(V)'       51     5    1.942   1.80  2.30   2.45
   'Sc'          21     3    1.849   1.76  2.23   2.32
   'Se(IV)'      34     4    1.811   1.73  2.22   2.25
   'Se(VI)'      34     6    1.788   1.69  2.16   2.25
   'Si'          14     4    1.624   1.58  2.03   2.13
   'Sm'          62     3    2.088   1.977 2.466  2.55
   'Sn(II)'      50     2    1.984   1.925 2.36   2.45
   'Sn(IV)'      50     4    1.905   1.84  2.28   2.45
   'Sr'          38     2    2.118   2.019 2.51   2.59
   'Ta'          73     5    1.920   1.88  2.30   2.39
   'Tb'          65     3    2.049   1.936 2.427  2.51
   'Te(IV)'      52     4    1.977   1.87  2.37   2.45
   'Te(VI)'      52     6    1.917   1.82  2.30   2.45
   'Th'          90     4    2.167   2.07  2.55   2.64
   'Ti(III)'     22     3    1.791   1.723 2.17   2.24
   'Ti(IV)'      22     4    1.815   1.76  2.19   2.24
   'Tl(I)'       81     1    2.172   2.15  2.56   2.63
   'Tl(III)'     81     3    2.003   1.88  2.32   2.63
   'Tm'          69     3    2.000   1.842 2.38   2.45
   'U(IV)'       92     4    2.112   2.034 2.48   2.56
   'U(VI)'       92     6    2.075   1.966 2.46   2.56
   'V(III)'      23     3    1.743   1.702 2.19   2.23
   'V(IV)'       23     4    1.784   1.70  2.16   2.23
   'V(V)'        23     5    1.803   1.71  2.16   2.23
   'W'           74     6    1.921   1.83  2.27   2.39
   'Y'           39     3    2.014   1.904 2.40   2.48
   'Yb'          70     3    1.985   1.875 2.371  2.43
   'Zn'          30     2    1.704   1.62  2.01   2.09
   'Zr'          40     4    1.937   1.854 2.33   2.41  };

% reshape the data
data=reshape([C{:,2:end}],[size(C,1) size(P,2)])';

% look up the parameter names
nw=length(whichP);
iw=false(nw,1);
ip=zeros(nw,1);
for i=1:nw
   m=find(strcmp(whichP(i),P));
   if ~isempty(m)
       iw(i)=1;
       ip(i)=m;
   end
end

% look up the cations
ncat=length(Cations);
icat=false(ncat,1);
ic=zeros(ncat,1);
for i=1:ncat
   m=find(strcmp(Cations(i),C(:,1)));
   if ~isempty(m)
       icat(i)=1;
       ic(i)=m;
   end
end

% extract the desired data
parm=nan(nw,ncat);
parm(iw,icat)=data(ip(iw),ic(icat));
end % OF FUNCTION parameter_lookup


function data=residual_calc(data)
% calculates residual values from current x
data.chem_calc=data.C*data.x';
data.Cx=data.x*diag(data.C); % site population
data.vac=1-sum(data.x);
data.cc_vac=data.C*data.vac'; % total vacancy
data.Cvac=data.vac.*data.C;   % vacancy population
data.r_calc=[
    data.Q*data.x
    data.Z*data.x
    diag(data.D*data.x)'
    ];
data.r_obs=[
    data.Qobs
    diag(data.V*data.x)'
    diag(data.dbar(:)*ones(1,data.n)*data.x)'
    ];
data.res=data.r_calc-data.r_obs;
data.cIonic=data.Z*data.chem_calc(:);
end % OF FUNCTION residual_calc


function write_Mfile(data,outfile)
% Writes all current data values in the file named by datafile
% get the user data

% open the file
fid=fopen(outfile,'w+t');
if fid<0
   dialog(sprintf(' Error opening file %s !',outfile))
   return
end
fprintf('Writing output file: %s\n',outfile);

if ~isempty(data.name)
   fprintf(fid,'Name = ''%s'';\n',data.name);
end

fprintf(fid,'\n%% REQUIRED INPUT PARAMETERS\n\n');

fprintf(fid,'Atoms = {');
for k=1:data.n
    fprintf(fid,' %9s',sprintf('''%s''',data.atoms{k}));
end
fprintf(fid,' };\n');

fprintf(fid,['Cobs  = [' sprintf(' %9.4f',data.Cobs) ' ];\n\n']);

fprintf(fid,'Sites = {');
for j=1:data.m
    fprintf(fid,' %9s',sprintf('''%s''',data.sites{j}));
end
fprintf(fid,' };\n');

fprintf(fid,['Qobs  = [' sprintf(' %9.4f',data.Qobs) ' ];\n']);
fprintf(fid,['C     = [' sprintf(' %9d',data.C)     ' ];\n']);
fprintf(fid,['Coord = [' sprintf(' %9d',data.Coord) ' ];\n\n']);

fprintf(fid,'dA = {');
ii=0;
for j=1:data.m
    for i=ii+(1:data.Coord(j))
       fprintf(fid,'%8.4f %4s',data.dA{i,1},...
                   sprintf('''%s''',data.dA{i,2}));
       if i==size(data.dA,1)
          fprintf(fid,' }; %%''%s''\n',data.sites{j});
       else
          fprintf(fid,'    %%''%s''\n      ',data.sites{j});
       end
    end
    ii=ii+data.Coord(j);
end

fprintf(fid,'\n%% OPTIONAL INPUT PARAMETERS\n\n');

fprintf(fid,'%%        ');
fprintf(fid,' %8s',data.sites{:});

fprintf(fid,'\nx     = [');
for k=1:data.n
    fprintf(fid,' %8.4f',data.x(k,:));
    if k==data.n
        fprintf(fid,' ]; %% %s',data.atoms{k});
    else
        fprintf(fid,'    %% %s\n         ',data.atoms{k});
    end
end

fprintf(fid,'\nfixed = [');
for k=1:data.n
    fprintf(fid,' %8d',data.fixed(k,:));
    if k==data.n
        fprintf(fid,' ]; %% %s',data.atoms{k});
    else
        fprintf(fid,'    %% %s\n         ',data.atoms{k});
    end
end

fprintf(fid,['\nfull  = [' sprintf(' %8d',data.full) ' ];']);
fprintf(fid,['\nW_Q   = [' sprintf(' %8.4f',data.W_res(1,:)) ' ];']);
fprintf(fid,['\nW_V   = [' sprintf(' %8.4f',data.W_res(2,:)) ' ];']);
fprintf(fid,['\nW_D   = [' sprintf(' %8.4f',data.W_res(3,:)) ' ];']);
if ~isempty(data.aIonic)
   fprintf(fid,['\n\nW_I   = ' sprintf(' %8.4f',data.W_I) ';']);
end

fprintf(fid,'\n\n%%        ');
fprintf(fid,' %8s',data.atoms{:});

fprintf(fid,['\nW_C   = [' sprintf(' %8.4f',data.W_C) ' ];']);
if isempty(data.Caggr) 
   fprintf(fid,['\nCaggr = [' sprintf(' %8d',1:data.n) ' ];']);
else
   fprintf(fid,['\nCaggr = [' sprintf(' %8d',data.Caggr) ' ];']);
end
fprintf(fid,['\nQ     = [' sprintf(' %8d',data.Q) ' ];']);
fprintf(fid,['\nZ     = [' sprintf(' %8d',data.Z) ' ];']);

if ~isempty(data.aIonic)
   fprintf(fid,'\naIonic  =%8.4f;',data.aIonic);
end

fprintf(fid,'\n        ');
for k=1:data.n, fprintf(fid,' %8s',' '); end
fprintf(fid,'  %% anion\nR     = {');
for l=1:size(data.R,1)
    fprintf(fid,' %8.4f',[data.R{l,1:data.n}]);
    fprintf(fid,' %5s',sprintf('''%s''',data.R{l,end}));
    if l<size(data.R,1)
        fprintf(fid,'\n         ');
    else
        fprintf(fid,' };\n\n');
    end
end

fprintf(fid,'\n%% CALCULATED VALUES\n\n');

fprintf(fid,'%% vacancy\n%%');
fprintf(fid,' %8s',data.sites{:});
fprintf(fid,['\n[' sprintf(' %8.4f',data.vac) ' ];\n\n']);

fprintf(fid,'%% calculated chemistry\n%%');
fprintf(fid,' %8s',data.atoms{:});
fprintf(fid,['\n[' sprintf(' %8.4f',data.chem_calc) ' ];\n\n']);

fprintf(fid,'%% chemistry vacancy:\n');
fprintf(fid,'%%   calc''d  obsrv''d\n[');
fprintf(fid,' %8.4f',data.cc_vac,sum(data.C)-sum(data.Cobs));
fprintf(fid,' ];\n\n');

fprintf(fid,'%% residuals\n%%');
fprintf(fid,' %8s',data.sites{:});
fprintf(fid,['\n[' sprintf(' %8.4f',data.res(1,:))...
             '     %% scattering']);
fprintf(fid,['\n ' sprintf(' %8.4f',data.res(2,:)) '     %% valence']);
fprintf(fid,['\n ' sprintf(' %8.4f',data.res(3,:)) ]);
fprintf(fid,' ];  %% length\n\n');

if ~isempty(data.aIonic)
   fprintf(fid,'%% ionic total for cations\n');
   fprintf(fid,' %8.4f;\n',data.cIonic);
end

fprintf(fid,'%% occupancy as a fraction of calc''d chemistry\n%%');
fprintf(fid,' %8s',data.sites{:});
fprintf(fid,'\n[');
for k=1:data.n
    if k~=1, fprintf(fid,' '); end
    fprintf(fid,' %8.4f',data.Cx(k,:));
    if k==data.n
        fprintf(fid,' ];');
    else
        fprintf(fid,'   ');
    end
    fprintf(fid,' %% %s\n',data.atoms{k});
end
fprintf(fid,['[' sprintf(' %8.4f',data.Cvac(:)) ' ]; %% vacancy\n\n']);


fprintf(fid,'%% bond-valence (as calculated from d and R)\n%%');
fprintf(fid,' %8s',data.sites{:});
fprintf(fid,'\n[');
for k=1:data.n
    if k~=1, fprintf(fid,' '); end
    fprintf(fid,' %8.4f',data.V(:,k));
    if k==data.n
        fprintf(fid,' ];');
    else
        fprintf(fid,'   ');
    end
    fprintf(fid,' %% %s\n',data.atoms{k});
end

fprintf(fid,'\n%% distance (as calculated from Z and R)\n%%');
fprintf(fid,' %8s',data.sites{:});
fprintf(fid,'\n[');
for k=1:data.n
    if k~=1, fprintf(fid,' '); end
    fprintf(fid,' %8.4f',data.D(:,k));
    if k==data.n
        fprintf(fid,' ];');
    else
        fprintf(fid,'   ');
    end
    fprintf(fid,' %% %s\n',data.atoms{k});
end
end % OF FUNCTION write_Mfile


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2. Optimization code: convert data to model and apply numerical solver.         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [X,errstr]=...
   occup(Q,Qobs,Z,V,C,Cobs,aIonic,D,dbar,W,fixed,xval,Caggr)
% Implementation, by the first author, of the core mathematical calculation
% described in the article "Optimization of site-occupancies in minerals
% using quadratic programming" by S.E. Wright, J.A. Foley, and J.M. Hughes,
% American Mineralogist 85 (2000), 524-531.
%
% This function sets up the QP problem and calls a solver. The larger
% OccQP program posted in p-file format at
%    http://www.users.miamioh.edu/wrightse/occqp/
% provides an outdated interface to assist with managing the inputs and
% outputs of the current function.

% Note that the OCCUP function relies on subfunctions later in this file.

[m,n]=size(V);

% Convert aggregation of atomic species
CC = eye(n);
if ~isempty(Caggr)
   [aggrs,agkey] = sort(Caggr);
   k0 = 1;
   while k0<=n
      k1 = find(aggrs==aggrs(k0),1,'last');
      CC(agkey(k0:k1),agkey(k0:k1)) = ones(k1-k0+1);
      k0 = k1+1;
   end
end

% Set up coefficient matrix and target value for goal equations
N = m*n;
M = 4*m+n+1;
Coeff = zeros(M,N);
rowres = 1+m*(0:3); % This is slid forward during the loop below.
rowC = 4*m+(1:n);
rowion = M;
col = 1:n;
for site=1:m
    Coeff([rowres,rowC,rowion],col) = [
        ones(1,n)           % occupancy
        Q                   % scattering
        V(site,:)-Z         % valence
        D(site,:)-dbar(site)% distance
        C(site)*CC          % chemistry
        Z*C(site)           % ionic balance
    ];
    rowres = rowres+1;
    col = col+n;
end
Target = [
    ones(m,1)  % occupancy
    Qobs'      % scattering
    zeros(m,1) % valence
    zeros(m,1) % distance
    CC*Cobs'   % chemistry
    aIonic     % ionic balance
];

% Identify constraints
infW = isinf(W);
finW = ~infW;
eqcon = find(infW);
ineqcon = find(finW(1:m));
obj = m+find(finW(m+1:end));

% Create data for a linear-equality constrained NNLS problem
A = diag(sqrt(W(obj))) * Coeff(obj,:);
b = sqrt(W(obj)) .* Target(obj);
E = Coeff(eqcon,:);
f = Target(eqcon);
G = Coeff(ineqcon,:);
h = Target(ineqcon);
lo = zeros(n*m,1);
% up = ones(n*m,1);
up = inf*ones(n*m,1);
lo(fixed) = xval(fixed);
up(fixed) = xval(fixed);
% x0 = xval(:); % can't pass x0 to LHnnls, so don't define it

[x,how] = use_LHnnls(A,b,E,f,G,h,lo,up);

% How did we do?

X = []; % Assume the worst.
if how==1
    errstr = sprintf(...
        ['\n  Solver exceeded its maximum iteration count.\n'...
        '  Please report this problem to the developers.\n'...
        '  Include a copy of an input data file causing'...
        ' the error.\n']);
    return
elseif any(size(x)~=[m*n,1])
    errstr = sprintf(...
        ['\n  Solver returned an invalid answer.\n'...
        '  Please report this problem to the developers.\n'...
        '  Include a copy of an input data file causing'...
        ' the error.\n']);
    return
elseif any(norm(E*x-f,inf)>8*eps(norm([E,f],inf))) ...
        || any(G*x-h>8*eps(norm([G,h],inf)))
   errstr = sprintf(...
    ['\n  Problem is overly restrictive! Try using fewer'...
          ' infinite weights,\n'...
       '  fewer fixed occupancies, or fewer full sites.']);
   return
else % Otherwise, we actually found something
    errstr = [];
    x=abs(x); % cheating to eliminate unsightly negative signs
              % Note: may also prevent detecting solver errors!
    X=reshape(x,[n,m]);
end
end % OF FUNCTION occup


function [x,how] = use_LHnnls(A,b,E,f,G,h,lo,up)
% Apply Lawson & Hanson's NNLS algorithm, but with weighted equations

% Build data for weighted NNLS problem
mA = size(A,1);
mE = size(E,1);
mG = size(G,1); % # slack variables
Ifix = lo==up; % fixed variables
J = ~Ifix;
% rescale factors for weighting
mu = norm(A,1) / sqrt(eps)/ (norm(E,1) + max(1,norm(G,1)));
% stack data for nnls (remove fixed variables, add slack variables)
C = [
    E(:,J)*mu  , zeros(mE,mG)
    G(:,J)*mu  , diag(mu(ones(mG,1))) % mu*[G,eye(mG)]
    A(:,J)     , zeros(mA,mG)
    ];
if any(Ifix)
    d = [
        (f-E(:,Ifix)*lo(Ifix)) *mu
        (h-G(:,Ifix)*lo(Ifix)) *mu
        b-A(:,Ifix)*lo(Ifix)
        ];
else
    d = [
        f*mu
        h*mu
        b
        ];
end

% Solve the NNLS problem
[w,~,lambda] = LHnnls(C,d);

% Report the outcome
x = lo;
x(J) = w(1:sum(J));
if isempty(lambda)
    how = 1;
else
    how = 0;
end
end % OF FUNCTION use_LHnnls


function [x,rnorm,w] = LHnnls(A,b,x0)
% [x,rnorm,w] = LHnnls(A,b)
%
% Solves the nonnegative least squares problem
%
%    min norm(A*x-b,2) over all x>=0 .
%
% The returned arguments include x, the residual magnitude rnorm==norm(A*x-b,2),
% and the KKT multiplier vector w.  The routine was successful if and only if
% the returned value of w is nonempty; if w is empty, then the maximum iteration
% count 3*size(A,2) was exceeded.
%
% An optional 3rd input argument x0 can be used to start the algorithm.

% The code in this function is a fairly "faithful" translation of Lawson & Hanson's
% FORTRAN subroutine NNLS.  The main differences are as follows:
% 1. FORTRAN "goto" structures have been replaced by logically equivalent "while true"
%    structures or inlined code.
% 2. FORTRAN "do" structures used for common array operations (i.e., vector copying,
%    array index listing, matrix multiplication, sums, norms, finds, maximums, minimums)
%    have been replaced by standard MATLAB operations and functions.
% 3. FORTRAN "subroutine" (if not inlined) are replaced by MATLAB functions, which
%    are necessarily pass-by-value rather than pass-by-reference.
% A general code based on the exposition in Lawson & Hanson's book is given in
% the function 'swnnls.m'.

[m,n] = size(A);
factor = 0.01;

iter = 0;
itmax = 3*n;

% INITIALIZE x AND Z (= complement of P)
x = zeros(n,1);
if nargin==3
    x(x0>0) = x0(x0>0);
end
nsetp = sum(x>0);
npp1 = nsetp+1;
iz1 = npp1;
iz2 = n;
index = zeros(n,1); % pre-allocate
index(1:nsetp) = find(x>0);
index(npp1:end) = find(x==0);

% REDUCE MATRIX IF STARTING GUESS PROVIDED
if nsetp>1
    if m<n+1 % introduces m-by-m Q
        [Q,A(:,index(1:nsetp))] = qr(A(:,index(1:nsetp)));
        A(:,index(iz1:iz2)) = Q'*A(:,index(iz1:iz2));
        b = Q'*b;
    else % introduces m-by-(n+1) Ab
        Ab = qr( [A(:,index([1:nsetp,iz1:iz2])) , b] );
        A(:,index([1:nsetp,iz1:iz2])) = triu(Ab(:,1:end-1));
        b = [Ab(1:min(m,n),end); zeros(m-min(m,n),1)];
    end
end

w = zeros(n,1); % pre-allocate
zz = zeros(m,1);

while true % MAIN LOOP BEGINS HERE
    % CHECK FOR SUCCESS OR OPPORTUNITY FOR IMPROVEMENT
    w(index(iz1:iz2)) = A(npp1:m,index(iz1:iz2))'*b(npp1:m);
    while true % exits with success or candidate for improvement
        [wmax,i] = max(w(index(iz1:iz2)));
        izmax = iz1-1+i;
        if isempty(wmax) || wmax<=0 % SUCCESSFUL TERMINATION
            rnorm = norm(b,2);
            return
        end
        iz = izmax;
        j = index(iz); % j is a possible candidate
        Asave = A(npp1,j);
        [A(npp1,j),upivot] = h1(npp1,npp1+1,A(:,j));
        unorm = norm(A(1:nsetp,j),2);
        if (unorm+abs(A(npp1,j))*factor)-unorm>0
            zz = b;
            zz = h2(npp1,npp1+1,A(:,j),upivot,zz);
            if zz(npp1)/A(npp1,j)>0
                break % yes, it is a candidate
            end
        end
        % otherwise, try again
        A(npp1,j) = Asave;
        w(j) = 0;
    end
    
    % MOVE INDEX j INTO SET P AND UPDATE FACTORIZATION
    index(iz) = index(iz1);
    index(iz1) = j;
    iz1 = iz1+1;
    nsetp = npp1;
    npp1 = npp1+1;
    A(:,index(iz1:iz2)) = h2(nsetp,npp1,A(:,j),upivot,A(:,index(iz1:iz2)));
    A(npp1:m,j) = 0;
    w(j) = 0;
    
    % SOLVE LLS OVER SET P
    b = zz;
    zz(1:nsetp) = A(1:nsetp,index(1:nsetp))\b(1:nsetp);
    
    % SECONDARY LOOP: TAKE LARGEST FEASIBLE STEP AND UPDATE P
    while true 
        iter = iter+1;
        if iter>itmax % ITERATION COUNT EXCEEDED
            rnorm = norm(b,2);
            w = [];
            return
        end
        
        % DETERMINE MAXIMUM FEASIBLE STEP
        alpha = 2;
        for ip=1:nsetp
            if zz(ip)<=0
                t = -x(index(ip))/(zz(ip)-x(index(ip)));
                if alpha>t
                    alpha = t;
                    jj = ip;
                end
            end
        end
        
        if alpha==2 % full step is allowed
            break % out of secondary loop
        end
        % TAKE MAXIMUM FEASIBLE STEP
        x(index(1:nsetp)) = x(index(1:nsetp)) + ...
            alpha*(zz(1:nsetp)-x(index(1:nsetp)));
        
        % REMOVE NEW NONPOSITIVES FROM SET P
        i = index(jj);
        while true % stops when all members of P are positive
            x(i) = 0;
            if jj~=nsetp
                jj = jj+1;
                for j=jj:nsetp
                    ii = index(j);
                    index(j-1) = ii;
                    % build rotation
                    aa = A(j-1,ii);
                    bb = A(j,ii);
                    if abs(aa)>abs(bb)
                        xr = bb/aa;
                        yr = sqrt(1+xr^2);
                        cc = sign(aa)/yr;
                        ss = cc*xr;
                        A(j-1,ii) = abs(aa)*yr;
                    elseif bb~=0
                        xr = aa/bb;
                        yr = sqrt(1+xr^2);
                        ss = sign(bb)/yr;
                        cc = ss*xr;
                        A(j-1,ii) = abs(bb)*yr;
                    else
                        cc = 0;
                        ss = 1;
                        A(j-1,ii) = zero;
                    end
                    cs = [cc,ss;-ss,cc];
                    % apply rotation
                    A(j,ii) = 0;
                    L=[1:ii-1,ii+1:n];
                    A([j-1,j],L) = cs*A(j-1:j,L);
                    b(j-1:j) = cs*b(j-1:j);
                end
            end
            npp1 = nsetp;
            nsetp = nsetp-1;
            iz1 = iz1-1;
            index(iz1) = i;
            jj = find(x(index(1:nsetp))<=0,1,'first');
            if isempty(jj), break, end
            i = index(jj);
        end
        
        % SOLVE LLS OVER SET P
        zz(1:nsetp) = A(1:nsetp,index(1:nsetp))\b(1:nsetp);
        
    end % SECONDARY LOOP
    
    x(index(1:nsetp)) = zz(1:nsetp);
    
end % MAIN LOOP
end % OF FUNCTION LHnnls

%%%% helper functions for 'LHnnls' are below %%%%

function [unew,upivot] = h1(ipivot,i1,u)
% Given u, returns scalar pair for u(ipivot) defining Householder
% transformation to zero out u(i1:end) and place mass in u(ipivot).
m = length(u);
unew = u(ipivot);
if ipivot<=0 || ipivot>=i1 || i1>m
    upivot = 0;
    return
end
c = norm( u([ipivot,i1:m]) ,2);
if c==0, return, end
if unew>0
    c = -c;
end
upivot = unew-c;
unew = c;
end % OF FUNCTION h1

function V=h2(ipivot,i1,u,upivot,V)
% Given V and arguments from prior call to h1, returns result of
% applying specified Householder transformation to V.
m = length(u);
if ipivot<=0 || ipivot>=i1 || i1>m, return, end
bb = upivot*u(ipivot);
if bb>=0, return, end
S = ( [upivot;u(i1:m)]'*V([ipivot,i1:m],:) )/bb;
V([ipivot,i1:m],:) = V([ipivot,i1:m],:) + [upivot;u(i1:m)] * S;
end % OF FUNCTION h2
