function [L] = liusimc(L,LX,configdd)
%LIUSIMC  Consolidation for running lihtslx.slx or liutslx.slx 
% - Simulink version of liuqe / lih.
% [L,LX] = liusimc(L,LX,configdd)
%
% configdd is an optional data dictionary name that contains the
% Simulink configuration. If empty, a default configuration is added in the
% data dictionary itself. This is used when operated as referenced
% subsystem with externaly supplied common configuration.
%
% [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved.

% don't use sldd if in GITLAB_CI environment due to concurrency issues
use_sldd = (isempty(getenv('GITLAB_CI')) ...
  || ~strcmp(getenv('CI_PROJECT_NAME'),'meq'));
use_sldd = use_sldd && ~(isfield(L.P,'nosldd') && L.P.nosldd); % optional override

if nargin==2
  configdd={};
end

assert(logical(L.P.slx),'L.P.slx must be 1 to call liusimc')
assert(~verLessThan('matlab','9.6.0'),'MEQ Simulink models now require MATLAB version>=9.6.0')

liuhsmodel = [L.code,'sim'];
liuhtmodel = [L.code,'tslx'];
liuhxmodel = 'liuxslx'; % common data source

%% Set variant subsystems: parameters for switches
switcher_param = setVariants(L);
Pswitchers = [fieldnames(switcher_param),struct2cell(switcher_param)]'; % turn into param-argument pairs

%% create buses
% LX bus
if ~isfield(LX,'Ip') || isempty(LX.Ip), LX.Ip = zeros(size(LX.t)); end
LX.Ie = [LX.Ia;LX.Iu];
LXfields = {'Ff','Bm','Ip','Ie','rBt','Ft'}; % LX fields for bus
for ifield = LXfields
  LXstruct.(ifield{:}) = single(LX.(ifield{:})(:,1));
end
businfo = Simulink.Bus.createObject(LXstruct);
LXbus = evalin('base',businfo.busName);
LXbus = vectorBusDimensions(LXbus);
evalin('base',sprintf('clear %s',businfo.busName)); % remove temp bus in base

% LY bus
LYstruct = getTemplateLYstruct(L,LX);
businfo = Simulink.Bus.createObject(LYstruct);
LYbus = evalin('base',businfo.busName);
LYbus = vectorBusDimensions(LYbus,{'FW'});
evalin('base',sprintf('clear %s',businfo.busName)); % remove temp bus in base

%% Data conversion
% cast L entries to singles to avoid precision loss warnings
L = struct2single(L);

% Data dictionary cannot store function handle. Cast fh to strings.
if isa(L.P.bfct,'function_handle'), L.P.bfct = func2str(L.P.bfct); end
if isa(  L.bfct,'function_handle'),   L.bfct = func2str(  L.bfct); end
if isa(L.P.zfct,'function_handle'), L.P.zfct = func2str(L.P.zfct); end
if ~isempty(L.P.infct) && isa(L.P.infct,'function_handle'), L.P.infct = func2str(L.P.infct); end

if strcmp(L.P.bfct,'bfgenD'), for ib = 1:size(L.P.bfp,1), L.P.bfp{ib,1} = func2str(L.P.bfp{ib,1});end;end
if strcmp(  L.bfct,'bfgenD'), for ib = 1:size(  L.bfp,1),   L.bfp{ib,1} = func2str(  L.bfp{ib,1});end;end

% Set pseudo time with constant time steps
if numel(LX.t)>1
  ok = numel(LX.t)<=2 || max(abs(diff(diff(LX.t))))<1e-6;
  if ~ok && L.nn, warning('LX.t not equispaced, loop voltage may yield wrong results'); end
  L.dt = round(mean(diff(LX.t)),4);
else
  L.dt = 1; % dummy
end

%% Configuration settings
configSet = Simulink.ConfigSet; % generate default configSet
setConfigTimes(configSet,L,LX); % set the appropriate values

%% Data dictionaries
if use_sldd
  % Generate data dictionaries
  
  % Dictionary names
  liutdictionaryName   = [liuhtmodel,'.sldd'];
  liuxdictionaryName   = [liuhxmodel,'.sldd'];
  liusdictionaryName   = [liuhsmodel,'.sldd'];
  
  % Close data dictionaries
  listopendictionaries = Simulink.data.dictionary.getOpenDictionaryPaths;
  for mydict = {liutdictionaryName,liuxdictionaryName,liusdictionaryName}
    if numel(listopendictionaries)>0 && any(contains(listopendictionaries, mydict{:}))
      Simulink.data.dictionary.closeAll(mydict{:}, '-discard');
    end
  end
  
  % Put data into model data dictionary
  % generate data dictionaries
  liutddobj = generateDD(liutdictionaryName);
  liuxddobj = generateDD(liuxdictionaryName);
  liusddobj = generateDD(liusdictionaryName);
  % add variables
  replaceOrCreateDDentry(getSection(liutddobj,'Design Data'),'L',L,'LXbus',LXbus,'LYbus',LYbus,Pswitchers{:});
  replaceOrCreateDDentry(getSection(liuxddobj,'Design Data'),'LX',LX);
  
  if isempty(configdd)
    % add configurationSettings in this data dictionary directly
    replaceOrCreateDDentry(getSection(liutddobj,'Configurations'),'configurationSettings',configSet)
  else
    configSection = getSection(liutddobj,'Configurations');
    if configSection.exist('configurationSettings') && ... % entry is present
        isequal(configSection.getEntry('configurationSettings').DataSource,liutddobj) % entry is local to the data dictionary
        % delete it locally
        deleteDDentry(configSection,'configurationSettings'); % delete local configuationSettings to avoid conflicts
    end
    setReferencedDD(liutddobj,configdd);
  end
  
  % add data sources
  setReferencedDD(liuxddobj,liutdictionaryName);
  setReferencedDD(liusddobj,liutdictionaryName,liuxdictionaryName)
  
  % Save & close
  try
  if liutddobj.HasUnsavedChanges, liutddobj.saveChanges; end
  if liuxddobj.HasUnsavedChanges, liuxddobj.saveChanges; end 
  if liusddobj.HasUnsavedChanges, liusddobj.saveChanges; end
  catch ME
    % if write permissions were denied, write .sldd to local folder
    if strcmp(ME.identifier,'MATLAB:save:permissionDenied')
      error('no permission to save data dictionaries, set ''nossldd''=true or use a meq instance where you have write permissions');
    else
      rethrow(ME);
    end
  end
else
  % Assign variables in base workspace
  assignin('base','L',L);
  assignin('base','LX',LX);
  assignin('base','LXbus',LXbus);
  assignin('base','LYbus',LYbus);
  for ii = 1:2:numel(Pswitchers), assignin('base',Pswitchers{ii:ii+1}); end

  % assign configurationSettings to base workspace
  assignin('base','configurationSettings',configSet);
 
  % remove sldd link in model & enable base workspace access

  w = warning('off');  % turn off warning for missing variables due to missing sldd
  args = {'EnableAccessToBaseWorkspace','on'};
  for file = {liuhtmodel,liuhxmodel,liuhsmodel}
    hh = load_system(file{:});
    set_param(hh,args{:},'DataDictionary','');
  end 
  warning(w); % reset warning state
end

end


function LYstruct = getTemplateLYstruct(L,LX)
    
switch L.code
  case 'liu'
    L.liurtemu = true;
    LY = liut(L,meqxk(LX,1),'slx',false); % one step of liut - must use iterq to assign fields
    assert(~isempty(LY),'liut default call failed')
    L.liurtemu = false;
    LY.Ie = [LY.Ia;LY.Iu];
    fieldsrm = {'Ia','Iu','Iv','Is',...
      'shot','nB','PpQg','TTpQg','rIp','zIp',...
      'qA','btD','VpD','bpD','Ft0D','FtD','VpD','liD',...
      'WND','Wt0D','WtD','WpD','WkD','FtPQ','VpQ','IpQ',...
      'rrmax','rrmin','zzmax','zzmin','rzmax','rzmin','zrmax','zrmin'}; 
    if L.P.nFW==0  , fieldsrm = [fieldsrm,{'FW'}]; end
    if L.P.iterq==0, fieldsrm = [fieldsrm,{'aq','aW'}]; end
    for ii=fieldsrm
      if isfield(LY,ii{:}), LY = rmfield(LY,ii{:}); end
    end
  case 'lih'
    % LY bus
    LY = liht(L,meqxk(LX,'last')); % one step of liht
    LY.Ie = [LY.Ia;LY.Iu];
    keepfields = {'t','Ip','zIp','rIp','IpD','zIpD','rIpD','Ia','Iu','Iy',...
      'Fx','lB','rA','zA','nA','Ie'};
    LY = rmfield(LY,setdiff(fieldnames(LY),keepfields));
end

LY.err = int32(0);
for ifield=fieldnames(LY)'
  if isempty(LY.(ifield{:}))
    LYstruct.(ifield{:}) = single(0); % placeholder
  elseif isa(LY.(ifield{:}),'double')
    LYstruct.(ifield{:}) = single(LY.(ifield{:}));
  else
    LYstruct.(ifield{:}) = LY.(ifield{:});
  end
end
% special sizes for rX,zX,FX
[LYstruct.rX,LYstruct.zX,LYstruct.FX] = deal(single(zeros(L.dimw{:},1))); % memory size for rX,zX,FX

end

function dictionaryObj = generateDD(dictionaryName)
% Add data in varargin to the dictionary
dictdir = fileparts(mfilename('fullpath'));
dictionaryPath = fullfile(dictdir,dictionaryName);
if exist(dictionaryPath, 'file')
  dictionaryObj = Simulink.data.dictionary.open(dictionaryPath);   % Open the existing data dictionary
else
  try 
    dictionaryObj = Simulink.data.dictionary.create(dictionaryPath); % create new
  catch ME
    % if write permissions were denied, write .sldd to local folder
    if strcmp(ME.identifier,'MATLAB:save:permissionDenied')
      error('no permission to write `.sldd` to %s, set ''nossldd'',true',dictionaryPath);
    else
      rethrow(ME);
    end
  end
end

end

function [out_struct] = setVariants(L)
% set variant subsystem based on parameters
  
% set variant subsystem based on parameters (for lihsim)

switch L.code
  case 'lih', variantList = {'HFCTMODE','IPMODE','DBLTMODE'};
  case 'liu', variantList = {'BFCTMODE','HFCTMODE','IPMODE','BSLVMODE','INFCTMODE','AWMODE',...
                             'LOCSMODE','LOCRMODE','ITERQMODE','KYWMODE'};
  otherwise, error('code %s not supported in liusimc',L.code);
end

switcher_values = struct();
SimVariant = struct();

for variant = variantList
  switch variant{1}
    case 'BFCTMODE'
      switch func2str(L.bfct)
        case 'bf3pmex'
          value=1;
        case 'bf3imex'
          value=2;
        case 'bfefmex'
          value=3;
        otherwise
          error('invalid bfct %s for liutsim',func2str(L.bfct));
      end
      valueList = [1,2,3];
    case 'HFCTMODE'
      if L.P.iterh==0
        value=1; % standard FE solver
      elseif isfield(L.P,'ipmhdprec') && ~L.P.ipmhdprec
        value=2; % constrained FE solver single precision
      else
        value=3; % constrained FE solver double precision
      end
      valueList = [1,2,3];
    case 'IPMODE'
      if L.P.Ipmeas
        value=1; % with Ip meas
      else
        value=2; % no Ip meas - trapeze estimator
      end
      valueList = [1,2];
    case 'INFCTMODE'
      if ~L.nn
        value = 1;
      elseif isequal(func2str(L.P.infct),'qintmex')
        value = 2;
      else
        error('unsupported infct for Simulink')
      end
      valueList = [1,2];
    case 'BSLVMODE'
      if isfield(L.P,'bslvdprec') && ~L.P.bslvdprec
        value=1; % single precision
      else
        value=2; % double precision for bslv & iata
      end
      valueList = [1,2];
    case 'AWMODE'
      if ~L.nW
        value = 0;
      else
        value = 1;
      end
      valueList = [0,1];
    case 'LOCSMODE'
      if isempty(L.P.raS)
        value = 0;
      else
        value = 1;
      end
      valueList = [0,1];
    case 'LOCRMODE'
      if isempty(L.P.iqR)
        value = 0;
      else
        value = 1;
      end
      valueList = [0,1];
    case 'ITERQMODE'
      if L.P.iterq == 0
        value = 0;
      else
        value = 1;
      end
      valueList = [0,1];
    case 'KYWMODE'
      if isempty(L.kyw)
        value = 0;
      else
        value = 1;
      end
      valueList = [0,1];
    case 'DBLTMODE'
      if L.P.idoublet
        value=2;
      else
        value=1;
      end
      valueList = [1,2];
    otherwise
      error('Unrecognized variant name in liusimc: %s',variant{1});
  end
  switcher_values.(variant{1}) = value;
  for v = valueList
    SimVariant.(sprintf('liu%svariant%d',variant{1},v)) = Simulink.Variant(sprintf('%s==%d',variant{1},v));
  end
end

% return outputs
if nargout > 0
    out_struct = switcher_values;
    for field=fieldnames(SimVariant).'
         out_struct.(field{1}) =  SimVariant.(field{1});
    end
end
end

function replaceOrCreateDDentry(ddSectionObj,varargin)
for k=1:2:numel(varargin)
  entry = varargin{k};
  value = varargin{k+1};
  if ddSectionObj.exist(entry)
    oldEntry = ddSectionObj.getEntry(entry);
    assert(numel(oldEntry)==1,'multiple entries found for %s',entry)
    if ~isequal(oldEntry.getValue,value)
      oldEntry.setValue(value); % replace
    end
  else
    ddSectionObj.addEntry(entry,value);
  end
end
end

function deleteDDentry(ddSectionObj,varargin)
for k=1:numel(varargin)
  entry = varargin{k};
  if ddSectionObj.exist(entry)
    ddSectionObj.deleteEntry(entry)
  end
end
end

function setReferencedDD(dictionaryObj,varargin)
for k = 1:numel(varargin)
  objval = varargin{k};
  % First remove previous DataSources which aren't needed anymore
  for iSource = reshape(setdiff(dictionaryObj.DataSources,objval),1,[])
    dictionaryObj.removeDataSource(iSource{:});
  end
  % Now add new sources
  for iSource = reshape(setdiff(objval,dictionaryObj.DataSources),1,[])
    dictionaryObj.addDataSource(iSource{:})
  end
end
end

function setConfigTimes(configSet,L,LX)
    set_param(configSet,...
      'name','configurationSettings',...
      'StartTime','0',...
      'StopTime' ,num2str((numel(LX.t)-1)*L.dt,'%.4f'),...
      'SolverType', 'Fixed-step',...
      'FixedStep','L.dt'...
      );
end

function myBus = vectorBusDimensions(myBus,exceptions)
if nargin==1, exceptions = {}; end

  % set matrix-sized dimensions [n 1] to array dimensions [n]
  for iElem = 1:numel(myBus.Elements)
    Elem = myBus.Elements(iElem);
    if any(contains(Elem.Name,exceptions)), continue; end
    if isequal(size(Elem.Dimensions),[1 2]) && Elem.Dimensions(2) == 1
      myBus.Elements(iElem).Dimensions = Elem.Dimensions(1);
    end
  end
end
