%%%%%%%%  ARTIFICIAL NEURAL NETWORK BASED ON QUALIKIZ REGRESSION
%%%%%%%%  KAREL VAN DE PLASSCHE JONATHAN CITRIN

%%%%%%%% IN PRACTICE THE REGULARIZED NETWORKS CAN OFTEN REASONABLY EXTRAPOLATE, BUT CAVEAT EMPTOR!

%%%%%%%% TRANSPORT MODEL DEVELOPED FOR USE IN RAPTOR (F. FELICI et al)

%%%%%%%% Karel 2015a version. Needs to be converted to output
%%%%%%%% RAPTOR-compatible outputs!

function varargout = QLKNN10D(Nin,net,pp)
% #codegen

coder.extrinsic('sprintf_general');


if nargout>1 % calculate derivatives?
  calcder = true;
else
  calcder = false;
end

if nargin==0
end

if nargin==0
  
  % call to set up QLKNN
  RAPTORpath = RAPTOR_path();
  submodulePath = fullfile(RAPTORpath(),'submodules');
  hyperNamelistPath = fullfile(submodulePath,'qlknn-hyper-namelists');
  qlknn_network.namelist_dir = hyperNamelistPath;
  qlknn_network.name = 'QLKNN10D';
  qlknnPath = fullfile(submodulePath,'QLKNN-fortran');

  % Initizalize QLKNN and namelists
  cmd = ['git -C ', submodulePath, ' submodule init QLKNN-fortran qlknn-hyper-namelists'];
  cmd_res = system(cmd);
  if cmd_res ~= 0
      error('Problem initializing QLKNN')
  end
  
  % Update QLKNN
  cmd = ['git -C ', submodulePath, ' submodule update QLKNN-fortran qlknn-hyper-namelists'];
  cmd_res = system(cmd);
  if cmd_res ~= 0
      error('Problem updating QLKNN')
  end
  
  % Initizalize build system
  cmd = ['git -C ' qlknnPath ' submodule init tubs'];
  cmd_res = system(cmd);
  if cmd_res ~= 0
    error('Problem initializing build system')
  end
  
  % Update build system
  cmd = ['git -C ' qlknnPath ' submodule update tubs'];
  cmd_res = system(cmd);
  if cmd_res ~= 0
    error('Problem updating build system')
  end
  
  success = qlknn_compile();
  if ~success
    error('Could not compile/MEX QLKNN')
  end
    
  qlknn_bin = fullfile(qlknnPath, 'bin');
  if ~exist(qlknn_bin, 'dir')
    error(['Could not find QLKNN bin dir here: ' qlknn_bin]) 
  elseif ~is_in_path(qlknn_bin)
    addpath(qlknn_bin)
  end

  % QLKNN-fortran has its own internal defaults. We copy that structure here
  qlknn_params.use_effective_diffusivity = true; % Use Gamma_e, not separate Ds and Vs
  qlknn_params.use_ion_diffusivity_networks = false; % Use D_e and V_e, not D_i and D_e (if using Ds and Vs)
  qlknn_params.apply_victor_rule = false; % Apply the 'victor rule' rotation model
  qlknn_params.calc_heat_transport = true; % Calculate q_e and q_i
  qlknn_params.calc_part_transport = true; % Calculate Gamma/D/V
  qlknn_params.use_ETG = true; % Calculate ETG mode transport
  qlknn_params.use_ITG = true; % Calculate ITG mode transport
  qlknn_params.use_TEM = true; % Calculate TEM mode transport

  qlknn_params.constrain_inputs = true; % Clip inputs to min_input and max_input, bound by margin_input
  qlknn_params.constrain_outputs = false; % Clip outputs to min_input and max_input, bound by margin_input

  qlknn_params.min_input = [1., 0., 0., -5., 0.66, -1., .09, 0.25, -5., -100., -100.];
  qlknn_params.max_input = [3., 14., 14., 6., 15., 5., .99, 2.5, 0., 100., 100.];
  qlknn_params.margin_input = 0.95;

  qlknn_params.min_output = -100;
  qlknn_params.max_output = 100;
  qlknn_params.margin_output = 1;
  
  qlknn_params.em_stab = 0; % Apply em stabilisation rule from JETTO. Currenty 0 (off) or 1 (reduce gradient by fast-ion pressure)
  
  qlknn_params.shear_minus_alpha = false;
  
  % Debugging parameters
  qlknn_params.apply_stability_clipping = true; % Clip stable modes, should always be true
  qlknn_params.merge_modes = true; % Merge ETG/ITG/TEM modes together to a single total flux
  qlknn_params.force_evaluate_all = false; % Debug flag: Force all networks to be evaluated
  qlknn_params.verbosity = 0;
  qlknn_params.channel_tag = 'qlknn'; % Tag in plots for this network
  
  qlknn_params = appendRAPTORInterfaceParameters10D(qlknn_params);

  varargout{2} = qlknn_params;
  varargout{1} = qlknn_network;
  return
end

if isfield(pp, 'norms')
    norms = pp.norms;
    opts = pp;
    opts = rmfield(opts, 'norms');
else
    if pp.apply_victor_rule
        error('pp.norms needs to be defined if applying victor rule')
    end
    opts = pp;
end
%% Assign outputs
if calcder
    if pp.apply_victor_rule
      [Nout, dNout_dNin] = qlknn_mex(net.namelist_dir, Nin, opts.verbosity, opts, norms);
    else
      [Nout, dNout_dNin] = qlknn_mex(net.namelist_dir, Nin, opts.verbosity, opts);
    end
else
    if pp.apply_victor_rule
      Nout = qlknn_mex(net.namelist_dir, Nin, opts.verbosity, opts, norms);
    else
      Nout = qlknn_mex(net.namelist_dir, Nin, opts.verbosity, opts);
    end
end

varargout{1} = Nout;
if calcder
    varargout{2} = dNout_dNin;
end


return

function is_in_path = is_in_path(Folder)
  pathCell = regexp(path, pathsep, 'split');
  if ispc  % Windows is not case-sensitive
    is_in_path = any(strcmpi(Folder, pathCell));
  else
    is_in_path = any(strcmp(Folder, pathCell));
  end
  return
