new changes

Author: nawazishkhan1-nk

Libs/Optimize/Optimize.cpp

@@ -359,6 +359,7 @@ void Optimize::SetVerbosity(int verbosity_level)
 void Optimize::SetDomainsPerShape(int domains_per_shape)
 {
   this->m_domains_per_shape = domains_per_shape;
+  std::cout << "Initially, dps set to = "<< domains_per_shape << std::endl;
   this->m_sampler->SetDomainsPerShape(this->m_domains_per_shape);
 }
 
@@ -586,9 +587,14 @@ void Optimize::InitializeSampler()
   m_sampler->GetEnsembleEntropyFunction()->SetRecomputeCovarianceInterval(1);
   m_sampler->GetEnsembleEntropyFunction()->SetHoldMinimumVariance(false);
 
-  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVariance(m_starting_regularization);
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVariance(m_starting_regularization); // For Initialization set start_reg_params
   m_sampler->GetEnsembleMlpcaEntropyFunction()->SetRecomputeCovarianceInterval(1);
   m_sampler->GetEnsembleMlpcaEntropyFunction()->SetHoldMinimumVariance(false);
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetDomainsPerShapeInfo(m_domains_per_shape);
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceWithin(m_starting_regularization_multilevel); // For Initialization set start_reg_params
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetHoldMinimumVarianceWithin(false);
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceBetween(m_starting_regularization_multilevel_between); // For Initialization set start_reg_params
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetHoldMinimumVarianceBetween(false);
 
   std::cout << "mlpca initialized" << std::endl;
   m_sampler->GetMeshBasedGeneralEntropyGradientFunction()->SetMinimumVariance(
@@ -621,15 +627,15 @@ void Optimize::InitializeSampler()
     ->SetRecomputeCovarianceInterval(m_recompute_regularization_interval);
   m_sampler->GetEnsembleMlpcaEntropyFunction()
     ->SetRecomputeCovarianceInterval(m_recompute_regularization_interval);
-  std::cout << "mlpca initialized 2" << std::endl;
+  // std::cout << "mlpca initialized 2" << std::endl;
 
 
   // These flags must be set before Initialize, since Initialize need to know which domains are fixed ahead of time
   for (unsigned int i = 0; i < this->m_domain_flags.size(); i++) {
     this->GetSampler()->GetParticleSystem()->FlagDomain(this->m_domain_flags[i]);
   }
   m_sampler->Initialize();
-  std::cout << "Sampler Init 1 done" << std::endl;
+  // std::cout << "Sampler Init 1 done" << std::endl;
 
   m_sampler->GetOptimizer()->SetTolerance(0.0);
 
@@ -728,18 +734,24 @@ void Optimize::Initialize()
     m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceDecay(m_starting_regularization,
                                                                      m_ending_regularization,
                                                                      m_iterations_per_split);
-
+    m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceDecayWithin(m_starting_regularization_multilevel,
+                                                                     m_ending_regularization_multilevel,
+                                                                     m_iterations_per_split);
+    m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceDecayBetween(m_starting_regularization_multilevel_between,
+                                                                     m_ending_regularization_multilevel_between,
+                                                                     m_iterations_per_split);
+                                                                     
     m_sampler->GetMeshBasedGeneralEntropyGradientFunction()->SetMinimumVarianceDecay(
       m_starting_regularization,
       m_ending_regularization,
       m_iterations_per_split);
   }
   else {
     // force to mean
-    if (m_use_mlpca_optimize){
-      m_sampler->SetCorrespondenceMode(shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropyMeanEnergy);
-    }
-    else if ((m_attributes_per_domain.size() > 0 &&
+    // if (m_use_mlpca_optimize){
+    //   m_sampler->SetCorrespondenceMode(shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropyMeanEnergy);
+    // }
+    if ((m_attributes_per_domain.size() > 0 &&
          *std::max_element(m_attributes_per_domain.begin(),
                            m_attributes_per_domain.end()) > 0) || m_mesh_based_attributes) {
       m_sampler->SetCorrespondenceMode(shapeworks::CorrespondenceMode::MeshBasedGeneralMeanEnergy);
@@ -1053,6 +1065,16 @@ void Optimize::RunOptimize()
     m_ending_regularization,
     m_optimization_iterations -
     m_optimization_iterations_completed);
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceDecayWithin(
+    m_starting_regularization_multilevel,
+    m_ending_regularization_multilevel,
+    m_optimization_iterations -
+    m_optimization_iterations_completed);
+  m_sampler->GetEnsembleMlpcaEntropyFunction()->SetMinimumVarianceDecayBetween(
+    m_starting_regularization_multilevel_between,
+    m_ending_regularization_multilevel_between,
+    m_optimization_iterations -
+    m_optimization_iterations_completed);
 
 
   m_sampler->SetCorrespondenceOn(); //  B - ON
@@ -2106,6 +2128,23 @@ void Optimize::SetStartingRegularization(double starting_regularization)
 void Optimize::SetEndingRegularization(double ending_regularization)
 { this->m_ending_regularization = ending_regularization; }
 
+//---------------------------------------------------------------------------
+void Optimize::SetStartingRegularizationMultilevelWithin(std::vector<double> reg_params_start)
+{ this->m_starting_regularization_multilevel = reg_params_start; }
+
+//---------------------------------------------------------------------------
+void Optimize::SetEndingRegularizationMultilevelWithin(std::vector<double> reg_params_end)
+{ this->m_ending_regularization_multilevel = reg_params_end; }
+
+//---------------------------------------------------------------------------
+void Optimize::SetStartingRegularizationMultilevelBetween(double starting_regularization)
+{ this->m_starting_regularization_multilevel_between = starting_regularization; }
+
+//---------------------------------------------------------------------------
+void Optimize::SetEndingRegularizationMultilevelBetween(double ending_regularization)
+{ this->m_ending_regularization_multilevel_between = ending_regularization; }
+
+
 //---------------------------------------------------------------------------
 void Optimize::SetRecomputeRegularizationInterval(int recompute_regularization_interval)
 { this->m_recompute_regularization_interval = recompute_regularization_interval; }

Libs/Optimize/Optimize.h

@@ -189,6 +189,15 @@ class Optimize {
   void SetStartingRegularization(double starting_regularization);
   //! Set the ending regularization (TODO: details)
   void SetEndingRegularization(double ending_regularization);
+
+  void SetStartingRegularizationMultilevelWithin(std::vector<double> reg_params_start);
+  //! Set the ending regularization (TODO: details)
+  void SetEndingRegularizationMultilevelWithin(std::vector<double> reg_params_end);
+  void SetStartingRegularizationMultilevelBetween(double reg_params_start);
+  //! Set the ending regularization (TODO: details)
+  void SetEndingRegularizationMultilevelBetween(double reg_params_end);
+
+
   //! Set the interval for recomputing regularization (TODO: details)
   void SetRecomputeRegularizationInterval(int recompute_regularization_interval);
   //! Set if initialization splits should be saved or not
@@ -395,6 +404,10 @@ class Optimize {
   double m_initial_relative_weighting = 0.05;
   double m_starting_regularization = 1000;
   double m_ending_regularization = 1.0;
+  std::vector<double> m_ending_regularization_multilevel;
+  std::vector<double> m_starting_regularization_multilevel;
+  double m_starting_regularization_multilevel_between = 1000;
+  double m_ending_regularization_multilevel_between = 1.0;
   int m_recompute_regularization_interval = 1;
   bool m_save_init_splits = false;
   unsigned int m_checkpointing_interval = 50;

Libs/Optimize/OptimizeParameters.cpp

@@ -149,6 +149,54 @@ bool OptimizeParameters::get_use_mlpca_optimize()
   return use_mlpca_optimize;
 }
 
+//---------------------------------------------------------------------------
+std::vector<double> OptimizeParameters::get_starting_regularization_multilevel()
+{
+  return this->params_.get("starting_regularization_multilevel", {1000.0});
+}
+
+//---------------------------------------------------------------------------
+void OptimizeParameters::set_starting_regularization_multilevel(std::vector<double> reg_params_start)
+{
+  return this->params_.set("starting_regularization_multilevel", reg_params_start);
+}
+
+//---------------------------------------------------------------------------
+std::vector<double> OptimizeParameters::get_ending_regularization_multilevel()
+{
+  return this->params_.get("ending_regularization_multilevel", {10.0});
+}
+
+//---------------------------------------------------------------------------
+void OptimizeParameters::set_ending_regularization_multilevel(std::vector<double> reg_params_end)
+{
+  return this->params_.set("ending_regularization_multilevel", reg_params_end);
+}
+
+//---------------------------------------------------------------------------
+double OptimizeParameters::get_starting_regularization_between()
+{
+  return this->params_.get("starting_regularization_between", 1000.0);
+}
+
+//---------------------------------------------------------------------------
+void OptimizeParameters::set_starting_regularization_between(double value)
+{
+  this->params_.set("starting_regularization_between", value);
+}
+
+//---------------------------------------------------------------------------
+double OptimizeParameters::get_ending_regularization_between()
+{
+  return this->params_.get("ending_regularization_between", 10.0);
+}
+
+//---------------------------------------------------------------------------
+void OptimizeParameters::set_ending_regularization_between(double value)
+{
+  this->params_.set("ending_regularization_between", value);
+}
+
 //---------------------------------------------------------------------------
 double OptimizeParameters::get_normals_strength()
 {
@@ -265,6 +313,16 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize)
   }
 
   optimize->SetMlpcaOptimize(use_mlpca_optimize);
+  if(use_mlpca_optimize)
+  {
+    //Set reg parameters for each organ indiviually
+
+    optimize->SetStartingRegularizationMultilevelWithin(this->get_starting_regularization_multilevel());
+    optimize->SetEndingRegularizationMultilevelWithin(this->get_ending_regularization_multilevel());
+    optimize->SetStartingRegularizationMultilevelBetween(this->get_starting_regularization_between());
+    optimize->SetEndingRegularizationMultilevelBetween(this->get_ending_regularization_between());
+
+  }
   optimize->SetUseNormals(use_normals);
   optimize->SetUseXYZ(use_xyz);
   optimize->SetUseMeshBasedAttributes(normals_enabled);

Libs/Optimize/OptimizeParameters.h

@@ -51,6 +51,18 @@ class OptimizeParameters {
   bool get_use_mlpca_optimize();
   void set_use_mlpca_optimize(bool value);
 
+  std::vector<double> get_starting_regularization_multilevel();
+  void set_starting_regularization_multilevel(std::vector<double> reg_params_start);
+
+  std::vector<double> get_ending_regularization_multilevel();
+  void set_ending_regularization_multilevel(std::vector<double> reg_params_end);
+
+  double get_starting_regularization_between();
+  void set_starting_regularization_between(double value);
+
+  double get_ending_regularization_between();
+  void set_ending_regularization_between(double value);
+
   double get_normals_strength();
   void set_normals_strength(double value);
 

Libs/Optimize/ParticleSystem/CorrespondenceMode.h

@@ -9,7 +9,7 @@ namespace shapeworks {
     EnsembleMixedEffectsEntropy = 4,
     MeshBasedGeneralEntropy = 5,
     MeshBasedGeneralMeanEnergy = 6,
-    MlpcaBasedEnsembleEntropy = 7,
-    MlpcaBasedEnsembleEntropyMeanEnergy = 8
+    MlpcaBasedEnsembleEntropy = 7
+    // MlpcaBasedEnsembleEntropyMeanEnergy = 8
   };
 }

Libs/Optimize/ParticleSystem/Sampler.cpp

@@ -36,7 +36,6 @@ Sampler::Sampler()
   m_EnsembleRegressionEntropyFunction = itk::ParticleEnsembleEntropyFunction<Dimension>::New();
   m_EnsembleMixedEffectsEntropyFunction = itk::ParticleEnsembleEntropyFunction<Dimension>::New();
   m_MeshBasedGeneralEntropyGradientFunction = itk::ParticleMeshBasedGeneralEntropyGradientFunction<Dimension>::New();
-  // TODO: Initialize Mlpca optimization function
   m_MlpcaBasedEnsembleEntropyFunction = itk::ParticleEnsembleMlpcaEntropyFunction<Dimension>::New();
   std::cout << "sampler constructor " << std::endl;
   m_ShapeMatrix = itk::ParticleShapeMatrixAttribute<double, Dimension>::New();
@@ -269,6 +268,7 @@ void Sampler::Execute()
 
   //this->GetOptimizer()->SetShapeMatrix(this->m_ShapeMatrix);
   if(m_CorrespondenceMode == shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropy){
+  // if(m_CorrespondenceMode == shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropy || m_CorrespondenceMode == shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropyMeanEnergy){
     std::cout << "Inside MLPCA sampler execute" << std::endl;
     this->GetOptimizer()->StartMlpcaOptimization(m_LinkingFunction->GetBOn());
   }

Libs/Optimize/ParticleSystem/Sampler.h

@@ -272,12 +272,12 @@ class Sampler {
       m_LinkingFunction->SetFunctionB(m_MlpcaBasedEnsembleEntropyFunction);
       m_MlpcaBasedEnsembleEntropyFunction->UseEntropy();
     }
-    else if (mode == shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropyMeanEnergy) {
-      std::cout << "mode IS --> H " << std::endl;
-      std::cout << "Linking set to mlpca ensemble Mean Energy entropy BEFORE" << std::endl;
-      m_LinkingFunction->SetFunctionB(m_MlpcaBasedEnsembleEntropyFunction);
-      m_MlpcaBasedEnsembleEntropyFunction->UseMeanEnergy();
-    }
+    // else if (mode == shapeworks::CorrespondenceMode::MlpcaBasedEnsembleEntropyMeanEnergy) {
+    //   std::cout << "mode IS --> H " << std::endl;
+    //   std::cout << "Linking set to mlpca ensemble Mean Energy entropy BEFORE" << std::endl;
+    //   m_LinkingFunction->SetFunctionB(m_MlpcaBasedEnsembleEntropyFunction);
+    //   m_MlpcaBasedEnsembleEntropyFunction->UseMeanEnergy();
+    // }
 
     m_CorrespondenceMode = mode;
 

Libs/Optimize/ParticleSystem/itkParticleDualVectorFunction.h

@@ -155,6 +155,32 @@ class ParticleDualVectorFunction : public ParticleVectorFunction<VDimension>
         return ansB;
     }
 
+    virtual double EnergyBBetween(unsigned int idx, unsigned int d, const ParticleSystemType *system) const
+    {
+        m_FunctionB->BeforeEvaluate(idx, d, system);
+        double ansB = 0.0;
+        if (m_BOn == true)
+        {
+            const itk::ParticleEnsembleMlpcaEntropyFunction<VDimension>* mlpca_function_b = dynamic_cast <const itk::ParticleEnsembleMlpcaEntropyFunction<VDimension>*> (m_FunctionB.GetPointer());
+            ansB = mlpca_function_b->EnergyFromWithinResiduals(idx, d, system);
+        }
+        ansB *= m_RelativeEnergyScaling;
+        return ansB;
+    }
+
+    virtual double EnergyBWithin(unsigned int idx, unsigned int d, const ParticleSystemType *system) const
+    {
+        m_FunctionB->BeforeEvaluate(idx, d, system);
+        double ansB = 0.0;
+        if (m_BOn == true)
+        {   
+            const itk::ParticleEnsembleMlpcaEntropyFunction<VDimension>* mlpca_function_b = dynamic_cast <const itk::ParticleEnsembleMlpcaEntropyFunction<VDimension>*> (m_FunctionB.GetPointer());
+            ansB = mlpca_function_b->EnergyFromBetweenResiduals(idx, d, system);
+        }
+        ansB *= m_RelativeEnergyScaling; // TODO: Ensure proper usage here 
+        return ansB;
+    }
+
     virtual double Energy(unsigned int idx, unsigned int d, const ParticleSystemType *system) const
     {
         double ansA = 0.0;
@@ -275,7 +301,7 @@ class ParticleDualVectorFunction : public ParticleVectorFunction<VDimension>
         {
             // B is active, A is not active
             finalEnergy = ansB;
-            return finalEnergy ;
+            return m_RelativeEnergyScaling * finalEnergy ;
 
         }
         return 0.0;
@@ -387,19 +413,25 @@ class ParticleDualVectorFunction : public ParticleVectorFunction<VDimension>
 
         // evaluate individual functions: A = surface energy, B = correspondence
         if (m_AOn == true)
-        {
+        {   
+            std::cout << "evaluating A" << std::endl;
             ansA = m_FunctionA->Evaluate(idx, d, system, maxA, energyA);
             const_cast<ParticleDualVectorFunction *>(this)->m_AverageGradMagA = m_AverageGradMagA + ansA.magnitude();
             const_cast<ParticleDualVectorFunction *>(this)->m_AverageEnergyA = m_AverageEnergyA + energyA;
+            std::cout << "evaluating A done" << std::endl;
+
         }
 
         if (m_BOn == true)
         {   
+            std::cout << "evaluating B" << std::endl;
             const itk::ParticleEnsembleMlpcaEntropyFunction<VDimension>* mlpca_function_b = dynamic_cast <const itk::ParticleEnsembleMlpcaEntropyFunction<VDimension>*> (m_FunctionB.GetPointer());
             ansB = mlpca_function_b->EvaluateFromBetweenResiduals(idx, d, system, maxB, energyB);
             const_cast<ParticleDualVectorFunction *>(this)->m_AverageWithinGradMagB = m_AverageWithinGradMagB + ansB.magnitude();
             const_cast<ParticleDualVectorFunction *>(this)->m_AverageWithinEnergyB = m_AverageWithinEnergyB + energyB;
-            std::cout << "WITHIN Subspace Grad Mag " << ansB.magnitude() << " Energy " << energyB << std::endl;
+            // std::cout << "WITHIN Subspace Grad Mag " << ansB.magnitude() << " Energy " << energyB << std::endl;
+            std::cout << "evaluating B done" << std::endl;
+
         }
 
         if( m_RelativeEnergyScaling == 0.0)
@@ -506,7 +538,7 @@ class ParticleDualVectorFunction : public ParticleVectorFunction<VDimension>
             // std::cout << " inside B between eval done" << std::endl;
             const_cast<ParticleDualVectorFunction *>(this)->m_AverageBetweenGradMagB = m_AverageBetweenGradMagB + ansB.magnitude();
             const_cast<ParticleDualVectorFunction *>(this)->m_AverageBetweenEnergyB = m_AverageBetweenEnergyB + energyB;
-            std::cout << "BETWEEN Subspace Grad Mag " << ansB.magnitude() << " Energy " << energyB << std::endl;
+            // std::cout << "BETWEEN Subspace Grad Mag " << ansB.magnitude() << " Energy " << energyB << std::endl;
         }
 
         if( m_RelativeEnergyScaling == 0.0)