Added per base scenario weights

src/main.cpp

@@ -4,6 +4,7 @@
 #include "reducedmodeloptimizer.hpp"
 #include "simulationhistoryplotter.hpp"
 #include "trianglepattterntopology.hpp"
+#include <boost/algorithm/string.hpp>
 #include <chrono>
 #include <filesystem>
 #include <iostream>
@@ -12,7 +13,6 @@
 #include <string>
 #include <string_view>
 #include <vcg/complex/algorithms/update/position.h>
-#include <boost/algorithm/string.hpp>
 
 #ifdef POLYSCOPE_DEFINED
 #include "polyscope/curve_network.h"
@@ -21,7 +21,7 @@
 #endif
 int main(int argc, char *argv[])
 {
-    if (argc <= 7) {
+    if (argc <= 6) {
         std::cerr << "Wrong number of input parameters. Expects at least 6 input parameters."
                      "Usage:\n"
                      "1)full pattern file path\n"
@@ -29,8 +29,8 @@ int main(int argc, char *argv[])
                      "3)Number of optimizaion function calls\n"
                      "4)Translational error weight\n"
                      "5)Optimization results directory path\n"
-                     "6)[optional]Optimization parameters\n"
+                     "6)[optional]Intermediate results directory path\n"
-                     "7)[optional]Intermediate results directory path\n"
+                     "7)[optional]Optimization parameters\n"
                      "Exiting.."
                   << std::endl;
         return 1;
@@ -49,6 +49,12 @@ int main(int argc, char *argv[])
     reducedPattern.scale(0.03, interfaceNodeIndex);
 
     // Set the optization settings
+    const std::filesystem::path optimizationSettingsFilePath = argv[3];
+    if (!std::filesystem::exists(optimizationSettingsFilePath)) {
+        std::cerr << "Input optimization settings file does not exist:"
+                  << optimizationSettingsFilePath << std::endl;
+    }
+
     ReducedPatternOptimization::xRange beamE{"E", 0.001, 1000};
     ReducedPatternOptimization::xRange beamA{"A", 0.001, 1000};
     ReducedPatternOptimization::xRange beamI2{"I2", 0.001, 1000};
@@ -57,65 +63,77 @@ int main(int argc, char *argv[])
     ReducedPatternOptimization::xRange innerHexagonSize{"R", 0.05, 0.95};
     ReducedPatternOptimization::xRange innerHexagonAngle{"Theta", -30.0, 30.0};
     ReducedPatternOptimization::Settings settings_optimization;
-    settings_optimization.parameterRanges
+    settings_optimization.variablesRanges
         = {beamE, beamA, beamI2, beamI3, beamJ, innerHexagonSize, innerHexagonAngle};
     settings_optimization.numberOfFunctionCalls = std::atoi(argv[3]);
 
     const bool input_optimizationParametersDefined = argc >= 8;
     if (input_optimizationParametersDefined) {
         const std::string optimizationParametersTag = argv[7];
-        std::vector <std::string> split=Utilities::split(optimizationParametersTag,",");
+        std::vector<std::string> split = Utilities::split(optimizationParametersTag, ",");
         //parse parameter tag
         std::vector<std::vector<ReducedPatternOptimization::OptimizationParameterIndex>>
             optimizationParameters;
         std::vector<ReducedPatternOptimization::OptimizationParameterIndex> parameterGroup;
         for (std::string &s : split) {
-//        std::cout<<s<<std::endl;
+            //        std::cout<<s<<std::endl;
-//        s=Utilities::trimLeftAndRightSpaces(s);
+            //        s=Utilities::trimLeftAndRightSpaces(s);
-            if(boost::algorithm::contains(s,"{")){
+            if (boost::algorithm::contains(s, "{")) {
-//            std::cout<<"{ detected"<<std::endl;
+                //            std::cout<<"{ detected"<<std::endl;
                 parameterGroup.clear();
-                }
+            }
-               if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::E)))){
+            if (boost::algorithm::contains(s,
+                                           std::to_string(
+                                               static_cast<int>(ReducedPatternOptimization::E)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::E);
-                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::A)))){
+            } else if (boost::algorithm::contains(s,
+                                                  std::to_string(static_cast<int>(
+                                                      ReducedPatternOptimization::A)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::A);
-                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::I2)))){
+            } else if (boost::algorithm::contains(s,
+                                                  std::to_string(static_cast<int>(
+                                                      ReducedPatternOptimization::I2)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::I2);
-                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::I3)))){
+            } else if (boost::algorithm::contains(s,
+                                                  std::to_string(static_cast<int>(
+                                                      ReducedPatternOptimization::I3)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::I3);
-                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::J)))){
+            } else if (boost::algorithm::contains(s,
+                                                  std::to_string(static_cast<int>(
+                                                      ReducedPatternOptimization::J)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::J);
-                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::R)))){
+            } else if (boost::algorithm::contains(s,
+                                                  std::to_string(static_cast<int>(
+                                                      ReducedPatternOptimization::R)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::R);
-                }else if(boost::algorithm::contains(s,std::to_string(static_cast<int>(ReducedPatternOptimization::Theta)))){
+            } else if (boost::algorithm::contains(s,
+                                                  std::to_string(static_cast<int>(
+                                                      ReducedPatternOptimization::Theta)))) {
                 parameterGroup.push_back(ReducedPatternOptimization::Theta);
-                }
+            } else {
-                else{
                 std::cerr << "Wrong optimization parameter input: " << optimizationParametersTag
                           << std::endl;
-
+            }
-                }
+            if (boost::algorithm::contains(s, "}")) {
-                if(boost::algorithm::contains(s,"}")){
                 optimizationParameters.push_back(parameterGroup);
-                }
+            }
-
         }
