Program Listing for File hierarchicalOptimization.h

↰ Return to documentation for file (include/parpeamici/hierarchicalOptimization.h)

#ifndef HIERARCHICALOPTIMIZATION_H
#define HIERARCHICALOPTIMIZATION_H

#include <parpeamici/hierarchicalOptimizationAnalyticalParameterProvider.h>
#include <parpeamici/multiConditionProblem.h>
#include <parpeoptimization/optimizationProblem.h>

#include <gsl/gsl-lite.hpp>

#include <cmath>
#include <memory>
#include <numeric>


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;

    virtual void fillInitialParameters(gsl::span<double> buffer) const override;

    virtual void fillParametersMin(gsl::span<double> buffer) const override;

    virtual 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
    virtual 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
    virtual bool iterationFinished(
        gsl::span<const double> parameters,
        double objectiveFunctionValue,
        gsl::span<const double> objectiveFunctionGradient) const override;

    virtual 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