iOptimizerADMMLim.cc

Go to the documentation of this file.

#include "iOptimizerADMMLim.hh"
#include "sOutputManager.hh"

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================


iOptimizerADMMLim::iOptimizerADMMLim() : vOptimizer()
{
  // ---------------------------
  // Mandatory member parameters
  // ---------------------------

  // Initial value at 1
  m_initialValue = 1.;
  // Two additional backward images in addition to default one
  m_nbBackwardImages = 3;
  // ADMMLim accepts penalties
  m_requiredPenaltyDerivativesOrder = 1;
  // ADMMLim is compatible with histogram data
  m_listmodeCompatibility = false;
  m_histogramCompatibility = true;
  // ADMMLim is specific to emission data
  m_emissionCompatibility = true;
  m_transmissionCompatibility = false;

  // ADMMLim needs pre process loop to initialize v (compute v^0)
  m_needPreIteration = true;
  // ADMMLim needs 3 sub steps: gradient computation, image update using the gradient, v and u update
  m_nbSubSteps = 3;

  // --------------------------
  // Specific member parameters
  // --------------------------

  m4p_firstDerivativePenaltyImage = NULL;
  mp_outputIterations = NULL;

  m2p_rBackgroundEvents = NULL;
  m2p_yData = NULL;
  
  m_alpha = -1.;
  m_previousAlpha = -1.;
  m_tau = -1.;
  m_mu = -1.;
  m_tauMax = -1.;
  m_xi = -1.;
  
  m4p_currentGrad = NULL;
  m4p_currentGradFrameGates = NULL;
  m2p_toWrite_uk = NULL;
  m2p_toWrite_vk = NULL;
  m2p_vectorAx = NULL;

  mp_gradSquareNorm = NULL;
  mp_projGradSquareNorm = NULL;
  mp_penaltyGradSquareNorm = NULL;
  m_relativePrimalSquareNorm = -1.;
  m_relativeDualSquareNorm = -1.;
  mp_primalSquareNorm = NULL;
  mp_squareNorm_Ax = NULL;
  mp_squareNorm_v = NULL;

  m_adaptiveParameters = ADMMLIM_NOT_DEFINED;
  m_stoppingCriterionValue = -1.;
  m_saveSinogramsUAndV = false;
  m_criterionReachedIteration = -1.;
  m_criterionReachedIsPassed = false;
  m_firstIteration = -1;
  m_firstIterationIsPassed = false;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

iOptimizerADMMLim::~iOptimizerADMMLim()
{
  // Note: there is no need to deallocate the images themselves as they are allocate using the
  //       miscellaneous image function from the image space, which automatically deals with
  //       memory deallocations.
  // Delete the first order derivative penalty image and other variables indexed on basis functions
  if (m4p_firstDerivativePenaltyImage)
  {
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      if (m4p_firstDerivativePenaltyImage[tbf])
      {
        // Loop over respiratory basis functions
        for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
        {
          if (m4p_firstDerivativePenaltyImage[tbf][rbf])
          {
            free(m4p_firstDerivativePenaltyImage[tbf][rbf]);
          }
        }
        free(m4p_firstDerivativePenaltyImage[tbf]);
      }
    }
    free(m4p_firstDerivativePenaltyImage);
  }
  // Deallocate all pointers
  if (m2p_rBackgroundEvents)
  {
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      if (m2p_rBackgroundEvents[th])
      {
        free(m2p_rBackgroundEvents[th]);
      }
    }
    free(m2p_rBackgroundEvents);
  }
  if (m2p_yData)
  {
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      if (m2p_yData[th])
      {
        free(m2p_yData[th]);
      }
    }
    free(m2p_yData);
  }
  if (m4p_currentGrad)
  {
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      if (m4p_currentGrad[tbf])
      {
        // Loop over respiratory basis functions
        for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
        {
          if (m4p_currentGrad[tbf][rbf])
          {
            free(m4p_currentGrad[tbf][rbf]);
          }
        }
        free(m4p_currentGrad[tbf]);
      }
    }
    free(m4p_currentGrad);
  }
  if (m4p_currentGradFrameGates)
  {
    // Loop over time frames
    for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
    {
      if (m4p_currentGradFrameGates[fr])
      {
        // Loop over respiratory gates
        for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++)
        {
          if (m4p_currentGradFrameGates[fr][rg])
          {
            free(m4p_currentGradFrameGates[fr][rg]);
          }
        }
        free(m4p_currentGradFrameGates[fr]);
      }
    }
    free(m4p_currentGradFrameGates);
  }
  if (m2p_toWrite_uk)
  {
    free(m2p_toWrite_uk);
  }
  if (m2p_toWrite_vk)
  {
    free(m2p_toWrite_vk);
  }
  if (m2p_vectorAx)
  {
    free(m2p_vectorAx);
  }
  if (mp_gradSquareNorm)
  {
    free(mp_gradSquareNorm);
  }
  if (mp_projGradSquareNorm)
  {
    free(mp_projGradSquareNorm);
  }
  if (mp_penaltyGradSquareNorm)
  {
    free(mp_penaltyGradSquareNorm);
  }
  if (mp_primalSquareNorm)
  {
    free(mp_primalSquareNorm);
  }
  if (mp_squareNorm_Ax)
  {
    free(mp_squareNorm_Ax);
  }
  if (mp_squareNorm_v)
  {
    free(mp_squareNorm_v);
  }
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