-        settings_optimization.optimizationVariables=optimizationParameters;
+        settings_optimization.optimizationStrategy = optimizationParameters;
-        #ifdef POLYSCOPE_DEFINED
+#ifdef POLYSCOPE_DEFINED
-        for(const std::vector<ReducedPatternOptimization::OptimizationParameterIndex>& v:settings_optimization.optimizationVariables){
+        for (const std::vector<ReducedPatternOptimization::OptimizationParameterIndex> &v :
-        for(const ReducedPatternOptimization::OptimizationParameterIndex& i:v ){
+             settings_optimization.optimizationStrategy) {
-        std::cout<<static_cast<int>(i)<<" ";
+            for (const ReducedPatternOptimization::OptimizationParameterIndex &i : v) {
+                std::cout << static_cast<int>(i) << " ";
+            }
+            std::cout << std::endl;
         }
-        std::cout<<std::endl;
+#endif
-        }
-        #endif
 
     } else {
         enum OptimizationParameterComparisonScenarioIndex {
             AllVar,
             GeoYM,
-            noYM,
+            MatGeo,
             YMMat_Geo,
             YM_MatGeo,
             MatGeo_YM,
@@ -134,7 +152,7 @@ int main(int argc, char *argv[])
                 using namespace ReducedPatternOptimization;
                 optimizationParameters[AllVar] = {{E, A, I2, I3, J, R, Theta}};
                 optimizationParameters[GeoYM] = {{R, Theta, E}};
-                optimizationParameters[noYM] = {{A, I2, I3, J, R, Theta}};
+                optimizationParameters[MatGeo] = {{A, I2, I3, J, R, Theta}};
                 optimizationParameters[YMMat_Geo] = {{E, A, I2, I3, J}, {R, Theta}};
                 optimizationParameters[YM_MatGeo] = {{E}, {A, I2, I3, J, R, Theta}};
                 optimizationParameters[MatGeo_YM] = {{A, I2, I3, J, R, Theta}, {E}};
@@ -145,111 +163,158 @@ int main(int argc, char *argv[])
                 return optimizationParameters;
             }();
 
-        constexpr OptimizationParameterComparisonScenarioIndex scenario = YMGeo_Mat;
+        constexpr OptimizationParameterComparisonScenarioIndex scenario = MatGeo;
-        settings_optimization.optimizationVariables = optimizationParameters[scenario];
+        settings_optimization.optimizationStrategy = optimizationParameters[scenario];
+        if (scenario == YMGeo_Mat) {
+            settings_optimization.optimizationVariablesGroupsWeights = {0.15, 0.85};
+        }
     }
-    if(settings_optimization.optimizationVariables.empty()){
+    if (settings_optimization.optimizationStrategy.empty()) {
-    std::cerr<<"No optimization variables. Exiting.."<<std::endl;
+        std::cerr << "No optimization variables. Exiting.." << std::endl;
-    std::terminate();
+        std::terminate();
-    }
-
-    const bool input_intermediateResultsDirectoryDefined = argc >= 7;
-    if (input_intermediateResultsDirectoryDefined) {
-        settings_optimization.intermediateResultsDirectoryPath = std::filesystem::path(argv[6]);
     }
 
     settings_optimization.normalizationStrategy
         = ReducedPatternOptimization::Settings::NormalizationStrategy::Epsilon;
-#ifdef POLYSCOPE_DEFINED
+    //#ifdef POLYSCOPE_DEFINED
-  settings_optimization.translationNormalizationParameter = 0;
+    //#else
-  settings_optimization.rotationNormalizationParameter = 0;
+    settings_optimization.translationNormalizationParameter = 4e-3;
-#else
+    settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
-  settings_optimization.translationNormalizationParameter = 3e-4;
-  settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
-#endif
-  settings_optimization.solverAccuracy = 1e-1;
-  settings_optimization.objectiveWeights.translational = std::atof(argv[4]);
-  settings_optimization.objectiveWeights.rotational = 2 - std::atof(argv[4]);
 
-  std::string xConcatNames;
+    //    settings_optimization.translationNormalizationParameter = 0;
-  for (const auto &x : settings_optimization.parameterRanges) {
+    //    settings_optimization.rotationNormalizationParameter = 0;
-      xConcatNames.append(x.label + "_");
-  }
-  xConcatNames.pop_back();
 
-  // Optimize pairthere
+    //    settings_optimization.translationNormalizationParameter = 1e-3;
-  const std::string pairName = fullPattern.getLabel(); // + "@" + reducedPattern.getLabel();
+    //    settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
-  const std::string optimizationName = pairName + "("
-                                       + std::to_string(settings_optimization.numberOfFunctionCalls)
-                                       + "_"
-                                       + to_string_with_precision(
-                                           settings_optimization.objectiveWeights.translational)
-                                       + ")" + "_" + xConcatNames;
-  const std::string optimizationResultsDirectory = argv[5];
-  std::string resultsOutputDir;
-  bool optimizationResultFolderExists = false;
-  const std::filesystem::path crashedJobsDirPath(std::filesystem::path(optimizationResultsDirectory)
-                                                     .append("CrashedJobs")
-                                                     .append(optimizationName));
-  if (std::filesystem::exists(crashedJobsDirPath)) {
-      resultsOutputDir = crashedJobsDirPath.string();
-      optimizationResultFolderExists = true;
-  }
-  const std::filesystem::path convergedJobsDirPath(
-      std::filesystem::path(optimizationResultsDirectory)
-          .append("ConvergedJobs")
-          .append(optimizationName));
-  if (std::filesystem::exists(convergedJobsDirPath)) {
-      resultsOutputDir = convergedJobsDirPath.string();
-      optimizationResultFolderExists = true;
-  }
 
-  ReducedPatternOptimization::Results optimizationResults;
+    //    settings_optimization.translationNormalizationParameter = 5e-2;
-  constexpr bool shouldReoptimize = true;
+    //    settings_optimization.rotationNormalizationParameter = vcg::math::ToRad(3.0);
-  bool optimizationAlreadyComputed = false;
-  if (!shouldReoptimize && optimizationResultFolderExists) {
-      const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir);
-      if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) {
-          optimizationAlreadyComputed = true;
-      }
-  }
 
-  if (!optimizationAlreadyComputed) {
+    //#endif
-      auto start = std::chrono::system_clock::now();
+    settings_optimization.solverAccuracy = 1e-2;
-      const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
+    //TODO refactor that
-      assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]);
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial]
-      ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
+        .translational
-      optimizer.initializePatterns(fullPattern, reducedPattern, settings_optimization.parameterRanges);
+        = 1.5;
-      optimizer.optimize(settings_optimization, optimizationResults);
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial].rotational
-      optimizationResults.label = optimizationName;
+        = 2
-      optimizationResults.baseTriangleFullPattern.copy(fullPattern);
+          - settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial]
-      optimizationResults.settings = settings_optimization;
+                .translational;
-      auto end = std::chrono::system_clock::now();
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear]
-      auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
+        .translational
-      optimizationResults.time = elapsed.count() / 1000.0;
+        = 1.5;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear].rotational
+        = 2
+          - settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Shear]
+                .translational;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending]
+        .translational
+        = 1.2;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending]
+        .rotational
+        = 2
+          - settings_optimization
+                .perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Bending]
+                .translational;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome]
+        .translational
+        = 1.95;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome].rotational
+        = 2
+          - settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Dome]
+                .translational;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle]
+        .translational
+        = 1.2;
+    settings_optimization.perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle]
+        .rotational
+        = 2
+          - settings_optimization
+                .perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Saddle]
+                .translational;
 
