function agconc = meqagconc(L)
% parser for ag constraint specifications.
%
% Reads constraint specifications from L.P.agcon, parses it, and stores the
% result in L.agconc
%
% Specify constraint either as strings (e.g. 'Ip') which correspond to
% functions in meqagconfun.m
%
% Several ways of input:
%
% 1) Simple constraints, e.g. {'Ip','bp','li'}
%    This directly specifies the basis function constraints to be used.
%    For doublets, this defaults to having these constraints for the two
%    lobes, and fixed coefficients for mantle.
%
% 2) Constraints per domain, e.g. {{'Ip','bp','li'},{'Ip','bt','qA'},'{'ag','ag','ag'}
%    This explicitly specifies the constraints for each domain.
%
% 3) Direct specification in final format e.g.
%    {'Ip',1;'Ip',2,'bp',1,'bp',2,'Ip',3,'bp',1};
%    This specifies the constraint to be used and the corresponding domain.
%    Empty domain means the constraint applies to the whole plasma.
%
% outputs:
%  agconc is a cell array with 3 columns:
%     1) function handle for a given constraint
%     2) domains to which the constraint applies
%     3) corresponding quantity in LX entering into that constraint.
%  This output is used by meqagcon.m to build the basis function constraints
%  
%
% See meqagcon.m, meqagconfun.m
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

ll = meqagconfun(); % get catalog of functions
in = L.P.agcon;
assert(iscell(in),'L.P.agcon must be a cell');

%#ok<*AGROW>
out = cell(0,2); % init
if size(in,2) ==2 && (isnumeric(in{1,2}) || isempty(in{1,2})) 
  % Input mode 3: direct input of cell of constraints
  % basic checks/conversions
  for ii=size(in,1)
    switch class(in{ii,1})
      case 'function_handle' % ok keep it
      case 'char',  in(ii,1) = {myfun(ll,in{ii})}; % convert to func
      otherwise, error('don''t know hwo to handle this type')
    end
    nval = in{ii,2};
    assert(isempty(nval) || isnumeric(nval),...
      'AGCONPARSE:DomainIndicatorType','domain indicator in second column must be empty or numeric')
    assert(isempty(nval) || (nval>0 && nval<=L.nD),...
      'AGCONPARSE:DomainIndicatorValue','domain indicator must be empty or between 0 and nD=%d',L.nD);
  end
  out = in;
elseif size(in,1)==1 % input a row of constraints
  if ~iscell(in{1}) % mode(1): single constraint set given: e.g. {'Ip','bp','qA'}
    if ~L.P.idoublet
      % singlet: only one set of constraints for one domain
      out = filliD(out,in,ll,1);
    else
      % doublet: assume same for two lobes, then fix ag for mantle
      ng = sum(L.TDg(3,:)); % Number of basis functions for mantle.
      in = {in,in,repmat({'ag'},1,ng)};
      for iD=1:numel(in)
        out = filliD(out,in{iD},ll,iD);
      end
    end
  else % mode (2) multiple constraints are given
    for iD=1:numel(in)
      out = filliD(out,in{iD},ll,iD);
    end
  end
end

% find corresponding LX field and its index that enters constraint, if any
for icon=1:size(out,1)
  [fieldname] = meqagconq(out,icon,L.P.idoublet);
  iD = out{icon,2}; % domain index of this constraint
  switch fieldname
    case 'ag'
      % special case: keep ag value for this constraint equal to the
      % previous one. Need to find index of correct ag to use for this
      % constraint in meqagcon.
      
      igD = find(L.TDg(iD,:));                % Indices of ag active in this domain
      ii = igD(find([out{:,2}]==iD) == icon); % index within constraints with this domain index
      
    otherwise
      ii = iD; % for other fields, take iDth index (for multidomain)
  end
  
  LXfieldname{icon,1} = fieldname;
  LXindex    {icon,1} = ii;
end

agconc = [out,LXfieldname,LXindex];

if size(agconc,1) ~= L.ng
  % check number of constraints
  disp(agconc)
  error('found %d constraints, expected %d\n',size(agconc,1),L.ng);
end
end

function out = filliD(out,in,ll,iD)
% fill some rows of table for the constraints of this domain
for ii = 1:numel(in)
  out = [out;{myfun(ll,in{ii}),iD}];
end
end

function fun = myfun(ll,myfunstr)
% find function handle from list
llstr = cellfun(@func2str,ll,'UniformOutput',false); % convert to string
ill = strcmp(llstr,myfunstr);
assert(any(ill),   'contraint function %s not found in meqagconfun',myfunstr)
assert(sum(ill)==1,'multiple contraint functions %s found in meqagconfun',myfunstr)
fun = ll{ill};
end

function [agconq] = meqagconq(agconc,icon,idoublet)
% return string of field in LX, that enters in constraint function agcon
% with agcon in format {function_handle,domain_index}

qstr = func2str(agconc{icon,1});
iD   = agconc{icon,2};

% fields that can be computed for each separate domain, that have 'D'
% counterpart
fieldsD = {'Ip','Wk','Wp','Wt','Wt0','WN','Vp','Ft','Ft0','bp','bp','Ini','li'};

if ~idoublet || isempty(iD)
  agconq = qstr; % fields like 'Ip'
else
  if contains(qstr,fieldsD)
    agconq = [qstr,'D'];   % fields like 'IpD'
  else
    agconq = qstr;
  end
end

end