void iOptimizerADMMLim::ShowHelpSpecific()
{
  cout << "This optimizer is the ADMM with non-negativity constraint on projection space." << endl;
  cout << "It was proposed by H. Lim, Y. K. Dewaraja, and J. A. Fessler, Physics " << endl;
  cout << "in Medicine & Biology, 2018. It is for now developed without subsets" << endl; 
  cout << "It is implemented using the same equations as in the original paper, also with one inner iteration." << endl;
  cout << "The algorithm can be used without penalty, with quadratic penalty or with Markov Random Field (MRF) penalty with quadratic potential." << endl;
  cout << "The following options can be used according to the paper by B. Wohlberg in 2017 on residual balancing:" << endl;
  cout << "  - alpha: to set the convergence speed of the ADMM with penalty parameter alpha." << endl;
  cout << "  - adaptive : to set which parameters are adaptive (must be set to nothing, alpha, or tau (which means alpha and tau))." << endl;
  cout << "  - mu: factor to balance primal and dual residual in adaptive alpha computation." << endl;
  cout << "  - tau: Factor to multiply alpha in adaptive alpha computation." << endl;
  cout << "  - xi: Factor to balance primal and dual residual convergence speed in adaptive tau computation." << endl;
  cout << "  - tau_max: Maximum value of the tau factor to multiply alpha in adaptive alpha and tau computation." << endl;
  cout << "  - stoppingCriterionValue: Error value for a stopping criterion based on small residuals. 0 means no stopping criterion." << endl;
  cout << "  - saveSinogramsUAndV: to save or not the computed sinograms u and v for same iterations as image " << endl;
  cout << "    (sinograms u and v can only be saved for non TOF data)" << endl;
  cout << "Adviced practice is to use adaptive tau with xi to 1 or to use adaptive alpha with mu to 10 and tau to 2." << endl;
  cout << "A stopping criterion can be used (both residuals smaller than 1e-4 for instance)." << endl;
  cout << "This optimizer was extended from Lim et al. paper to dynamic frame by frame reconstruction, and TOF reconstruction. However, " << endl;
  cout << "only one alpha value can be derived from ADMM as the likelihood is maximized over all the frames/TOF bins at the same time." << endl;
  cout << "Adviced practice is to launch a frame by frame reconstruction with one CASToR command line per frame, so that alpha is adapted to the frame." << endl;
  cout << "Moreover, the use of time basis functions may not make sense as, for some of them, the image can have negative voxels (as the optimizer" << endl;
  cout << "allows for it), which can be meaningless from a physiological or physical point of view." << endl; 
  cout << "TOF reconstruction is not expected to converge quicker than non-TOF reconstruction as alpha is related to the constraint Ax=v" << endl;
  cout << "which is harder to satisfy for several TOF bins." << endl;
  cout << "N.B.: The first iteration will not update the image if no penalty is used or if a uniform initialization image is used with MRF penalty." << endl;
  cout << "In these cases, it will only change sinograms u and v." << endl;}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::ReadConfigurationFile(const string& a_configurationFile)
{
  string key_word = "";
  string buffer = "";
  // Read the penalty parameter alpha
  key_word = "alpha";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_alpha, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read which parameters are adaptive
  key_word = "adaptive";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &buffer, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  if (buffer == "nothing")
  {
    m_adaptiveParameters = ADMMLIM_NOT_ADAPTIVE;
  }
  else if (buffer == "alpha")
  {
    m_adaptiveParameters = ADMMLIM_ADAPTIVE_ALPHA;
  }
  else if (buffer == "both")
  {
    m_adaptiveParameters = ADMMLIM_ADAPTIVE_ALPHA_AND_TAU;
  }
  // Read the factor mu
  key_word = "mu";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_mu, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the factor tau
  key_word = "tau";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_tau, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the factor xi
  key_word = "xi";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_xi, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the factor tau_max
  key_word = "tau_max";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_tauMax, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read the value of the stopping criterion
  key_word = "stoppingCriterionValue";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &m_stoppingCriterionValue, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  // Read if the sinograms u and v must be saved
  key_word = "saveSinogramsUAndV";
  if (ReadDataASCIIFile(a_configurationFile, key_word, &buffer, 1, KEYWORD_MANDATORY))
  {
    Cerr("***** iOptimizerADMMLim::ReadConfigurationFile() -> Failed to get the '" << key_word << "' keyword !" << endl);
    return 1;
  }
  if (buffer == "true")
  {
    m_saveSinogramsUAndV = ADMMLIM_SAVE_SINOGRAMS;
  }
  else if (buffer == "false")
  {
    m_saveSinogramsUAndV = ADMMLIM_DO_NOT_SAVE_SINOGRAMS;
  }
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::ReadOptionsList(const string& a_optionsList)
{
  // There are 8 floating point variables as option
  const int nb_options = 8;
  FLTNB options[nb_options];
  // Read them
  if (ReadStringOption(a_optionsList, options, nb_options, ",", "ADMMLim configuration"))
  {
    Cerr("***** iOptimizerADMMLim::ReadOptionsList() -> Failed to correctly read the list of options !" << endl);
    return 1;
  }
  // Affect options
  m_alpha = options[0];
  m_adaptiveParameters = options[1];
  m_mu = options[2];
  m_tau = options[3];
  m_xi = options[4];
  m_tauMax = options[5];
  m_stoppingCriterionValue = options[6];
  m_saveSinogramsUAndV = options[7];
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::CheckSpecificParameters()
{
  // Check adaptive parameters
  if (m_adaptiveParameters == ADMMLIM_NOT_DEFINED || (m_adaptiveParameters != ADMMLIM_NOT_ADAPTIVE && m_adaptiveParameters != ADMMLIM_ADAPTIVE_ALPHA && m_adaptiveParameters != ADMMLIM_ADAPTIVE_ALPHA_AND_TAU))
  {
    Cerr("***** iOptimizerADMMLim::CheckSpecificParameters() -> Should provide which parameters are adaptive (nothing, alpha or tau)" << endl);
    return 1;
  }
  // Check the user provided a boolean to save sinograms or not
  if (m_saveSinogramsUAndV != ADMMLIM_DO_NOT_SAVE_SINOGRAMS && m_saveSinogramsUAndV != ADMMLIM_SAVE_SINOGRAMS)
  {
    Cerr("***** iOptimizerADMMLim::CheckSpecificParameters() -> Please specify if sinograms u and v need to be stored (put 0 or 1)" << endl);
    return 1;
  }
  //
  if (m_saveSinogramsUAndV && !(mp_ImageDimensionsAndQuantification->GetNbTimeFrames() == 1 && mp_ImageDimensionsAndQuantification->GetNbRespGates() == 1 && mp_ImageDimensionsAndQuantification->GetNbCardGates() == 1))
  {
    Cerr("!!!!! iOptimizerADMMLim::CheckSpecificParameters() -> There are several frames or gates, thus sinograms u and v will not be stored even if asked by the user" << endl);
  }
  // Check that penalty parameter alpha is positive
  if (m_alpha<=0.)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Provided alpha (" << m_alpha << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that factor mu is positive
  if (m_mu<=0.)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Provided mu (" << m_mu << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that factor tau is greater than 1
  if (m_tau<=1.)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Provided tau (" << m_tau << ") must be strictly greater than 1 !" << endl);
    return 1;
  }
  // Check that factor xi is positive
  if (m_xi<=0.)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Provided xi (" << m_xi << ") must be strictly positive !" << endl);
    return 1;
  }
  // Check that factor tau_max is greater than 1
  if (m_tauMax<=1.)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Provided tau (" << m_tauMax << ") must be strictly greater than 1 !" << endl);
    return 1;
  }
    // Check that value of the stopping criterion is positive
  if (m_stoppingCriterionValue<0.)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Provided value for the stopping criterion (" << m_stoppingCriterionValue << ") must be positive !" << endl);
    return 1;
  }
  // Check that number of additional data is 0 (if no initialization of v and u) or 2 (if user initialized v and u).
  if (mp_DataFile->GetNbAdditionalData()==2)
  {
    Cerr("!!!!! iOptimizerADMMLim::CheckSpecificParameters() -> ADMMLim does not require additional data, except if you want to initialize u^k and v^k !" << endl);
  }
  else if (mp_DataFile->GetNbAdditionalData()!=0 && mp_DataFile->GetNbAdditionalData()!=2)
  {
    Cerr("***** iOptimizerADMMLim::CheckSpecificParameters() -> ADMMLim requires 2 additional data for u^k and v^k if they need to be initialized, or nothing !" << endl);
    return 1;
  }
  // Cannot deal with list-mode or transmission data
  if (m_dataSpec==SPEC_TRANSMISSION || m_dataMode==MODE_LIST)
  {
    Cerr("***** iOptimizerADMMLim->CheckSpecificParameters() -> Cannot reconstruct list-mode or transmission data !" << endl);
    return 1;
  }
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::InitializeSpecific()
{
  // -------------------------------------------------------------------
  // First, check if a Penalty has been initialized
  if (mp_Penalty != NULL)
  {
    // Check that Penalty is MRF
    if (mp_Penalty->GetPenaltyID() == "MRF")
    {
      // If Penalty is MRF, now check potential function is quadratic
      if ((dynamic_cast<iPenaltyMarkovRandomField*>(mp_Penalty))->GetPotentialType() != MRF_POTENTIAL_QUADRATIC)
      {
        Cerr("***** iOptimizerADMMLim->InitializeSpecific() -> This optimizer is only compatible with Quadratic potential function for the Markov Random Field penalty (MRF) penalty !" << endl);
        Cerr("                                                             Please use 'potential function: quadratic' in the penalty initialization file!" << endl);
        return 1;
      }
    }
    else if (mp_Penalty->GetPenaltyID() != "QUAD")
    {  
      Cerr("***** iOptimizerADMMLim->InitializeSpecific() -> This optimizer is only compatible with the Markov Random Field (MRF) penalty, with quadratic penalty, or without any penalty !" << endl);
      return 1;
    }  
  }
  // Allocate and create the penalty image
  m4p_firstDerivativePenaltyImage = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()*sizeof(FLTNB***));
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  { 
    m4p_firstDerivativePenaltyImage[tbf] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB**));
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      m4p_firstDerivativePenaltyImage[tbf][rbf] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions()*sizeof(FLTNB*));
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        // Get a pointer to a newly allocated image
        m4p_firstDerivativePenaltyImage[tbf][rbf][cbf] = mp_ImageSpace -> AllocateMiscellaneousImage();
        // Initialize to 0, in case the penalty is not used
        for (INTNB v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++) m4p_firstDerivativePenaltyImage[tbf][rbf][cbf][v] = 0.;
      }
    }
  }
  
  // Allocate and initialize useful variables for u and v computation (as many as threads and as many as TOF bins)
  m2p_rBackgroundEvents = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(FLTNB*));
  m2p_yData = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(FLTNB*));
  for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
  {
    m2p_rBackgroundEvents[th] = (FLTNB*)malloc(m_nbTOFBins*sizeof(FLTNB));
    m2p_yData[th] = (FLTNB*)malloc(m_nbTOFBins*sizeof(FLTNB));
    for (int b=0; b<m_nbTOFBins; b++)
    {
      m2p_rBackgroundEvents[th][b] = 0.;
      m2p_yData[th][b] = 0.;
    }
  }

  // Allocate norms related to gradient and penalty strength (as many as threads)
  mp_gradSquareNorm = (HPFLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(HPFLTNB));
  mp_projGradSquareNorm = (HPFLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(HPFLTNB));
  mp_penaltyGradSquareNorm = (HPFLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(HPFLTNB));
  mp_primalSquareNorm = (HPFLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(HPFLTNB));
  mp_squareNorm_Ax = (HPFLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(HPFLTNB));
  mp_squareNorm_v = (HPFLTNB*)malloc(mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection()*sizeof(HPFLTNB));

  // Allocate and create the whole gradient at previous and current iteration, which need to be stored to compute forward projection and norms of them
  m4p_currentGrad = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()*sizeof(FLTNB***));
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  { 
    m4p_currentGrad[tbf] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()*sizeof(FLTNB**));
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      m4p_currentGrad[tbf][rbf] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions()*sizeof(FLTNB*));
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        // Get a pointer to a newly allocated image
        m4p_currentGrad[tbf][rbf][cbf] = mp_ImageSpace -> AllocateMiscellaneousImage();
        // Initialize to 0
        for (INTNB v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++) m4p_currentGrad[tbf][rbf][cbf][v] = 0.;
      }
    }
  }

  // Allocate gradient indexed on frames and gates
  m4p_currentGradFrameGates = (FLTNB****)malloc(mp_ImageDimensionsAndQuantification->GetNbTimeFrames()*sizeof(FLTNB***));
  // Loop over time frames
  for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
  { 
    m4p_currentGradFrameGates[fr] = (FLTNB***)malloc(mp_ImageDimensionsAndQuantification->GetNbRespGates()*sizeof(FLTNB**));
    // Loop over respiratory gates
    for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++)
    {
      m4p_currentGradFrameGates[fr][rg] = (FLTNB**)malloc(mp_ImageDimensionsAndQuantification->GetNbCardGates()*sizeof(FLTNB*));
      // Loop over cardiac gates
      for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++)
      {
        // Get a pointer to a newly allocated image
        m4p_currentGradFrameGates[fr][rg][cg] = mp_ImageSpace -> AllocateMiscellaneousImage();
      }
    }
  }
  //*/

  // Allocate and create the sinograms u and v if no are provided. Otherwise use pointers to additional data memory space
  m2p_toWrite_uk = (FLTNB**)malloc(m_nbTOFBins*sizeof(FLTNB*));
  m2p_toWrite_vk = (FLTNB**)malloc(m_nbTOFBins*sizeof(FLTNB*));
  // Check if sinograms u and v were given by user
  if (mp_DataFile->GetNbAdditionalData()==2)  u_and_v_from_user = true;
  else  u_and_v_from_user = false;
  // Loop on TOF bins
  for (int b=0; b<m_nbTOFBins; b++)
  {
    if (u_and_v_from_user)
    {
      m2p_toWrite_uk[b] = mp_DataFile->m2p_additionalData[0];
      m2p_toWrite_vk[b] = mp_DataFile->m2p_additionalData[1];
    }
    else
    {
      m2p_toWrite_uk[b] = mp_DataFile->GetNewAdditionalDataMatrix(1);
      m2p_toWrite_vk[b] = mp_DataFile->GetNewAdditionalDataMatrix(1);
    }
  }
  // Allocate and create the gradient projection for each TOF bin
  m2p_vectorAx = (FLTNB**)malloc(m_nbTOFBins*sizeof(FLTNB*));
  for (int b=0; b<m_nbTOFBins; b++)
  {
    m2p_vectorAx[b] = mp_DataFile->GetNewAdditionalDataMatrix(1);
  }

  // Verbose
  if (m_verbose>=VERBOSE_NORMAL)
  {
    Cout("iOptimizerADMMLim::InitializeSpecific() -> Use the ADMMLim optimizer" << endl);
    if (m_verbose>=VERBOSE_DETAIL)
    {
      Cout("  --> Initial image value: " << m_initialValue << endl);
      if (m_adaptiveParameters == ADMMLIM_NOT_ADAPTIVE)
      {
        Cout("  --> Adaptive mode: " << "not adaptive" << endl);
      }
      else if (m_adaptiveParameters == ADMMLIM_ADAPTIVE_ALPHA)
      {
        Cout("  --> Adaptive mode: " << "adaptive alpha" << endl);
      }
      else if (m_adaptiveParameters == ADMMLIM_ADAPTIVE_ALPHA_AND_TAU)
      {
        Cout("  --> Adaptive mode: " << "adaptive alpha and tau" << endl);
      }
      Cout("  --> Initial penalty parameter alpha: " << m_alpha << endl);
      Cout("  --> Initial penalty parameter tau: " << m_tau << endl);
      if (m_adaptiveParameters != ADMMLIM_NOT_ADAPTIVE)
      {
        Cout("  --> Factor mu: " << m_mu << endl);
        Cout("  --> Factor xi: " << m_xi << endl); 
      }
      if (m_adaptiveParameters == ADMMLIM_ADAPTIVE_ALPHA_AND_TAU)
      {
        Cout("  --> Maximum tau value (if adaptive alpha and tau) tau_max: " << m_tau << endl);
      }
      if (m_stoppingCriterionValue > 0.)
      {
        Cout("  --> Stopping criterion is used" << endl);
      }
    }
  }
  // Normal end
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::PreDataUpdateSpecificStep()
{
  // Check the user asked only for one subset (check only first subset of mp_nbSubsets for the sake of simplicity)
  // CASToR's structure does not enable to access the number of subsets in CheckSpecificParameters(), so we do it here
  if (mp_nbSubsets[0] != 1)
  {
    Cerr("***** iOptimizerADMMLim::PreDataUpdateSpecificStep() -> Should provide only one subset !" << endl);
    // Leave computation
    return 1;
  }
  // Initialize previous value of alpha to alpha
  if (m_isInPreIteration) m_previousAlpha = m_alpha;
  // Initialize norm of gradient
  if (m_currentSubStep == 0)
  {
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      mp_gradSquareNorm[th] = 0.;
    }
  }
  // Initialize norms for relative primal residual computation
  if (m_isInPreIteration || m_currentSubStep == 2) 
  if (m_currentSubStep == 2) 
  {
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      mp_primalSquareNorm[th] = 0.;
      mp_squareNorm_Ax[th] = 0.;
      mp_squareNorm_v[th] = 0.;
    }
  }

  // Check if the stopping criterion is reached
  if (sqrt(m_relativePrimalSquareNorm) < m_stoppingCriterionValue && sqrt(m_relativeDualSquareNorm) < m_stoppingCriterionValue && !m_criterionReachedIsPassed)
  {
    m_criterionReachedIteration = m_currentIteration - 1;
    m_criterionReachedIsPassed = true;
  }
  // Write file with residuals and iteration and stop the computation if stopping criterion is reached
  if (m_currentIteration == m_criterionReachedIteration + 1 && m_criterionReachedIteration != -1)
  {
    // Get the path for the file for the stopping criterion
    sOutputManager* p_outputManager = sOutputManager::GetInstance();
    string temps_ss_stop = p_outputManager->GetPathName() + p_outputManager->GetBaseName();
    // Show both relative residuals and the final iteration
    temps_ss_stop +=  "_adaptive_stopping_criteria.log";
    ofstream outfile;
    outfile.open(temps_ss_stop);
    outfile << setprecision(numeric_limits<FLTNB>::digits10 + 1);
    outfile << "final outer iteration        : " << endl;
    outfile << m_currentIteration << endl;
    outfile << "relPrimal      : " << endl;
    outfile << sqrt(m_relativePrimalSquareNorm) << endl;
    outfile << "relDual        : " << endl;
    outfile << sqrt(m_relativeDualSquareNorm) << endl;
    outfile.close();
    // Leave computation
    Cerr("!!!!! Stopping criterion reached ! Leaving computation." << endl);
    return 1;
  }
  // Check if we are in first computed iteration, will be used to initialize sinogram u and v if user provieded some
  if (!m_isInPreIteration && !m_firstIterationIsPassed) 
  {
    m_firstIteration = m_currentIteration;
    m_firstIterationIsPassed = true;
  }

  if (m_currentSubStep == 1)
  {
    // Loop over time frames
    for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
    {
      // Loop over respiratory gates
      for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++)
      {
        // Loop over cardiac gates
        for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++)
        {
          // Voxel index
          INTNB v;
          // Multi-threading over voxels
          #pragma omp parallel for private(v) schedule(guided)
          for (v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
          { 
            // Clean buffer (gradient over frames and gates)
            m4p_currentGradFrameGates[fr][rg][cg][v] = 0.;
          }
        }
      }
    }
  }

  if (m_currentSubStep == 1) 
  {
    // Initialize norm of gradient projection before x update
    for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      mp_projGradSquareNorm[th] = 0.;
    }
    // Apply convolver to gradient image because image-based PSF not inside ForwardProject function
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      // Loop over respiratory basis functions
      for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
      {
        // Loop over cardiac basis functions
        for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
        {
          if (mp_ImageConvolverManager->ApplyConvolution(m4p_currentGrad[tbf][rbf][cbf]))
          {
            Cerr("***** iOptimizerADMMLim::PreDataUpdateSpecificStep() -> A problem occurred while applying image convolver to gradient image !" << endl);
            return 1;
          }
        }
      }
    }

    // Loop over time frames
    for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
    {
      // Loop over respiratory gates
      for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++)
      {
        // Loop over cardiac gates
        for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++)
        {
          // Loop over time basis functions
          for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
          {
            // Get frame/basis coefficient and continue if null
            FLTNB time_basis_coef = mp_ImageDimensionsAndQuantification->GetTimeBasisCoefficient(tbf,fr);
            if (time_basis_coef==0.) continue;
            // Loop over respiratory basis functions
            for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
            {
              // Get resp_gate/basis coefficient and continue if null
              FLTNB resp_basis_coef = mp_ImageDimensionsAndQuantification->GetRespBasisCoefficient(rbf,rg);
              if (resp_basis_coef==0.) continue;
              // Loop over cardiac basis functions
              for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
              {
                // Get card_gate_basis coefficient and continue if null
                FLTNB card_basis_coef = mp_ImageDimensionsAndQuantification->GetCardBasisCoefficient(cbf,cg);
                if (card_basis_coef==0.) continue;
                // Compute global coefficient
                FLTNB global_basis_coef = time_basis_coef * resp_basis_coef * card_basis_coef;
                
                // Voxel index
                INTNB v;
                // Multi-threading over voxels
                #pragma omp parallel for private(v) schedule(guided)
                for (v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
                { 
                  // Apply dynamic coefficient conversion
                  m4p_currentGradFrameGates[fr][rg][cg][v] += m4p_currentGrad[tbf][rbf][cbf][v] * global_basis_coef;
                }
              }
            }
          }
        }
      }
    }
  
    // If penalty is MRF with quadratic potential, compute finite difference matrix times gradient vector : Cg
    // Its norm is needed for stepsize computation to update image x with quadratic MRF penalty
    if (mp_Penalty != NULL)
    {
      if (mp_Penalty->GetPenaltyID() == "MRF")
      {
        // Initialize norm of penalty gradient
        for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
        {
          mp_penaltyGradSquareNorm[th] = 0.;
        }
        // Loop over time frames
        for (int fr=0; fr<mp_ImageDimensionsAndQuantification->GetNbTimeFrames(); fr++)
        {
          // Loop over respiratory gates
          for (int rg=0; rg<mp_ImageDimensionsAndQuantification->GetNbRespGates(); rg++)
          {
            // Loop over cardiac gates
            for (int cg=0; cg<mp_ImageDimensionsAndQuantification->GetNbCardGates(); cg++)
            {
              // Loop over time basis functions
              for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
              {
                // Get frame/basis coefficient and continue if null
                FLTNB time_basis_coef = mp_ImageDimensionsAndQuantification->GetTimeBasisCoefficient(tbf,fr);
                if (time_basis_coef==0.) continue;
                // Loop over respiratory basis functions
                for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
                {
                  // Get resp_gate/basis coefficient and continue if null
                  FLTNB resp_basis_coef = mp_ImageDimensionsAndQuantification->GetRespBasisCoefficient(rbf,rg);
                  if (resp_basis_coef==0.) continue;
                  // Loop over cardiac basis functions
                  for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
                  {
                    // Get card_gate_basis coefficient and continue if null
                    FLTNB card_basis_coef = mp_ImageDimensionsAndQuantification->GetCardBasisCoefficient(cbf,cg);
                    if (card_basis_coef==0.) continue;
                    // Compute global coefficient
                    FLTNB global_basis_coef = time_basis_coef * resp_basis_coef * card_basis_coef;
                    // In order to detect problems in the multi-threaded loop
                    bool problem = false;
                    // Voxel index
                    INTNB v;
                    // Multi-threading over voxels
                    #pragma omp parallel for private(v) schedule(guided)
                    for (v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
                    { 
                      // Get the thread index
                      int th = 0;
                      #ifdef CASTOR_OMP
                      th = omp_get_thread_num();
                      #endif
                      // Local precomputation step if needed by the penalty
                      if (mp_Penalty->LocalPreProcessingStep(tbf,rbf,cbf,v,th))
                      {
                        Cerr("***** iOptimizerADMMLim::PreDataUpdateSpecificStep() -> A problem occurred while computing the penalty local pre-processing step for voxel " << v << " !" << endl);
                        problem = true;
                      }
                      // Compute finite difference matrix times gradient vector Cg (its norm will be used in the image update step).
                      // We need to loop over basis functions to do the local pre processing step in MRF penalty
                      // Apply dynamic coefficient conversion
                      mp_penaltyGradSquareNorm[th] += 2. * (HPFLTNB)mp_Penalty->ComputePenaltyValue(m4p_currentGrad[tbf][rbf][cbf],v,th) * global_basis_coef;
                      //mp_penaltyGradSquareNorm[fr][rg][cg] += 2. * (HPFLTNB)mp_Penalty->ComputePenaltyValue(m4p_currentGradFrameGates[fr][rg][cg],v,th) * global_basis_coef;
                    }
                    // Check for problems
                    if (problem)
                    {
                      Cerr("***** iOptimizerADMMLim::PreImageUpdateSpecificStep() -> A problem occurred inside the multi-threaded loop, stop now !" << endl);
                      return 1;
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  // Normal end  
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::DataStep4Optional( oProjectionLine* ap_Line, vEvent* ap_Event,
                                   int a_bed, int a_timeFrame, int a_respGate, int a_cardGate,
                                   int a_th )
{
  // Loop on TOF bins
  for (int b=0; b<ap_Line->GetNbTOFBins(); b++)
  {
    // Set the current TOF bin of the projection-line
    ap_Line->SetCurrentTOFBin(b);
    // Scale u if alpha has just changed (cf Boyd et al. 2011 paper on ADMM, Foundations and Trends in Machine Learning, paragraph 3.4.1)
    if (m_currentSubStep == 0)
    {
      FLTNB coeff_alpha = m_alpha / m_previousAlpha;
      m2p_toWrite_uk[b][ap_Line->GetEventIndex()] /= coeff_alpha;
    }
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::SensitivitySpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_weight,
                                                   FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue,
                                                   FLTNB a_quantificationFactor, oProjectionLine* ap_Line )
{
  // Line weight here is simply 1
  *ap_weight = 1;
  // That's all
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::DataStep5ComputeCorrections( oProjectionLine* ap_Line, vEvent* ap_Event,
                                             int a_bed, int a_timeFrame, int a_respGate, int a_cardGate,
                                             int a_th )
{
  // Override dataStep5 to get Ax, y, r_bar (depending on thread) for u and v computation
  // Loop on TOF bins
  for (int b=0; b<ap_Line->GetNbTOFBins(); b++)
  {
    // Set the current TOF bin of the projection-line
    ap_Line->SetCurrentTOFBin(b);
    // Compute Ax = (Ax+r_bar) - r_bar after x update, and store data and background events without backprojection for v computation
    if (m_isInPreIteration || m_currentSubStep == 2)
    {
      m2p_yData[a_th][b] = (FLTNB)ap_Event->GetEventValue(b);
      // Truncate data to 0 if negative
      if (m2p_yData[a_th][b]<0.) m2p_yData[a_th][b] = 0.;
      m2p_rBackgroundEvents[a_th][b] = (FLTNB)ap_Event->GetAdditiveCorrections(b)*mp_ImageDimensionsAndQuantification->GetFrameDurationInSec(a_bed, a_timeFrame);
      m2p_vectorAx[b][ap_Line->GetEventIndex()] = (HPFLTNB)m2p_forwardValues[a_th][b] - (HPFLTNB)m2p_rBackgroundEvents[a_th][b];
    }

    // Backward project (Ax - v + u) to get the gradient
    if (m_currentSubStep == 0 && !m_isInPreIteration) m3p_backwardValues[a_th][b][0] = m2p_vectorAx[b][ap_Line->GetEventIndex()] - (HPFLTNB)m2p_toWrite_vk[b][ap_Line->GetEventIndex()] + (HPFLTNB)m2p_toWrite_uk[b][ap_Line->GetEventIndex()];
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::DataSpaceSpecificOperations( FLTNB a_data, FLTNB a_forwardModel, FLTNB* ap_backwardValues,
                                                 FLTNB a_multiplicativeCorrections, FLTNB a_additiveCorrections, FLTNB a_blankValue,
                                                 FLTNB a_quantificationFactor, oProjectionLine* ap_Line )
{
  // Not useful in this optimizer, this function is directly written in DataStep5ComputeCorrections function as used variables depend on the current thread
  // Still need to define it as it is a pure virtual function in vOptimizer
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::DataStep6Optional( oProjectionLine* ap_Line, vEvent* ap_Event,
                                   int a_bed, int a_timeFrame, int a_respGate, int a_cardGate,
                                   int a_th )
{
  // Initialize v^0 and u^0 if optimizer is at pre iteration and no additional data was given (using Ax computed in dataStep5)
  // Loop on TOF bins
  for (int b=0; b<ap_Line->GetNbTOFBins(); b++)
  {
    // Set the current TOF bin of the projection-line
    ap_Line->SetCurrentTOFBin(b);
    // Initialize v^0 and u^0 if optimizer is at pre iteration and no additional data was given (using Ax computed in dataStep5)
    if (m_isInPreIteration && !u_and_v_from_user)
    {
      m2p_toWrite_vk[b][ap_Line->GetEventIndex()] = (FLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()];
      m2p_toWrite_uk[b][ap_Line->GetEventIndex()] = 0.;
    }
  }
  // Compute norm of projection of gradient for stepsize in x update
  HPFLTNB proj_grad_i_b = 0.;
  // Loop over TOF bins
  for (int b=0; b<ap_Line->GetNbTOFBins(); b++)
  {
    // Set the current TOF bin
    ap_Line->SetCurrentTOFBin(b);
    // Project the gradient and apply dynamic coefficient conversion
    if (m_currentSubStep == 1)
    {
      // Compute norm of projection of gradient for stepsize in x update
      proj_grad_i_b = (HPFLTNB)ForwardProject(ap_Line,m4p_currentGradFrameGates[a_timeFrame][a_respGate][a_cardGate]);
      mp_projGradSquareNorm[a_th] += proj_grad_i_b * proj_grad_i_b;
    }
    // Compute v and u after x update
    else if (m_currentSubStep == 2 && !m_isInPostIteration)
    {
      FLTNB previous_vk_current_event = m2p_toWrite_vk[b][ap_Line->GetEventIndex()];

      FLTNB beta = 0.5 * (1 / m_alpha + m2p_rBackgroundEvents[a_th][b] - m2p_toWrite_uk[b][ap_Line->GetEventIndex()] - (FLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()]);
      FLTNB gamma = m2p_rBackgroundEvents[a_th][b] * (m2p_toWrite_uk[b][ap_Line->GetEventIndex()] + (FLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()]) - (m2p_rBackgroundEvents[a_th][b] - m2p_yData[a_th][b]) / m_alpha;
      FLTNB v_hat = 0.;
      FLTNB beta_square_gamma = beta * beta + gamma;
      // Compute v_hat solution given the value of y
      if (m2p_yData[a_th][b] == 0.) 
      {
        v_hat = (FLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()] + m2p_toWrite_uk[b][ap_Line->GetEventIndex()] - 1 / m_alpha;
      }
      else
      {
        // Check for numerical instabilities. beta*beta+gamma should theoretically be positive. Set it to 0 otherwise
        if (beta_square_gamma < 0.)
        {
          Cerr("!!!!! beta * beta + gamma = " << beta_square_gamma << " , should not be negative ! Perhaps it is a case of divergence." << endl);
          beta_square_gamma = 0.;
        }
        // Compute v_hat using formula proposed in Lim et al. to avoid numerical instability of the quadratic solution
        v_hat = (beta < 0.) ? sqrt(beta_square_gamma) - beta :gamma / (sqrt(beta_square_gamma) + beta);
      }
      // Correct v if it does not respect the constraint (Ax+r>0) and store its value in sinogram
      m2p_toWrite_vk[b][ap_Line->GetEventIndex()] = ((v_hat + m2p_rBackgroundEvents[a_th][b]) > 0.) ? v_hat :  -m2p_rBackgroundEvents[a_th][b];

      // Update sinogram u (u^{k+1} := u^k + Ax^{k+1} - v^{k+1}) and store it
      m2p_toWrite_uk[b][ap_Line->GetEventIndex()] = m2p_toWrite_uk[b][ap_Line->GetEventIndex()] + (FLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()] - m2p_toWrite_vk[b][ap_Line->GetEventIndex()];

      // Store sinograms u^{k+1} and v^{k+1}-v^k to be backward projected for residuals computation
      m3p_backwardValues[a_th][b][1] = m2p_toWrite_uk[b][ap_Line->GetEventIndex()];
      m3p_backwardValues[a_th][b][2] = m2p_toWrite_vk[b][ap_Line->GetEventIndex()] - previous_vk_current_event;

      HPFLTNB diff = (HPFLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()] - (HPFLTNB)m2p_toWrite_vk[b][ap_Line->GetEventIndex()];
      mp_primalSquareNorm[a_th] += diff * diff;
      mp_squareNorm_Ax[a_th] += (HPFLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()]*(HPFLTNB)m2p_vectorAx[b][ap_Line->GetEventIndex()];
      mp_squareNorm_v[a_th] += (HPFLTNB)m2p_toWrite_vk[b][ap_Line->GetEventIndex()]*(HPFLTNB)m2p_toWrite_vk[b][ap_Line->GetEventIndex()];
    }
  }
  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::PreImageUpdateSpecificStep()
{
  // Merge multi-threads results to get final norms
  if (m_currentSubStep == 1)
  {
    for (int th=1; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      mp_gradSquareNorm[0] += mp_gradSquareNorm[th];
      mp_projGradSquareNorm[0] += mp_projGradSquareNorm[th];
      mp_penaltyGradSquareNorm[0] += mp_penaltyGradSquareNorm[th];
    }
  }
  else if (m_currentSubStep == 2)
  { 
    for (int th=1; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForProjection(); th++)
    {
      mp_primalSquareNorm[0] += mp_primalSquareNorm[th];
      mp_squareNorm_Ax[0] += mp_squareNorm_Ax[th];
      mp_squareNorm_v[0] += mp_squareNorm_v[th];
    }
  }

  // Residuals computation
  if (m_currentSubStep == 2)
  {
    int no_thread = 0;
    // Find the biggest norm between the 2 previously computed
    FLTNB primalMax = (mp_squareNorm_Ax[no_thread] > mp_squareNorm_v[no_thread]) ? mp_squareNorm_Ax[no_thread] : mp_squareNorm_v[no_thread];
    // Norm of relative primal residual
    m_relativePrimalSquareNorm = mp_primalSquareNorm[no_thread] / primalMax;
    FLTNB dual_norm = 0.;
    FLTNB AtuSquareNorm = 0.;
    FLTNB AtvvSquareNorm = 0.;
    // Loop over time basis functions
    for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
    {
      // Loop over respiratory basis functions
      for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
      {
        // Loop over cardiac basis functions
        for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
        {
          int no_thread = 0.;
          // Loop over voxels
          for (int v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
          {
            // Calculate norm of Atu and Atvv
            AtuSquareNorm += mp_ImageSpace->m6p_backwardImage[1][no_thread][tbf][rbf][cbf][v]*mp_ImageSpace->m6p_backwardImage[1][no_thread][tbf][rbf][cbf][v];
            AtvvSquareNorm += mp_ImageSpace->m6p_backwardImage[2][no_thread][tbf][rbf][cbf][v]*mp_ImageSpace->m6p_backwardImage[2][no_thread][tbf][rbf][cbf][v];
          }
        }
      }
    }
    // Compute norm of dual residual and relative one
    dual_norm = m_alpha * sqrt(AtvvSquareNorm);
    m_relativeDualSquareNorm = AtvvSquareNorm / AtuSquareNorm;

    // Define and compute relative dual residual which will be computed for current frame and gates
    if (m_adaptiveParameters != ADMMLIM_NOT_ADAPTIVE) 
    {
      if (m_adaptiveParameters == ADMMLIM_ADAPTIVE_ALPHA_AND_TAU)
      {
        // Compute adaptive tau (if asked for) to compute adaptive alpha for next iteration
        FLTNB adaptiveTau = sqrt((1/m_xi)*sqrt(sqrt(m_relativePrimalSquareNorm)/ sqrt(m_relativeDualSquareNorm)));
        if (adaptiveTau >= 1 && adaptiveTau < m_tauMax)
        {
          m_tau = adaptiveTau;
        }
        else if (adaptiveTau >= (1/m_tauMax) && adaptiveTau < 1)
        {
          m_tau = 1/adaptiveTau;
        }
        else
        {
          m_tau = m_tauMax;
        }
      }
      // Store current alpha value
      m_previousAlpha = m_alpha;
      // Compute adaptive alpha for next iteration (if asked for)
      if (m_adaptiveParameters != ADMMLIM_NOT_ADAPTIVE)
      {
        if (sqrt(m_relativePrimalSquareNorm) > m_xi*m_mu*sqrt(m_relativeDualSquareNorm))
        {
          m_alpha *= m_tau;
        }
        else if (sqrt(m_relativeDualSquareNorm) > (1/m_xi)*m_mu*sqrt(m_relativePrimalSquareNorm))
        {
          m_alpha /= m_tau;
        } 
        else
        {
          m_alpha = m_alpha;
        }
      }
    }
  
    if (m_optimizerFOMFlag)
    {
      // Get the path for the FOMs file
      sOutputManager* p_outputManager = sOutputManager::GetInstance();
      string temps_ss_alpha = p_outputManager->GetPathName() + p_outputManager->GetBaseName();
      // Show both alpha and the obtained adaptive alpha for next iteration, and residuals
      temps_ss_alpha +=  "_adaptive_it" + std::to_string(m_currentIteration + 1) + ".log";
      ofstream outfile;
      outfile.open(temps_ss_alpha);
      outfile << setprecision(numeric_limits<FLTNB>::digits10 + 1);
      outfile << "adaptive alpha : " << endl;
      outfile << m_alpha << endl;
      outfile << "adaptive tau" << endl;
      outfile << m_tau << endl;
      outfile << "relPrimal" << endl;
      outfile << sqrt(m_relativePrimalSquareNorm) << endl;
      outfile << "relDual" << endl;
      outfile << sqrt(m_relativeDualSquareNorm) << endl;
      outfile << "primal residual" << endl;
      int no_thread = 0;
      outfile << sqrt(mp_primalSquareNorm[no_thread]) << endl;
      outfile << "dual residual" << endl;
      outfile << dual_norm << endl;
      outfile.close();
    }

    // Compute true iteration (without subiterations)
    sOutputManager* p_outputManager = sOutputManager::GetInstance();
    if (m_saveSinogramsUAndV == ADMMLIM_SAVE_SINOGRAMS && m_nbTOFBins == 1)
    {
      if (mp_outputIterations[m_currentTotalSubIteration])
      {
        int first_TOF_bin = 0; 
        if (m2p_toWrite_vk)
        {
          // Write sinogram v^{k+1} with provided directory and filename
          stringstream temp_ss_v;
          temp_ss_v << p_outputManager->GetPathName() << p_outputManager->GetBaseName() << "_v_it" << m_currentIteration+1 <<  ".img";
          string a_pathToImg_v = temp_ss_v.str();
          IntfWriteImage(a_pathToImg_v, m2p_toWrite_vk[first_TOF_bin], mp_DataFile->GetNbEvents(), m_verbose);
        }
        if (m2p_toWrite_uk)
        {
          // Write sinogram u^{k+1} with provided directory and filename, to be used in next Python iteration
          stringstream temp_ss_u;
          temp_ss_u << p_outputManager->GetPathName() << p_outputManager->GetBaseName() << "_u_it" << m_currentIteration+1 <<  ".img";
          string a_pathToImg_u = temp_ss_u.str();
          IntfWriteImage(a_pathToImg_u, m2p_toWrite_uk[first_TOF_bin], mp_DataFile->GetNbEvents(), m_verbose);
        }
      }
    }


  }

  // ==========================================================================================
  // If no penalty, then exit (the penalty image term has been initialized to 0)
  if (mp_Penalty==NULL) return 0;
  // Set the number of threads
  #ifdef CASTOR_OMP
  omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation());
  #endif
  // Verbose
  if (m_verbose>=1) Cout("iOptimizerADMMLim::PreImageUpdateSpecificStep() -> Compute penalty term" << endl);
  // ==========================================================================================
  // Global precomputation step if needed by the penalty
  if (mp_Penalty->GlobalPreProcessingStep())
  {
    Cerr("***** iOptimizerADMMLim::PreImageUpdateSpecificStep() -> A problem occurred while computing the penalty pre-processing step !" << endl);
    return 1;
  }
  // ==========================================================================================
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  {
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        // In order to detect problems in the multi-threaded loop
        bool problem = false;
        // Voxel index
        INTNB v;
        // Multi-threading over voxels
        #pragma omp parallel for private(v) schedule(guided)
        for (v=0; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ(); v++)
        { 
          // Get the thread index
          int th = 0;
          #ifdef CASTOR_OMP
          th = omp_get_thread_num();
          #endif
          // Local precomputation step if needed by the penalty
          if (mp_Penalty->LocalPreProcessingStep(tbf,rbf,cbf,v,th))
          {
            Cerr("***** iOptimizerADMMLim::PreImageUpdateSpecificStep() -> A problem occurred while computing the penalty local pre-processing step for voxel " << v << " !" << endl);
            problem = true;
          }
          // Compute first derivative order penalty terms
          m4p_firstDerivativePenaltyImage[tbf][rbf][cbf][v] = mp_Penalty->ComputeFirstDerivative(mp_ImageSpace->m4p_image[tbf][cbf][rbf],v,th);
        }
        // Check for problems
        if (problem)
        {
          Cerr("***** iOptimizerADMMLim::PreImageUpdateSpecificStep() -> A problem occurred inside the multi-threaded loop, stop now !" << endl);
          return 1;
        }
      }
    }
  }
  // Normal end  
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::ImageUpdateStep()
{
  // Verbose
  if (m_verbose>=2 || m_optimizerImageStatFlag) Cout("iOptimizerADMMLim::ImageUpdateStep() -> Start image update" << endl);
  // Number of spaces for nice printing
  string spaces_tbf = ""; if (mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()>1) spaces_tbf += "  ";
  string spaces_rbf = ""; if (mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()>1) spaces_rbf += "  ";
  string spaces_cbf = ""; if (mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions()>1) spaces_cbf += "  ";
  // Loop over time basis functions
  for (int tbf=0; tbf<mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions(); tbf++)
  {
    // Verbose
    if (mp_ImageDimensionsAndQuantification->GetNbTimeBasisFunctions()>1 && (m_verbose>=3 || m_optimizerImageStatFlag))
    {
      if (mp_ImageDimensionsAndQuantification->GetTimeStaticFlag()) Cout(spaces_tbf << "--> Time frame " << tbf+1 << endl);
      else Cout(spaces_tbf << "--> Time basis function: " << tbf+1 << endl);
    }
    // Loop over respiratory basis functions
    for (int rbf=0; rbf<mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions(); rbf++)
    {
      // Verbose
      if (mp_ImageDimensionsAndQuantification->GetNbRespBasisFunctions()>1 && (m_verbose>=3 || m_optimizerImageStatFlag))
      {
        if (mp_ImageDimensionsAndQuantification->GetRespStaticFlag()) Cout(spaces_tbf << spaces_rbf << "--> Respiratory gate " << rbf+1 << endl);
        else Cout(spaces_tbf << spaces_rbf << "--> Respiratory basis function: " << rbf+1 << endl);
      }
      // Loop over cardiac basis functions
      for (int cbf=0; cbf<mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions(); cbf++)
      {
        // Verbose
        if (mp_ImageDimensionsAndQuantification->GetNbCardBasisFunctions()>1 && (m_verbose>=3 || m_optimizerImageStatFlag))
        {
          if (mp_ImageDimensionsAndQuantification->GetCardStaticFlag()) Cout(spaces_tbf << spaces_rbf << spaces_cbf << "--> Cardiac gate " << cbf+1 << endl);
          else Cout(spaces_tbf << spaces_rbf << spaces_cbf << "--> Cardiac basis function: " << cbf+1 << endl);
        }
        // Set the number of threads for image computation
        #ifdef CASTOR_OMP
        omp_set_num_threads(mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation());
        #endif
        // Reset counter for image stats
        if (m_optimizerImageStatFlag)
        {
          for (int th=0; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation(); th++)
          {
            mp_imageStatNbVox[th] = 0;
            mp_imageStatMin[th] = mp_ImageSpace->m4p_image[tbf][rbf][cbf][0];
            mp_imageStatMax[th] = mp_ImageSpace->m4p_image[tbf][rbf][cbf][0];
            mp_imageStatMean[th] = 0.;
            mp_imageStatVariance[th] = 0.;
            mp_correctionStatMean[th] = 0.;
            mp_correctionStatVariance[th] = 0.;
          }
        }
        // OMP loop over voxels
        INTNB v;
        bool error = false;
        #pragma omp parallel for private(v) schedule(guided)
        for (v=0 ; v<mp_ImageDimensionsAndQuantification->GetNbVoxXYZ() ; v++) 
        {
          // Get the thread index
          int th = 0;
          #ifdef CASTOR_OMP
          th = omp_get_thread_num();
          #endif
          // Compute sensitivity
          FLTNB sensitivity = ComputeSensitivity(mp_ImageSpace->m5p_sensitivity[0],tbf,rbf,cbf,v);
          // If the sensitivity is negative or null, we skip this voxel
          if (sensitivity<=0.) continue;
          // Keep the current image value in a buffer for update statistics
          FLTNB old_image_value = mp_ImageSpace->m4p_image[tbf][rbf][cbf][v];
          // Fill class-local buffer of backward values (thread-safe)
          for (int nb=0; nb<m_nbBackwardImages; nb++)
            m3p_backwardValues[th][0][nb] = mp_ImageSpace->m6p_backwardImage[nb][0][tbf][rbf][cbf][v];

          FLTNB a_currentImageValue = mp_ImageSpace->m4p_image[tbf][rbf][cbf][v];
          FLTNB* apImageValue = &mp_ImageSpace->m4p_image[tbf][rbf][cbf][v];
          //FLTNB a_sensitivity = sensitivity;
          FLTNB* ap_correctionValues = m3p_backwardValues[th][0];
          INTNB a_voxel = v;
          int a_tbf = tbf;
          int a_rbf = rbf;
          int a_cbf = cbf;
          int a_th = th;

          // Store gradient for next sub step
          int no_thread = 0;
          if (m_currentSubStep == 0 && !m_isInPreIteration)
          {
            // Scale penalty with respect to the number of subsets to get correct balance between likelihood and penalty
            HPFLTNB penalty = ((HPFLTNB)(m4p_firstDerivativePenaltyImage[a_tbf][a_rbf][a_cbf][a_voxel])) / ((HPFLTNB)(mp_nbSubsets[m_currentTotalSubIteration]));
            // Compute gradient
            m4p_currentGrad[a_tbf][a_rbf][a_cbf][a_voxel] = -m_alpha * *ap_correctionValues - penalty;
            mp_gradSquareNorm[a_th] += (HPFLTNB)m4p_currentGrad[a_tbf][a_rbf][a_cbf][a_voxel]*(HPFLTNB)m4p_currentGrad[a_tbf][a_rbf][a_cbf][a_voxel];
          }
          // Compute x update as gradient was updated
          if (m_currentSubStep == 1)
          {
            // Compute part of stepsize depending on penalty type
            HPFLTNB added_term = 0.;
            // If no penalty, then exit (the penalty image term has been initialized to 0)
            if (mp_Penalty==NULL)
            {
              added_term = 0.;
            }
            // If penalty is quadratic MRF
            else if (mp_Penalty->GetPenaltyID() == "MRF")
            {
              added_term = mp_penaltyGradSquareNorm[no_thread];
            }
            else if (mp_Penalty->GetPenaltyID() == "QUAD") // If penalty is quadratic
            {
              added_term = ((HPFLTNB)mp_Penalty->GetPenaltyStrength()*mp_gradSquareNorm[no_thread]);
            }
            // Compute stepsize avoiding division by 0
            HPFLTNB stepsize = (((HPFLTNB)m_alpha * mp_projGradSquareNorm[no_thread]) + added_term == 0.) ? 0. : mp_gradSquareNorm[no_thread] / (((HPFLTNB)m_alpha * mp_projGradSquareNorm[no_thread]) + added_term);
            //stepsize = 0.0001;
            // Compute additive image update factor
            HPFLTNB additive_image_update_factor = stepsize * (HPFLTNB)m4p_currentGrad[a_tbf][a_rbf][a_cbf][a_voxel];
            // Update image value
            *apImageValue = (HPFLTNB)a_currentImageValue + additive_image_update_factor;
          }

          // Do some image statistics
          if (m_optimizerImageStatFlag)
          {
            if (mp_imageStatMin[th]>mp_ImageSpace->m4p_image[tbf][rbf][cbf][v]) mp_imageStatMin[th] = mp_ImageSpace->m4p_image[tbf][rbf][cbf][v];
            if (mp_imageStatMax[th]<mp_ImageSpace->m4p_image[tbf][rbf][cbf][v]) mp_imageStatMax[th] = mp_ImageSpace->m4p_image[tbf][rbf][cbf][v];
            // The one pass algorithm is used to compute the variance here for both the image and correction step.
            // Do not change anyhting to the order of the operations !
            mp_imageStatNbVox[th]++;
            HPFLTNB sample = ((HPFLTNB)(mp_ImageSpace->m4p_image[tbf][rbf][cbf][v]));
            HPFLTNB delta = sample - mp_imageStatMean[th];
            mp_imageStatMean[th] += delta / ((HPFLTNB)(mp_imageStatNbVox[th]));
            mp_imageStatVariance[th] += delta * (sample - mp_imageStatMean[th]);
            // Same one-pass variance algorithm for the correction step
            sample = ((HPFLTNB)(mp_ImageSpace->m4p_image[tbf][rbf][cbf][v] - old_image_value));
            delta = sample - mp_correctionStatMean[th];
            mp_correctionStatMean[th] += delta / ((HPFLTNB)(mp_imageStatNbVox[th]));
            mp_correctionStatVariance[th] += delta * (sample - mp_correctionStatMean[th]);
          }
        }
        // Check for errors (we cannot stop the multi-threaded loop so we do it here)
        if (error) return 1;
        // Finish image update statistics computations: these operations are specific to the merging of the
        // multi-threaded one-pass variance estimations used above.
        if (m_optimizerImageStatFlag)
        {
          for (int th=1; th<mp_ImageDimensionsAndQuantification->GetNbThreadsForImageComputation(); th++)
          {
            // Image minimum and maximum
            if (mp_imageStatMin[0]>mp_imageStatMin[th]) mp_imageStatMin[0] = mp_imageStatMin[th];
            if (mp_imageStatMax[0]<mp_imageStatMax[th]) mp_imageStatMax[0] = mp_imageStatMax[th];
            // Cast the counts into HPFLTNB
            HPFLTNB count1 = ((HPFLTNB)(mp_imageStatNbVox[0]));
            HPFLTNB count2 = ((HPFLTNB)(mp_imageStatNbVox[th]));
            // Update the count
            mp_imageStatNbVox[0] += mp_imageStatNbVox[th];
            HPFLTNB count12 = ((HPFLTNB)(mp_imageStatNbVox[0]));
            // Compute the delta mean
            HPFLTNB delta = mp_imageStatMean[0] - mp_imageStatMean[th];
            // Update the variance
            mp_imageStatVariance[0] += mp_imageStatVariance[th] + delta * delta * count1 * count2 / count12;
            // Update the mean
            mp_imageStatMean[0] = (count1*mp_imageStatMean[0] + count2*mp_imageStatMean[th]) / count12;
            // Do the same for the correction step
            delta = mp_correctionStatMean[0] - mp_correctionStatMean[th];
            mp_correctionStatVariance[0] += mp_correctionStatVariance[th] + delta * delta * count1 * count2 / count12;
            mp_correctionStatMean[0] = (count1*mp_correctionStatMean[0] + count2*mp_correctionStatMean[th]) / count12;
          }
          // Finish variance computation
          mp_imageStatVariance[0] /= ((HPFLTNB)(mp_imageStatNbVox[0]));
          mp_correctionStatVariance[0] /= ((HPFLTNB)(mp_imageStatNbVox[0]));
          // Print out
          Cout(spaces_tbf << spaces_rbf << spaces_cbf << "  --> Image update step "
                                                      << " | mean: " << mp_correctionStatMean[0]
                                                      << " | stdv: " << sqrt(mp_correctionStatVariance[0]) << endl);
          Cout(spaces_tbf << spaces_rbf << spaces_cbf << "  --> New image estimate"
                                                      << " | mean: " << mp_imageStatMean[0]
                                                      << " | stdv: " << sqrt(mp_imageStatVariance[0])
                                                      << " | min/max: " << mp_imageStatMin[0] << "/" << mp_imageStatMax[0] << endl);
        }
      }
    }
  }

  // End
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================

int iOptimizerADMMLim::ImageSpaceSpecificOperations( FLTNB a_currentImageValue, FLTNB* apImageValue,
                                                  FLTNB a_sensitivity, FLTNB* ap_correctionValues,
                                                  INTNB a_voxel, int a_tbf, int a_rbf, int a_cbf )
{
  // Do not call pure virtual function for specific data space operations ImageSpaceSpecificOperations because we need current thread for gradient norm computation
  return 0;
}

// =====================================================================
// ---------------------------------------------------------------------
// ---------------------------------------------------------------------
// =====================================================================