-      if (!optimizationResults.wasSuccessful) {
+    // Optimize pairthere
-          //          resultsOutputDir = crashedJobsDirPath.string();
+    const std::string pairName = fullPattern.getLabel(); // + "@" + reducedPattern.getLabel();
-          return 1;
+    const std::string optimizationName
-      }
+        = pairName + "(" + std::to_string(settings_optimization.numberOfFunctionCalls) + "_"
-      resultsOutputDir = convergedJobsDirPath.string();
+          + to_string_with_precision(
-      optimizationResults.save(resultsOutputDir, true);
+              settings_optimization
-      csvFile csv_resultsLocalFile(std::filesystem::path(resultsOutputDir).append("results.csv"),
+                  .perBaseScenarioObjectiveWeights[ReducedPatternOptimization::Axial]
-                                   true);
+                  .translational)
-      csvFile csv_results({}, false);
+          + ")";
-      std::vector<csvFile *> csvVector{&csv_resultsLocalFile, &csv_results};
+    const std::string optimizationResultsDirectory = argv[5];
-      csv_results << "Name";
+    std::string resultsOutputDir;
-      csv_resultsLocalFile << "Name";
+    bool optimizationResultFolderExists = false;
-      optimizationResults.writeHeaderTo(csvVector);
+    const std::filesystem::path crashedJobsDirPath(
-      settings_optimization.writeHeaderTo(csv_results);
+        std::filesystem::path(optimizationResultsDirectory)
-      csv_results << endrow;
+            .append("CrashedJobs")
-      csv_resultsLocalFile << endrow;
+            .append(optimizationName));
-      csv_results << std::to_string(fullPattern.EN()) + "#" + pairName;
+    if (std::filesystem::exists(crashedJobsDirPath)) {
-      optimizationResults.writeResultsTo(csvVector);
+        resultsOutputDir = crashedJobsDirPath.string();
-      settings_optimization.writeSettingsTo(csv_results);
+        optimizationResultFolderExists = true;
-      csv_results << endrow;
+    }
-      csv_resultsLocalFile << endrow;
+    const std::filesystem::path convergedJobsDirPath(
-  }
+        std::filesystem::path(optimizationResultsDirectory)
+            .append("ConvergedJobs")
+            .append(optimizationName));
+    if (std::filesystem::exists(convergedJobsDirPath)) {
+        resultsOutputDir = convergedJobsDirPath.string();
+        optimizationResultFolderExists = true;
+    }
+
+    ReducedPatternOptimization::Results optimizationResults;
+    constexpr bool shouldReoptimize = true;
+    bool optimizationAlreadyComputed = false;
+    if (!shouldReoptimize && optimizationResultFolderExists) {
+        const bool resultsWereSuccessfullyLoaded = optimizationResults.load(resultsOutputDir);
+        if (resultsWereSuccessfullyLoaded && optimizationResults.settings == settings_optimization) {
+            optimizationAlreadyComputed = true;
+        }
+    }
+
+    if (!optimizationAlreadyComputed) {
+        auto start = std::chrono::system_clock::now();
+        const std::vector<size_t> numberOfNodesPerSlot{1, 0, 0, 2, 1, 2, 1};
+        assert(interfaceNodeIndex == numberOfNodesPerSlot[0] + numberOfNodesPerSlot[3]);
+        ReducedModelOptimizer optimizer(numberOfNodesPerSlot);
+        const bool input_intermediateResultsDirectoryDefined = argc >= 7;
+        if (input_intermediateResultsDirectoryDefined) {
+            optimizer.setIntermediateResultsDirectoryPath(std::filesystem::path(argv[6]));
+        }
+        optimizer.initializePatterns(fullPattern,
+                                     reducedPattern,
+                                     settings_optimization.variablesRanges);
+        optimizer.optimize(settings_optimization, optimizationResults);
+        optimizationResults.label = optimizationName;
+        optimizationResults.baseTriangleFullPattern.copy(fullPattern);
+        optimizationResults.settings = settings_optimization;
+        auto end = std::chrono::system_clock::now();
+        auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(end - start);
+        optimizationResults.time = elapsed.count() / 1000.0;
+
+        if (!optimizationResults.wasSuccessful) {
+            //          resultsOutputDir = crashedJobsDirPath.string();
+            return 1;
+        }
+        resultsOutputDir = convergedJobsDirPath.string();
+        optimizationResults.save(resultsOutputDir, true);
+        csvFile csv_resultsLocalFile(std::filesystem::path(resultsOutputDir).append("results.csv"),
+                                     true);
+        csvFile csv_results({}, false);
+        std::vector<csvFile *> csvVector{&csv_resultsLocalFile, &csv_results};
+        csv_results << "Name";
+        csv_resultsLocalFile << "Name";
+        optimizationResults.writeHeaderTo(csvVector);
+        settings_optimization.writeHeaderTo(csv_results);
+        csv_results << endrow;
+        csv_resultsLocalFile << endrow;
+        csv_results << std::to_string(fullPattern.EN()) + "#" + pairName;
+        csv_resultsLocalFile << std::to_string(fullPattern.EN()) + "#" + pairName;
+        optimizationResults.writeResultsTo(csvVector);
+        settings_optimization.writeSettingsTo(csv_results);
+        csv_results << endrow;
+        csv_resultsLocalFile << endrow;
+    }
 
 #ifdef POLYSCOPE_DEFINED
   std::vector<std::string> scenarioLabels(
@@ -257,8 +322,8 @@ int main(int argc, char *argv[])
   const double colorAxial = 1;
   const double colorShear = 3;
   const double colorBending = 5;
-  const double colorDome = 7;
+  const double colorDome = 0.1;
-  const double colorSaddle = 9;
+  const double colorSaddle = 0;
   std::vector<double> colors(optimizationResults.objectiveValue.perSimulationScenario_total.size());
   for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
       scenarioLabels[scenarioIndex]
@@ -283,16 +348,24 @@ int main(int argc, char *argv[])
   std::vector<double> y(optimizationResults.objectiveValue.perSimulationScenario_total.size());
   for (int scenarioIndex = 0; scenarioIndex < scenarioLabels.size(); scenarioIndex++) {
       y[scenarioIndex]
-          = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
+          //          = optimizationResults.objectiveValue.perSimulationScenario_rawTranslational[scenarioIndex]
-            + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
+          //            + optimizationResults.objectiveValue.perSimulationScenario_rawRotational[scenarioIndex];
-      //      y[scenarioIndex] = optimizationResults.objectiveValue
+          = optimizationResults.objectiveValue.perSimulationScenario_total_unweighted[scenarioIndex];
-      //                             .perSimulationScenario_total[scenarioIndex];
   }
   std::vector<double> x = matplot::linspace(0, y.size() - 1, y.size());
   std::vector<double> markerSizes(y.size(), 5);
-  SimulationResultsReporter::createPlot("scenario index", "error", x, y, markerSizes, colors);
+  SimulationResultsReporter::createPlot("scenario index",
+                                        "error",
+                                        x,
+                                        y,
+                                        markerSizes,
+                                        colors,
+                                        std::filesystem::path(resultsOutputDir)
+                                            .append("perScenarioObjectiveValues.png"));
   //    optimizationResults.saveMeshFiles();
+  std::cout << "Saved results to:" << resultsOutputDir << std::endl;
   optimizationResults.draw();
+  
 #endif
 
   return 0;

src/reducedmodeloptimizer.cpp

@@ -17,6 +17,7 @@ using namespace ReducedPatternOptimization;
 
 struct GlobalOptimizationVariables
 {
+    //    int totalNumberOfSimulationScenarios;
     std::vector<SimulationResults> fullPatternResults;
     std::vector<double> translationalDisplacementNormalizationValues;
     std::vector<double> rotationalDisplacementNormalizationValues;
@@ -30,10 +31,10 @@ struct GlobalOptimizationVariables
     std::vector<double> objectiveValueHistory;
     std::vector<double> plotColors;
     std::array<double,
-               ReducedPatternOptimization::OptimizationParameterIndex::NumberOfOptimizationParameters>
+               ReducedPatternOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
         parametersInitialValue;
     std::array<double,
-               ReducedPatternOptimization::OptimizationParameterIndex::NumberOfOptimizationParameters>
+               ReducedPatternOptimization::OptimizationParameterIndex::NumberOfOptimizationVariables>
         optimizationInitialValue;
     std::vector<int> simulationScenarioIndices;
     double minY{DBL_MAX};
@@ -60,6 +61,7 @@ struct GlobalOptimizationVariables
     std::vector<double> xMin;
     std::vector<double> xMax;
     std::vector<double> scenarioWeights;
+    std::vector<ReducedPatternOptimization::Settings::ObjectiveWeights> objectiveWeights;
 } global;
 
 double ReducedModelOptimizer::computeDisplacementError(
@@ -136,7 +138,8 @@ double ReducedModelOptimizer::computeError(
         &reducedToFullInterfaceViMap,
     const double &normalizationFactor_translationalDisplacement,
     const double &normalizationFactor_rotationalDisplacement,
-    const double &scenarioWeight)
+    const double &scenarioWeight,
+    const Settings::ObjectiveWeights &objectiveWeights)
 {
     const double translationalError
         = computeDisplacementError(simulationResults_fullPattern.displacements,
@@ -149,21 +152,21 @@ double ReducedModelOptimizer::computeError(
     //                                       reducedToFullInterfaceViMap);
 
     //    std::cout << "normalization factor:" << normalizationFactor_rotationalDisplacement << std::endl;
-//        std::cout << "trans error:" << translationalError << std::endl;
+    //        std::cout << "trans error:" << translationalError << std::endl;
     const double rotationalError
         = computeRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
                                  simulationResults_reducedPattern.rotationalDisplacementQuaternion,
                                  reducedToFullInterfaceViMap,
                                  normalizationFactor_rotationalDisplacement);
-    //        computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
+    //    const double rotationalError
-    //                                  simulationResults_reducedPattern.rotationalDisplacementQuaternion,
+    //        = computeRawRotationalError(simulationResults_fullPattern.rotationalDisplacementQuaternion,
-    //                                  reducedToFullInterfaceViMap);
+    //                                    simulationResults_reducedPattern.rotationalDisplacementQuaternion,
+    //                                    reducedToFullInterfaceViMap);
     //        std::cout << "rot error:" << rotationalError<< std::endl;
-    const double simulationError = global.optimizationSettings.objectiveWeights.translational
+    const double simulationError = objectiveWeights.translational * translationalError
-                                       * translationalError
+                                   + objectiveWeights.rotational * rotationalError;
-                                   + global.optimizationSettings.objectiveWeights.rotational
+    assert(!std::isnan(simulationError));
-                                         * rotationalError;
+    return simulationError * simulationError * scenarioWeight * scenarioWeight;
-    return simulationError * simulationError * scenarioWeight;
 }
 
 //double ReducedModelOptimizer::objective(double E,double A,double J,double I2,double I3,
@@ -236,7 +239,7 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
     //    global.reducedPatternSimulationJobs[0]->pMesh->unregister();
 
     // run simulations
-    std::vector<double> simulationErrorsPerScenario(global.simulationScenarioIndices.size());
+    std::vector<double> simulationErrorsPerScenario(totalNumberOfSimulationScenarios, 0);
     LinearSimulationModel simulator;
     simulator.setStructure(pReducedPatternSimulationMesh);
     //    simulator.initialize();
@@ -287,21 +290,22 @@ double ReducedModelOptimizer::objective(const std::vector<double> &x)
                     global.reducedToFullInterfaceViMap,
                     global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
                     global.rotationalDisplacementNormalizationValues[simulationScenarioIndex],
-                    scenarioWeight);
+                    scenarioWeight,
+                    global.objectiveWeights[simulationScenarioIndex]);
 
                 simulationErrorsPerScenario[simulationScenarioIndex] = simulationScenarioError;
                 //            }
