.. _program_listing_file_include_parpeamici_hierarchicalOptimization.h:

Program Listing for File hierarchicalOptimization.h
===================================================

|exhale_lsh| :ref:`Return to documentation for file <file_include_parpeamici_hierarchicalOptimization.h>` (``include/parpeamici/hierarchicalOptimization.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #ifndef HIERARCHICALOPTIMIZATION_H
   #define HIERARCHICALOPTIMIZATION_H
   
   #include <parpeamici/hierarchicalOptimizationAnalyticalParameterProvider.h>
   #include <parpeamici/multiConditionProblem.h>
   #include <parpeoptimization/optimizationProblem.h>
   
   #include <gsl/gsl-lite.hpp>
   
   #include <memory>
   
   
   namespace parpe {
   
   enum class ErrorModel
   {
       normal
   }; // TODO logNormal, Laplace
   
   class AnalyticalParameterProvider;
   class AnalyticalParameterHdf5Reader;
   class HierarchicalOptimizationProblemWrapper;
   class HierarchicalOptimizationWrapper;
   
   class HierarchicalOptimizationWrapper : public GradientFunction
   {
     public:
       HierarchicalOptimizationWrapper(
           AmiciSummedGradientFunction *wrapped_function,
           int numConditions = 0,
           int numObservables = 0);
   
       HierarchicalOptimizationWrapper(
           AmiciSummedGradientFunction *wrapped_function,
           const H5::H5File& file,
           const std::string& hdf5RootPath,
           int numConditions,
           int numObservables,
           ErrorModel errorModel);
   
       HierarchicalOptimizationWrapper(
           AmiciSummedGradientFunction *wrapped_function,
           std::unique_ptr<AnalyticalParameterProvider> scalingReader,
           std::unique_ptr<AnalyticalParameterProvider> offsetReader,
           std::unique_ptr<AnalyticalParameterProvider> sigmaReader,
           int numConditions,
           int numObservables,
           ErrorModel errorModel);
   
       using GradientFunction::evaluate;
   
       FunctionEvaluationStatus evaluate(gsl::span<double const> parameters,
                                         double& fval,
                                         gsl::span<double> gradient,
                                         Logger* logger,
                                         double* cpuTime) const override;
   
       FunctionEvaluationStatus evaluate(gsl::span<double const> reducedParameters,
                                         double& fval,
                                         gsl::span<double> gradient,
                                         std::vector<double>& fullParameters,
                                         std::vector<double>& fullGradient,
                                         Logger* logger,
                                         double* cpuTime) const;
   
       [[nodiscard]] std::vector<double> getDefaultScalingFactors() const;
   
       [[nodiscard]] std::vector<double> getDefaultOffsetParameters() const;
   
       [[nodiscard]] std::vector<double> getDefaultSigmaParameters() const;
   
       [[nodiscard]] std::tuple<std::vector<std::vector<double>>,
                  std::vector<std::vector<double>>>
       getUnscaledModelOutputsAndSigmas(
           const gsl::span<double const> reducedParameters,
           Logger* logger,
           double* cpuTime) const;
   
       [[nodiscard]] std::vector<double> computeAnalyticalScalings(
           std::vector<std::vector<double>> const& measurements,
           std::vector<std::vector<double>> const& modelOutputsUnscaled) const;
   
       void applyOptimalScalings(
           std::vector<double> const& proportionalityFactors,
           std::vector<std::vector<double>>& modelOutputs) const;
   
       [[nodiscard]] std::vector<double> computeAnalyticalOffsets(
           const std::vector<std::vector<double>>& measurements,
           std::vector<std::vector<double>>& modelOutputsUnscaled) const;
   
       [[nodiscard]] std::vector<double> computeAnalyticalSigmas(
           std::vector<std::vector<double>> const& measurements,
           const std::vector<std::vector<double>>& modelOutputsScaled) const;
   
       void applyOptimalOffsets(
           std::vector<double> const& offsetParameters,
           std::vector<std::vector<double>>& modelOutputs) const;
   
       void fillInAnalyticalSigmas(
           std::vector<std::vector<double>>& allSigmas,
           const std::vector<double>& analyticalSigmas) const;
   
       virtual FunctionEvaluationStatus evaluateWithOptimalParameters(
           std::vector<double> const& fullParameters,
           std::vector<double> const& sigmas,
           std::vector<std::vector<double>> const& measurements,
           std::vector<std::vector<double>> const& modelOutputsScaled,
           std::vector<std::vector<double> > &fullSigmaMatrices,
           double& fval,
           const gsl::span<double> gradient,
           std::vector<double>& fullGradient,
           Logger* logger,
           double* cpuTime) const;
   
       [[nodiscard]] int numParameters() const override;
   
       [[nodiscard]] int numProportionalityFactors() const;
   
       [[nodiscard]] std::vector<int> const& getProportionalityFactorIndices() const;
   
       [[nodiscard]] int numOffsetParameters() const;
   
       [[nodiscard]] int numSigmaParameters() const;
   
       [[nodiscard]] std::vector<int> const& getOffsetParameterIndices() const;
   
       [[nodiscard]] std::vector<int> const& getSigmaParameterIndices() const;
   
       [[nodiscard]] std::vector<int> getAnalyticalParameterIndices() const;
   
       [[nodiscard]] AmiciSummedGradientFunction* getWrappedFunction() const;
   
       [[nodiscard]] std::vector<std::string> getParameterIds() const override;
   
     private:
       void init();
       AmiciSummedGradientFunction *wrapped_function_;
   
       std::unique_ptr<AnalyticalParameterProvider> scalingReader;
       std::unique_ptr<AnalyticalParameterProvider> offsetReader;
       std::unique_ptr<AnalyticalParameterProvider> sigmaReader;
   
       std::vector<int> proportionalityFactorIndices;
       std::vector<int> offsetParameterIndices;
       std::vector<int> sigmaParameterIndices;
   
       int numConditions;
       int numObservables;
   
       ErrorModel errorModel = ErrorModel::normal;
   };
   
   
   class HierarchicalOptimizationProblemWrapper : public OptimizationProblem
   {
     public:
       HierarchicalOptimizationProblemWrapper() = default;
   
       HierarchicalOptimizationProblemWrapper(
           std::unique_ptr<OptimizationProblem> problemToWrap,
           const MultiConditionDataProviderHDF5* dataProvider);
   
       HierarchicalOptimizationProblemWrapper(
           std::unique_ptr<OptimizationProblem> problemToWrap,
           std::unique_ptr<HierarchicalOptimizationWrapper> costFun,
           std::unique_ptr<Logger> logger);
   
       HierarchicalOptimizationProblemWrapper(
           HierarchicalOptimizationProblemWrapper const& other) = delete;
   
       void fillInitialParameters(gsl::span<double> buffer) const override;
   
       void fillParametersMin(gsl::span<double> buffer) const override;
   
       void fillParametersMax(gsl::span<double> buffer) const override;
   
       void fillFilteredParams(std::vector<double> const& fullParams,
                               gsl::span<double> buffer) const;
   
       OptimizationOptions const& getOptimizationOptions() const override
       {
           return wrapped_problem_->getOptimizationOptions();
       }
       void setOptimizationOptions(OptimizationOptions const& options) override
       {
           wrapped_problem_->setOptimizationOptions(options);
       }
   
       // TODO: need to ensure that this will work with the reduced number of
       // parameters
       std::unique_ptr<OptimizationReporter> getReporter() const override;
   
     private:
       std::unique_ptr<OptimizationProblem> wrapped_problem_;
   };
   
   class HierarchicalOptimizationReporter : public OptimizationReporter
   {
     public:
       HierarchicalOptimizationReporter(
           HierarchicalOptimizationWrapper* gradFun,
           std::unique_ptr<OptimizationResultWriter> rw,
           std::unique_ptr<Logger> logger);
   
       using GradientFunction::evaluate;
   
       FunctionEvaluationStatus evaluate(gsl::span<const double> parameters,
                                         double& fval,
                                         gsl::span<double> gradient,
                                         Logger* logger = nullptr,
                                         double* cpuTime = nullptr) const override;
   
       // bool starting(gsl::span<const double> initialParameters) const override;
   
       // TODO: always update final parameters
       bool iterationFinished(
           gsl::span<const double> parameters,
           double objectiveFunctionValue,
           gsl::span<const double> objectiveFunctionGradient) const override;
   
       bool afterCostFunctionCall(
           gsl::span<const double> parameters,
           double objectiveFunctionValue,
           gsl::span<double const> objectiveFunctionGradient) const override;
   
       void finished(double optimalCost,
                     gsl::span<const double> parameters,
                     int exitStatus) const override;
   
       std::vector<double> const& getFinalParameters() const override;
   
       HierarchicalOptimizationWrapper* hierarchical_wrapper_ = nullptr;
   
       /* In addition to the vectors for the outer optimization problem,
        * we also keep the complete ones.
        */
       mutable std::vector<double> cached_full_parameters_;
       mutable std::vector<double> cached_full_gradient_;
       // TODO should override other functions as well
       // TODO: in all functions, we need to check of the provided parameters or
       // function values match To cached ones, if we want to provide all together
       // to downstream
   };
   
   void
   fillFilteredParams(std::vector<double> const& valuesToFilter,
                      const std::vector<int>& sortedIndicesToExclude,
                      gsl::span<double> result);
   
   double
   getDefaultScalingFactor(amici::ParameterScaling scaling);
   
   double
   getDefaultOffsetParameter(amici::ParameterScaling scaling);
   
   double
   computeAnalyticalScalings(
       int scalingIdx,
       const std::vector<std::vector<double>>& modelOutputsUnscaled,
       const std::vector<std::vector<double>>& measurements,
       const AnalyticalParameterProvider& scalingReader,
       int numObservables);
   
   double
   computeAnalyticalOffsets(
       int offsetIdx,
       const std::vector<std::vector<double>>& modelOutputsUnscaled,
       const std::vector<std::vector<double>>& measurements,
       AnalyticalParameterProvider& offsetReader,
       int numObservables);
   
   double
   computeAnalyticalSigmas(
       int sigmaIdx,
       const std::vector<std::vector<double>>& modelOutputsScaled,
       const std::vector<std::vector<double>>& measurements,
       const AnalyticalParameterProvider& sigmaReader,
       int numObservables,
       double epsilonAbs = 1e-12,
       double epsilonRel = 0.01);
   
   void
   applyOptimalScaling(int scalingIdx,
                       double scalingLin,
                       std::vector<std::vector<double>>& modelOutputs,
                       AnalyticalParameterProvider const& scalingReader,
                       int numObservables);
   
   void
   applyOptimalOffset(int offsetIdx,
                      double offsetLin,
                      std::vector<std::vector<double>>& modelOutputs,
                      AnalyticalParameterProvider const& offsetReader,
                      int numObservables);
   
   std::vector<double>
   spliceParameters(const gsl::span<double const> reducedParameters,
                    const std::vector<int>& proportionalityFactorIndices,
                    const std::vector<int>& offsetParameterIndices,
                    const std::vector<int>& sigmaParameterIndices,
                    const std::vector<double>& scalingFactors,
                    const std::vector<double>& offsetParameters,
                    const std::vector<double>& sigmaParameters);
   
   template <typename T>
   std::vector<T>
   removeInnerParameters(const gsl::span<T const> allParameters,
                         const std::vector<int>& proportionalityFactorIndices,
                         const std::vector<int>& offsetParameterIndices,
                         const std::vector<int>& sigmaParameterIndices)
   {
       std::vector<T> outerParameters(
           allParameters.size() - proportionalityFactorIndices.size() -
           offsetParameterIndices.size() - sigmaParameterIndices.size());
   
       int nextOuterIdx = 0;
       for(int idxFull = 0; idxFull < static_cast<int>(allParameters.size());
            ++idxFull) {
   
           // skip if current parameter is scaling/offset/sigma
           if(std::find(proportionalityFactorIndices.begin(),
                         proportionalityFactorIndices.end(), idxFull)
               != std::end(proportionalityFactorIndices))
               continue;
   
           if(std::find(offsetParameterIndices.begin(),
                         offsetParameterIndices.end(), idxFull)
               != std::end(offsetParameterIndices))
               continue;
   
           if(std::find(sigmaParameterIndices.begin(),
                         sigmaParameterIndices.end(), idxFull)
               != std::end(sigmaParameterIndices))
               continue;
   
           // otherwise copy
           outerParameters[nextOuterIdx] = allParameters[idxFull];
           ++nextOuterIdx;
       }
   
       Ensures(nextOuterIdx == static_cast<int>(outerParameters.size()));
       return outerParameters;
   }
   
   std::vector<double>
   getOuterParameters(std::vector<double> const& fullParameters,
                      H5::H5File const& parameterFile,
                      std::string const& parameterPath);
   
   
   double
   computeNegLogLikelihood(
       std::vector<std::vector<double>> const& measurements,
       std::vector<std::vector<double>> const& modelOutputsScaled,
       const std::vector<std::vector<double>>& sigmas);
   
   double
   computeNegLogLikelihood(std::vector<double> const& measurements,
                           std::vector<double> const& modelOutputsScaled,
                           const std::vector<double>& sigmas);
   
   void
   checkGradientForAnalyticalParameters(std::vector<double> const& gradient,
                                        std::vector<int> const& analyticalIndices,
                                        double threshold);
   
   } // namespace parpe
   
   #endif // HIERARCHICALOPTIMIZATION_H