-//                #ifdef POLYSCOPE_DEFINED
+                //                #ifdef POLYSCOPE_DEFINED
-//                            reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
+                //                            reducedJob->pMesh->registerForDrawing(Colors::reducedInitial);
-//                            reducedModelResults.registerForDrawing(Colors::reducedDeformed);
+                //                            reducedModelResults.registerForDrawing(Colors::reducedDeformed);
-//                            global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
+                //                            global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
-//                            global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
+                //                            global.fullPatternResults[simulationScenarioIndex].registerForDrawing(
-//                                Colors::fullDeformed);
+                //                                Colors::fullDeformed);
-//                            polyscope::show();
+                //                            polyscope::show();
-//                            reducedModelResults.unregister();
+                //                            reducedModelResults.unregister();
-//                            global.pFullPatternSimulationMesh->unregister();
+                //                            global.pFullPatternSimulationMesh->unregister();
-//                            global.fullPatternResults[simulationScenarioIndex].unregister();
+                //                            global.fullPatternResults[simulationScenarioIndex].unregister();
-//                #endif
+                //                #endif
             }
         });
     //    double totalError = std::accumulate(simulationErrorsPerScenario.begin(),
@@ -361,7 +365,7 @@ void ReducedModelOptimizer::createSimulationMeshes(
 {
     if (typeid(CrossSectionType) != typeid(RectangularBeamDimensions)) {
         std::cerr << "Error: A rectangular cross section is expected." << std::endl;
-        terminate();
+        std::terminate();
     }
     // Full pattern
     pFullPatternSimulationMesh = std::make_shared<SimulationMesh>(fullModel);
@@ -542,9 +546,10 @@ ReducedModelOptimizer::ReducedModelOptimizer(const std::vector<size_t> &numberOf
     scenarioIsSymmetrical[BaseSimulationScenario::Saddle] = true;
 }
 
-void ReducedModelOptimizer::initializePatterns(PatternGeometry &fullPattern,
+void ReducedModelOptimizer::initializePatterns(
-                                               PatternGeometry &reducedPattern,
+    PatternGeometry &fullPattern,
-                                               const std::vector<xRange> &optimizationParameters)
+    PatternGeometry &reducedPattern,
+    const std::array<xRange, NumberOfOptimizationVariables> &optimizationParameters)
 {
     assert(fullPattern.VN() == reducedPattern.VN() && fullPattern.EN() >= reducedPattern.EN());
     fullPatternNumberOfEdges = fullPattern.EN();
@@ -606,6 +611,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
 
     global.functions_updateReducedPatternParameter[E] =
         [](const double &newE, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+            //            std::cout << "Updating E with new value:" << newE << std::endl;
             for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
                 Element &e = pReducedPatternSimulationMesh->elements[ei];
                 e.setMaterial(ElementMaterial(e.material.poissonsRatio, newE));
@@ -613,6 +619,7 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
         };
     global.functions_updateReducedPatternParameter[A] =
         [](const double &newA, std::shared_ptr<SimulationMesh> &pReducedPatternSimulationMesh) {
+            assert(pReducedPatternSimulationMesh != nullptr);
             const double beamWidth = std::sqrt(newA);
             const double beamHeight = beamWidth;
             for (EdgeIndex ei = 0; ei < pReducedPatternSimulationMesh->EN(); ei++) {
@@ -648,15 +655,18 @@ void ReducedModelOptimizer::initializeUpdateReducedPatternFunctions()
 }
 
 void ReducedModelOptimizer::initializeOptimizationParameters(
-    const std::shared_ptr<SimulationMesh> &mesh, const std::vector<xRange> &optimizationParameters)
+    const std::shared_ptr<SimulationMesh> &mesh,
+    const std::array<ReducedPatternOptimization::xRange,
+                     ReducedPatternOptimization::NumberOfOptimizationVariables>
+        &optimizationParameters)
 {
-    global.numberOfOptimizationParameters = NumberOfOptimizationParameters;
+    global.numberOfOptimizationParameters = NumberOfOptimizationVariables;
 
     for (int optimizationParameterIndex = 0;
          optimizationParameterIndex < optimizationParameters.size();
          optimizationParameterIndex++) {
         for (int optimizationParameterIndex = E;
-             optimizationParameterIndex != NumberOfOptimizationParameters;
+             optimizationParameterIndex != NumberOfOptimizationVariables;
              optimizationParameterIndex++) {
             //            const xRange &xOptimizationParameter = optimizationParameters[optimizationParameterIndex];
             switch (optimizationParameterIndex) {
@@ -764,24 +774,19 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
                   << std::endl;
     }
     //Optimal X values
-    results.optimalXNameValuePairs.resize(NumberOfOptimizationParameters);
+    results.optimalXNameValuePairs.resize(NumberOfOptimizationVariables);
     for (int optimizationParameterIndex = E;
-         optimizationParameterIndex != NumberOfOptimizationParameters;
+         optimizationParameterIndex != NumberOfOptimizationVariables;
          optimizationParameterIndex++) {
-        if (settings.parameterRanges[optimizationParameterIndex].max
-                - settings.parameterRanges[optimizationParameterIndex].min
-            < 1e-10) {
-            continue;
-        }
         if (optimizationParameterIndex != OptimizationParameterIndex::R && optimizationParameterIndex != OptimizationParameterIndex::Theta ) {
             results.optimalXNameValuePairs[optimizationParameterIndex]
-                = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label,
+                = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label,
                                  optimalObjective.x[optimizationParameterIndex]
                                      * global.parametersInitialValue[optimizationParameterIndex]);
         } else {
             //Hex size and angle are pure values (not multipliers of the initial values)
             results.optimalXNameValuePairs[optimizationParameterIndex]
-                = std::make_pair(settings.parameterRanges[optimizationParameterIndex].label,
+                = std::make_pair(settings.variablesRanges[optimizationParameterIndex].label,
                                  optimalObjective.x[optimizationParameterIndex]);
         }
 
@@ -815,6 +820,8 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
         global.simulationScenarioIndices.size());
     results.objectiveValue.perSimulationScenario_total.resize(
         global.simulationScenarioIndices.size());
+    results.objectiveValue.perSimulationScenario_total_unweighted.resize(
+        global.simulationScenarioIndices.size());
     //#ifdef POLYSCOPE_DEFINED
     //    global.pFullPatternSimulationMesh->registerForDrawing(Colors::fullDeformed);
     //#endif
@@ -840,7 +847,17 @@ void ReducedModelOptimizer::getResults(const FunctionEvaluation &optimalObjectiv
             global.reducedToFullInterfaceViMap,
             global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
             global.rotationalDisplacementNormalizationValues[simulationScenarioIndex],
-            scenarioWeight);
+            scenarioWeight,
+            global.objectiveWeights[simulationScenarioIndex]);
+
+        results.objectiveValue.perSimulationScenario_total_unweighted[i] = computeError(
+            global.fullPatternResults[simulationScenarioIndex],
+            reducedModelResults,
+            global.reducedToFullInterfaceViMap,
+            global.translationalDisplacementNormalizationValues[simulationScenarioIndex],
+            global.rotationalDisplacementNormalizationValues[simulationScenarioIndex],
+            1,
+            global.objectiveWeights[simulationScenarioIndex]);
 
         //        results.objectiveValue.total += results.objectiveValue.perSimulationScenario_total[i];
         //Raw translational
@@ -995,9 +1012,9 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
 #if POLYSCOPE_DEFINED
 //    global.plotColors.reserve(settings.numberOfFunctionCalls);
 #endif
-    assert(!settings.optimizationVariables.empty());
+    assert(!settings.optimizationStrategy.empty());
     const std::vector<std::vector<OptimizationParameterIndex>> &optimizationParametersGroups
-        = settings.optimizationVariables;
+        = settings.optimizationStrategy;
 
 #ifndef USE_ENSMALLEN
     const int totalNumberOfOptimizationParameters
@@ -1011,15 +1028,16 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
 #endif
 
     FunctionEvaluation optimization_optimalResult;
-    optimization_optimalResult.x.resize(NumberOfOptimizationParameters,0);
+    optimization_optimalResult.x.resize(NumberOfOptimizationVariables, 0);
-        for (int optimizationParameterIndex = E;
+    for (int optimizationParameterIndex = E;
-             optimizationParameterIndex != NumberOfOptimizationParameters;
+         optimizationParameterIndex != NumberOfOptimizationVariables;
-             optimizationParameterIndex++) {
+         optimizationParameterIndex++) {
-            optimization_optimalResult.x[optimizationParameterIndex]
+        optimization_optimalResult.x[optimizationParameterIndex]
-                = global.parametersInitialValue[optimizationParameterIndex];
+            = global.optimizationInitialValue[optimizationParameterIndex];
-        }
+    }
-        for (const std::vector<OptimizationParameterIndex> &parameterGroup :
+        for (int parameterGroupIndex=0;parameterGroupIndex<optimizationParametersGroups.size();parameterGroupIndex++) {
-             optimizationParametersGroups) {
+        const std::vector<OptimizationParameterIndex> &parameterGroup =
+             optimizationParametersGroups[parameterGroupIndex];
             FunctionEvaluation parameterGroup_optimalResult;
             //Set update function. TODO: Make this function immutable by defining it once and using the global variable to set parameterGroup
             function_updateReducedPattern = [&](const std::vector<double> &x,
@@ -1036,8 +1054,8 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
                             = global.parametersInitialValue[parameterIndex];
                         return x[xIndex] * parameterInitialValue;
                     }();
-                    //                std::cout << "Optimization parameter:" << parameterIndex << std::endl;
+                    //                                    std::cout << "Optimization parameter:" << parameterIndex << std::endl;
-                    //                std::cout << "New value:" << parameterNewValue << std::endl;
+                    //                                    std::cout << "New value:" << parameterNewValue << std::endl;
                     global.functions_updateReducedPatternParameter[parameterIndex](parameterNewValue,
                                                                                    pMesh);
                 }
@@ -1051,8 +1069,8 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
             for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
                 const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
 
-                xMin[xIndex] = settings.parameterRanges[parameterIndex].min;
+                xMin[xIndex] = settings.variablesRanges[parameterIndex].min;
-                xMax[xIndex] = settings.parameterRanges[parameterIndex].max;
+                xMax[xIndex] = settings.variablesRanges[parameterIndex].max;
             }
 #ifdef USE_ENSMALLEN
         arma::mat x(parameterGroup.size(), 1);
@@ -1070,14 +1088,21 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
         // that is able to handle arbitrary functions.
         EnsmallenOptimizationObjective optimizationFunction;
         //Set min max values
-        //        ens::SA optimizer;
+        //                ens::SA optimizer;
-        //        ens::CNE optimizer;
+        //                ens::CNE optimizer;
-        //        ens::DE optimizer;
+        //                ens::DE optimizer;
-        //        ens::SPSA optimizer;
+        //                ens::SPSA optimizer;
         //        arma::mat xMin_arma(global.xMin);
         //        arma::mat xMax_arma(global.xMax);
-        //                ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000);
+        //        ens::LBestPSO optimizer(64, xMin_arma, xMax_arma, 1000);
-        ens::LBestPSO optimizer(64, 1, 1, 3000, 350, 1e-1, 2.05, 2.15);
+        ens::LBestPSO optimizer(64,
+                                1,
+                                1,
+                                settings.numberOfFunctionCalls,
+                                500,
+                                settings.solverAccuracy,
+                                2.05,
+                                2.35);
         //        ens::LBestPSO optimizer;
         const double minima = optimizer.Optimize(optimizationFunction, x);
         //        for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
@@ -1106,12 +1131,22 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
         for (int xIndex = 0; xIndex < parameterGroup.size(); xIndex++) {
             const OptimizationParameterIndex parameterIndex = parameterGroup[xIndex];
 
-            xMin_dlib(xIndex) = settings.parameterRanges[parameterIndex].min;
+            xMin_dlib(xIndex) = settings.variablesRanges[parameterIndex].min;
-            xMax_dlib(xIndex) = settings.parameterRanges[parameterIndex].max;
+            xMax_dlib(xIndex) = settings.variablesRanges[parameterIndex].max;
         }
-        const int optimizationNumberOfFunctionCalls = static_cast<int>(
+        const double portionOfTotalFunctionCallsBudget = [&]() {
-            (static_cast<double>(parameterGroup.size()) / totalNumberOfOptimizationParameters)
+            if (settings.optimizationVariablesGroupsWeights.has_value()) {
-            * settings.numberOfFunctionCalls);
+                assert(settings.optimizationVariablesGroupsWeights->size()
+                       == settings.optimizationStrategy.size());
+                return settings.optimizationVariablesGroupsWeights.value()[parameterGroupIndex];
+
+            } else {
+                return static_cast<double>(parameterGroup.size())
+                       / totalNumberOfOptimizationParameters;
+            }
+        }();
+        const int optimizationNumberOfFunctionCalls = portionOfTotalFunctionCallsBudget
+                                                      * settings.numberOfFunctionCalls;
         const dlib::function_evaluation optimalResult_dlib = [&]() {
             if (parameterGroup.size() == 1) {
                 dlib::function_evaluation result;
@@ -1181,9 +1216,16 @@ void ReducedModelOptimizer::runOptimization(const Settings &settings,
 
 #endif
 
+        const auto firstNonNullReducedSimulationJobIt
+            = std::find_if(global.reducedPatternSimulationJobs.begin(),
+                           global.reducedPatternSimulationJobs.end(),
+                           [](const std::shared_ptr<SimulationJob> &pJob) {
+                               return pJob != nullptr;
+                           });
         function_updateReducedPattern(
-            std::vector<double>(parameterGroup_optimalResult.x.begin(), parameterGroup_optimalResult.x.end()),
+            std::vector<double>(parameterGroup_optimalResult.x.begin(),
-            global.reducedPatternSimulationJobs[0]
+                                parameterGroup_optimalResult.x.end()),
+            (*firstNonNullReducedSimulationJobIt)
                 ->pMesh); //TODO: Check if its ok to update only the changed parameters
 #ifdef POLYSCOPE_DEFINED
         std::cout << "Optimal x:"
@@ -1203,7 +1245,7 @@ optimization_optimalResult.y=parameterGroup_optimalResult.y;
             optimization_optimalResult.x[parameterIndex] = parameterGroup_optimalResult.x[xIndex];
         }
         }
-    getResults(optimization_optimalResult, settings, results);
+        getResults(optimization_optimalResult, settings, results);
 }
 
 void ReducedModelOptimizer::constructAxialSimulationScenario(
@@ -1359,7 +1401,7 @@ double fullPatternMaxSimulationForceRotationalObjective(const double &forceMagni
     const double &desiredRotationAngle = global.desiredMaxRotationAngle;
     double error;
 #ifdef POLYSCOPE_DEFINED
-    job.pMesh->setLabel("initial");
+    //    job.pMesh->setLabel("initial");
     //    job.pMesh->registerForDrawing();
     //    results.registerForDrawing();
     //    polyscope::show();
@@ -1587,55 +1629,66 @@ double ReducedModelOptimizer::computeFullPatternMaxSimulationForce(
 }
 
 void ReducedModelOptimizer::computeScenarioWeights(
-    const std::vector<BaseSimulationScenario> &baseSimulationScenarios,
+    const std::vector<BaseSimulationScenario> &baseSimulationScenarios)
-    const std::array<double, 5> &maxForceMagnitudePerBaseScenario)
 {
-    assert(maxForceMagnitudePerBaseScenario.size() == NumberOfBaseSimulationScenarios);
+    std::array<double, NumberOfBaseSimulationScenarios> baseScenario_equalizationWeights;
-    const double userWeightsSum = std::accumulate(baseScenarioWeights.begin(),
+    baseScenario_equalizationWeights.fill(0);
-                                                  baseScenarioWeights.end(),
+    std::array<double, NumberOfBaseSimulationScenarios> baseScenario_userWeights;
-                                                  0);
+    baseScenario_userWeights.fill(0);
-    const int numberOfBaseScenarios = baseSimulationScenarios.size();
+    //Fill only used base scenario weights
-
-    Eigen::VectorXd baseScenario_equalizationWeights(numberOfBaseScenarios);
-    Eigen::VectorXd baseScenario_userWeights(numberOfBaseScenarios);
-    Eigen::Matrix<double, NumberOfBaseSimulationScenarios, 1> baseScenario_maxMagnitudeWeights(
-        maxForceMagnitudePerBaseScenario.data());
     for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) {
-        const double baseScenarioEqualizationWeight
+        baseScenario_equalizationWeights[baseScenario] = static_cast<double>(
-            = totalNumberOfSimulationScenarios
+            simulationScenariosResolution[baseScenario]);
-              / static_cast<double>(simulationScenariosResolution[baseScenario]);
+        baseScenario_userWeights[baseScenario] = baseScenarioWeights[baseScenario];
-        baseScenario_equalizationWeights(baseScenario) = baseScenarioEqualizationWeight;
-        const double baseScenarioUserWeight = static_cast<double>(baseScenarioWeights[baseScenario])
-                                              / userWeightsSum;
-        baseScenario_userWeights(baseScenario) = baseScenarioUserWeight;
     }
 
-    baseScenario_maxMagnitudeWeights.normalize();
+    Utilities::normalize(baseScenario_userWeights.begin(), baseScenario_userWeights.end());
-    baseScenario_equalizationWeights.normalize();
-    baseScenario_userWeights.normalize();
 
+    Utilities::normalize(baseScenario_equalizationWeights.begin(),
+                         baseScenario_equalizationWeights.end());
+    std::transform(baseScenario_equalizationWeights.begin(),
+                   baseScenario_equalizationWeights.end(),
+                   baseScenario_equalizationWeights.begin(),
+                   [](const double &d) {
+                       if (d == 0) {
+                           return d;
+                       }
+                       return 1 / d;
+                   });
+    Utilities::normalize(baseScenario_equalizationWeights.begin(),
+                         baseScenario_equalizationWeights.end());
+
+    std::array<double, NumberOfBaseSimulationScenarios> baseScenarios_weights;
+    baseScenarios_weights.fill(0);
+    for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) {
+        baseScenarios_weights[baseScenario] = baseScenario_equalizationWeights[baseScenario]
+                                              + 2 * baseScenario_userWeights[baseScenario];
+    }
+    Utilities::normalize(baseScenarios_weights.begin(), baseScenarios_weights.end());
+
+    global.objectiveWeights.resize(totalNumberOfSimulationScenarios);
     global.scenarioWeights.resize(
-        totalNumberOfSimulationScenarios); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization
+        totalNumberOfSimulationScenarios,
+        0); //This essentially holds the base scenario weights but I use totalNumberOfSimulationScenarios elements instead in order to have O(1) access time during the optimization
     for (const BaseSimulationScenario &baseScenario : baseSimulationScenarios) {
         const int baseSimulationScenarioIndexOffset
             = std::accumulate(simulationScenariosResolution.begin(),
                               simulationScenariosResolution.begin() + baseScenario,
                               0);
-        const double baseScenarioEqualizationWeight = baseScenario_equalizationWeights(baseScenario);
-        //        const double baseScenarioUserWeight = baseScenario_userWeights(baseScenario);
-        const double baseScenarioUserWeight = baseScenario_userWeights(baseScenario);
-
-        const double baseScenarioMaxMagnitudeWeight = baseScenario_maxMagnitudeWeights(baseScenario);
         for (int simulationScenarioIndex = 0;
              simulationScenarioIndex < simulationScenariosResolution[baseScenario];
              simulationScenarioIndex++) {
             const int globalScenarioIndex = baseSimulationScenarioIndexOffset
                                             + simulationScenarioIndex;
-            const double scenarioWeight = /*baseScenarioUserWeight +*/ baseScenarioEqualizationWeight
+            global.scenarioWeights[globalScenarioIndex] = baseScenarios_weights[baseScenario];
-                                          + baseScenarioMaxMagnitudeWeight;
+
-            global.scenarioWeights[globalScenarioIndex] = scenarioWeight;
+            global.objectiveWeights[globalScenarioIndex]
+                = global.optimizationSettings.perBaseScenarioObjectiveWeights[baseScenario];
         }
     }
+
+    //Dont normalize since we want the base scenarios to be normalized not the "global scenarios"
+    //    Utilities::normalize(global.scenarioWeights.begin(), global.scenarioWeights.end());
 }
 
 std::vector<std::shared_ptr<SimulationJob>> ReducedModelOptimizer::createFullPatternSimulationJobs(
@@ -1771,12 +1824,25 @@ void ReducedModelOptimizer::computeObjectiveValueNormalizationFactors()
     }
 }
 
+void ReducedModelOptimizer::setIntermediateResultsDirectoryPath(
+    const std::filesystem::path &newIntermediateResultsDirectoryPath)
+{
+    intermediateResultsDirectoryPath = newIntermediateResultsDirectoryPath;
+}
+
 void ReducedModelOptimizer::optimize(
     const Settings &optimizationSettings,
     ReducedPatternOptimization::Results &results,
     const std::vector<BaseSimulationScenario> &desiredBaseSimulationScenarioIndices)
 {
-    assert(!optimizationSettings.optimizationVariables.empty());
+    assert(!optimizationSettings.optimizationStrategy.empty());
+    assert(!optimizationSettings.optimizationVariablesGroupsWeights.has_value()
+           || (optimizationSettings.optimizationStrategy.size()
+                   == optimizationSettings.optimizationVariablesGroupsWeights->size()
+               && std::accumulate(optimizationSettings.optimizationVariablesGroupsWeights->begin(),
+                                  optimizationSettings.optimizationVariablesGroupsWeights->end(),
+                                  0)));
+
     for (int baseSimulationScenarioIndex : desiredBaseSimulationScenarioIndices) {
         //Increase the size of the vector holding the simulation scenario indices
         global.simulationScenarioIndices.resize(
@@ -1792,6 +1858,7 @@ void ReducedModelOptimizer::optimize(
                                   0));
     }
 
+    //    global.totalNumberOfSimulationScenarios = totalNumberOfSimulationScenarios;
     global.reducedPatternSimulationJobs.resize(totalNumberOfSimulationScenarios);
     global.fullPatternResults.resize(totalNumberOfSimulationScenarios);
     global.translationalDisplacementNormalizationValues.resize(totalNumberOfSimulationScenarios);
@@ -1802,9 +1869,6 @@ void ReducedModelOptimizer::optimize(
     global.optimizationSettings = optimizationSettings;
     global.pFullPatternSimulationMesh = m_pFullPatternSimulationMesh;
 
-    const std::filesystem::path intermediateResultsDirectoryPath(
-        optimizationSettings.intermediateResultsDirectoryPath);
-
     std::array<double, NumberOfBaseSimulationScenarios> fullPatternSimulationScenarioMaxMagnitudes
         = getFullPatternMaxSimulationForces(desiredBaseSimulationScenarioIndices,
                                             intermediateResultsDirectoryPath);
@@ -1813,14 +1877,14 @@ void ReducedModelOptimizer::optimize(
                                           fullPatternSimulationScenarioMaxMagnitudes);
     //  polyscope::removeAllStructures();
 
-    computeScenarioWeights(desiredBaseSimulationScenarioIndices,
+    //    ComputeScenarioWeights({Axial, Shear, Bending, Dome, Saddle},
-                           fullPatternSimulationScenarioMaxMagnitudes);
+    computeScenarioWeights(desiredBaseSimulationScenarioIndices);
 
     results.baseTriangle = global.baseTriangle;
 
     DRMSimulationModel::Settings simulationSettings;
     simulationSettings.maxDRMIterations = 200000;
-    simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-8;
+    simulationSettings.totalTranslationalKineticEnergyThreshold = 1e-10;
     simulationSettings.viscousDampingFactor = 5e-3;
     simulationSettings.useKineticDamping = true;
 
@@ -1848,8 +1912,9 @@ void ReducedModelOptimizer::optimize(
                 .append(m_pFullPatternSimulationMesh->getLabel())
                 .append(pFullPatternSimulationJob->getLabel()));
         //                .append(pFullPatternSimulationJob->getLabel() + ".json")
+        constexpr bool recomputeFullPatternResults = false;
         SimulationResults fullPatternResults;
-        if (std::filesystem::exists(jobResultsDirectoryPath)) {
+        if (!recomputeFullPatternResults && std::filesystem::exists(jobResultsDirectoryPath)) {
             fullPatternResults.load(std::filesystem::path(jobResultsDirectoryPath).append("Results"),
                                     pFullPatternSimulationJob);
         } else {

src/reducedmodeloptimizer.hpp

@@ -41,6 +41,7 @@ class ReducedModelOptimizer
     std::vector<bool> scenarioIsSymmetrical;
     int fullPatternNumberOfEdges;
     constexpr static double youngsModulus{1 * 1e9};
+    std::filesystem::path intermediateResultsDirectoryPath;
 
 public:
     struct FunctionEvaluation
@@ -64,13 +65,12 @@ public:
     //    };
     //    inline constexpr static ParameterLabels parameterLabels();
 
-    inline static std::array<std::string, ReducedPatternOptimization::NumberOfOptimizationParameters>
+    inline static std::array<std::string, ReducedPatternOptimization::NumberOfOptimizationVariables>
         parameterLabels = {"R", "A", "I2", "I3", "J", "Theta", "R"};
-    constexpr static std::array<int, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
+    constexpr static std::array<double, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
         simulationScenariosResolution = {11, 11, 20, 20, 20};
-    constexpr static std::array<int, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
+    constexpr static std::array<double, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
-        baseScenarioWeights = {1, 1, 100, 50, 50};
+        baseScenarioWeights = {1, 1, 2, 3, 2};
-    //    constexpr static std::array<int, 5> simulationScenariosResolution = {3, 3, 3, 3, 3};
     inline static int totalNumberOfSimulationScenarios
         = std::accumulate(simulationScenariosResolution.begin(),
                           simulationScenariosResolution.end(),
@@ -98,11 +98,6 @@ public:
 
     SimulationJob getReducedSimulationJob(const SimulationJob &fullModelSimulationJob);
 
-    void initializePatterns(
-        PatternGeometry &fullPattern,
-        PatternGeometry &reducedPatterm,
-        const std::vector<ReducedPatternOptimization::xRange> &optimizationParameters);
-
     static void runSimulation(const std::string &filename, std::vector<double> &x);
 
     static std::vector<std::shared_ptr<SimulationJob>> createFullPatternSimulationJobs(
@@ -170,7 +165,8 @@ public:
             &reducedToFullInterfaceViMap,
         const double &normalizationFactor_translationalDisplacement,
         const double &normalizationFactor_rotationalDisplacement,
-        const double &scenarioWeight);
+        const double &scenarioWeight,
+        const ReducedPatternOptimization::Settings::ObjectiveWeights &objectiveWeights);
     static void constructAxialSimulationScenario(
         const double &forceMagnitude,
         const std::vector<std::pair<FullPatternVertexIndex, FullPatternVertexIndex>>
@@ -225,6 +221,16 @@ public:
     static void initializeUpdateReducedPatternFunctions();
     static double objective(const double &xValue);
 
+    void initializePatterns(
+        PatternGeometry &fullPattern,
+        PatternGeometry &reducedPattern,
+        const std::array<ReducedPatternOptimization::xRange,
+                         ReducedPatternOptimization::NumberOfOptimizationVariables>
+            &optimizationParameters);
+
+    void setIntermediateResultsDirectoryPath(
+        const std::filesystem::path &newIntermediateResultsDirectoryPath);
+
 private:
     static void computeDesiredReducedModelDisplacements(
         const SimulationResults &fullModelResults,
@@ -236,7 +242,9 @@ private:
     void createSimulationMeshes(PatternGeometry &fullModel, PatternGeometry &reducedModel);
     static void initializeOptimizationParameters(
         const std::shared_ptr<SimulationMesh> &mesh,
-        const std::vector<ReducedPatternOptimization::xRange> &optimizationParamters);
+        const std::array<ReducedPatternOptimization::xRange,
+                         ReducedPatternOptimization::NumberOfOptimizationVariables>
+            &optimizationParamters);
 
     DRMSimulationModel simulator;
     void computeObjectiveValueNormalizationFactors();
@@ -256,8 +264,7 @@ private:
             &desiredBaseSimulationScenarioIndices,
         const std::filesystem::path &intermediateResultsDirectoryPath);
     void computeScenarioWeights(const std::vector<ReducedPatternOptimization::BaseSimulationScenario>
-                                    &baseSimulationScenarios,
+                                    &baseSimulationScenarios);
-                                const std::array<double, 5> &maxForceMagnitudePerBaseScenario);
     std::array<double, ReducedPatternOptimization::NumberOfBaseSimulationScenarios>
     computeFullPatternMaxSimulationForces(
         const std::vector<ReducedPatternOptimization::BaseSimulationScenario>