text:
/* roots/test.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <config.h>
#include <stdlib.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_roots.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_ieee_utils.h>

#include "roots.h"
#include "test.h"

/* stopping parameters */
const double EPSREL = (10 * GSL_DBL_EPSILON);
const double EPSABS = (10 * GSL_DBL_EPSILON);
const unsigned int MAX_ITERATIONS = 150;

void my_error_handler (const char *reason, const char *file,
                       int line, int err);

#define WITHIN_TOL(a, b, epsrel, epsabs) \
 ((fabs((a) - (b)) < (epsrel) * GSL_MIN(fabs(a), fabs(b)) + (epsabs)))

int
main (void)
{
  gsl_function F_sin, F_cos, F_func1, F_func2, F_func3, F_func4,
    F_func5, F_func6;

  gsl_function_fdf FDF_sin, FDF_cos, FDF_func1, FDF_func2, FDF_func3, FDF_func4,
    FDF_func5, FDF_func6;

  const gsl_root_fsolver_type * fsolver[4] ;
  const gsl_root_fdfsolver_type * fdfsolver[4] ;

  const gsl_root_fsolver_type ** T;
  const gsl_root_fdfsolver_type ** S;

  gsl_ieee_env_setup();

  fsolver[0] = gsl_root_fsolver_bisection;
  fsolver[1] = gsl_root_fsolver_brent;
  fsolver[2] = gsl_root_fsolver_falsepos;
  fsolver[3] = 0;

  fdfsolver[0] = gsl_root_fdfsolver_newton;
  fdfsolver[1] = gsl_root_fdfsolver_secant;
  fdfsolver[2] = gsl_root_fdfsolver_steffenson;
  fdfsolver[3] = 0;

  F_sin = create_function (sin_f) ;
  F_cos = create_function (cos_f) ;
  F_func1 = create_function (func1) ;
  F_func2 = create_function (func2) ;
  F_func3 = create_function (func3) ;
  F_func4 = create_function (func4) ;
  F_func5 = create_function (func5) ;
  F_func6 = create_function (func6) ;

  FDF_sin = create_fdf (sin_f, sin_df, sin_fdf) ;
  FDF_cos = create_fdf (cos_f, cos_df, cos_fdf) ;
  FDF_func1 = create_fdf (func1, func1_df, func1_fdf) ;
  FDF_func2 = create_fdf (func2, func2_df, func2_fdf) ;
  FDF_func3 = create_fdf (func3, func3_df, func3_fdf) ;
  FDF_func4 = create_fdf (func4, func4_df, func4_fdf) ;
  FDF_func5 = create_fdf (func5, func5_df, func5_fdf) ;
  FDF_func6 = create_fdf (func6, func6_df, func6_fdf) ;

  gsl_set_error_handler (&my_error_handler);

  for (T = fsolver ; *T != 0 ; T++)
    {
      test_f (*T, "sin(x) [3, 4]", &F_sin, 3.0, 4.0, M_PI);
      test_f (*T, "sin(x) [-4, -3]", &F_sin, -4.0, -3.0, -M_PI);
      test_f (*T, "sin(x) [-1/3, 1]", &F_sin, -1.0 / 3.0, 1.0, 0.0);
      test_f (*T, "cos(x) [0, 3]", &F_cos, 0.0, 3.0, M_PI / 2.0);
      test_f (*T, "cos(x) [-3, 0]", &F_cos, -3.0, 0.0, -M_PI / 2.0);
      test_f (*T, "x^20 - 1 [0.1, 2]", &F_func1, 0.1, 2.0, 1.0);
      test_f (*T, "sqrt(|x|)*sgn(x)", &F_func2, -1.0 / 3.0, 1.0, 0.0);
      test_f (*T, "x^2 - 1e-8 [0, 1]", &F_func3, 0.0, 1.0, sqrt (1e-8));
      test_f (*T, "x exp(-x) [-1/3, 2]", &F_func4, -1.0 / 3.0, 2.0, 0.0);
      test_f (*T, "(x - 1)^7 [0.9995, 1.0002]", &F_func6, 0.9995, 1.0002, 1.0);

      test_f_e (*T, "invalid range check [4, 0]", &F_sin, 4.0, 0.0, M_PI);
      test_f_e (*T, "invalid range check [1, 1]", &F_sin, 1.0, 1.0, M_PI);
      test_f_e (*T, "invalid range check [0.1, 0.2]", &F_sin, 0.1, 0.2, M_PI);
    }

  for (S = fdfsolver ; *S != 0 ; S++)
    {
      test_fdf (*S,"sin(x) {3.4}", &FDF_sin, 3.4, M_PI);
      test_fdf (*S,"sin(x) {-3.3}", &FDF_sin, -3.3, -M_PI);
      test_fdf (*S,"sin(x) {0.5}", &FDF_sin, 0.5, 0.0);
      test_fdf (*S,"cos(x) {0.6}", &FDF_cos, 0.6, M_PI / 2.0);
      test_fdf (*S,"cos(x) {-2.5}", &FDF_cos, -2.5, -M_PI / 2.0);
      test_fdf (*S,"x^{20} - 1 {0.9}", &FDF_func1, 0.9, 1.0);
      test_fdf (*S,"x^{20} - 1 {1.1}", &FDF_func1, 1.1, 1.0);
      test_fdf (*S,"sqrt(|x|)*sgn(x) {1.001}", &FDF_func2, 0.001, 0.0);
      test_fdf (*S,"x^2 - 1e-8 {1}", &FDF_func3, 1.0, sqrt (1e-8));
      test_fdf (*S,"x exp(-x) {-2}", &FDF_func4, -2.0, 0.0);
      test_fdf_e (*S,"max iterations x -> +Inf, x exp(-x) {2}", &FDF_func4, 2.0, 0.0);
      test_fdf_e (*S,"max iterations x -> -Inf, 1/(1 + exp(-x)) {0}", &FDF_func5, 0.0, 0.0);
    }

  test_fdf (gsl_root_fdfsolver_steffenson,
            "(x - 1)^7 {0.9}", &FDF_func6, 0.9, 1.0);

  /* now summarize the results */

  exit (gsl_test_summary ());
}

/* Using gsl_root_bisection, find the root of the function pointed to by f,
   using the interval [lower_bound, upper_bound]. Check if f succeeded and
   that it was accurate enough. */

void
test_f (const gsl_root_fsolver_type * T, const char * description, gsl_function *f,
        double lower_bound, double upper_bound, double correct_root)
{
  int status;
  size_t iterations = 0;
  double r, a, b;
  double x_lower, x_upper;
  gsl_root_fsolver * s;

  x_lower = lower_bound;
  x_upper = upper_bound;

  s = gsl_root_fsolver_alloc(T);
  gsl_root_fsolver_set(s, f, x_lower, x_upper) ;

  do
    {
      iterations++ ;

      gsl_root_fsolver_iterate (s);

      r = gsl_root_fsolver_root(s);

      a = gsl_root_fsolver_x_lower(s);
      b = gsl_root_fsolver_x_upper(s);

      if (a > b)
        gsl_test (GSL_FAILURE, "interval is invalid (%g,%g)", a, b);

      if (r < a || r > b)
        gsl_test (GSL_FAILURE, "r lies outside interval %g (%g,%g)", r, a, b);

      status = gsl_root_test_interval (a,b, EPSABS, EPSREL);
    }
  while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);

  gsl_test (status, "%s, %s (%g obs vs %g expected) ",
            gsl_root_fsolver_name(s), description,
            gsl_root_fsolver_root(s), correct_root);

  if (iterations == MAX_ITERATIONS)
    {
      gsl_test (GSL_FAILURE, "exceeded maximum number of iterations");
    }

  /* check the validity of the returned result */

  if (!WITHIN_TOL (r, correct_root, EPSREL, EPSABS))
    {
      gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)",
                r, correct_root);

    }

  gsl_root_fsolver_free(s);
}

void
test_f_e (const gsl_root_fsolver_type * T,
          const char * description, gsl_function *f,
          double lower_bound, double upper_bound, double correct_root)
{
  int status;
  size_t iterations = 0;
  double x_lower, x_upper;
  gsl_root_fsolver * s;

  x_lower = lower_bound;
  x_upper = upper_bound;

  s = gsl_root_fsolver_alloc(T);
  status = gsl_root_fsolver_set(s, f, x_lower, x_upper) ;

  gsl_test (status != GSL_EINVAL, "%s (set), %s", T->name, description);

  if (status == GSL_EINVAL)
    {
      gsl_root_fsolver_free(s);
      return ;
    }

  do
    {
      iterations++ ;
      gsl_root_fsolver_iterate (s);
      x_lower = gsl_root_fsolver_x_lower(s);
      x_upper = gsl_root_fsolver_x_lower(s);
      status = gsl_root_test_interval (x_lower, x_upper,
                                      EPSABS, EPSREL);
    }
  while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);

  gsl_test (!status, "%s, %s", gsl_root_fsolver_name(s), description,
            gsl_root_fsolver_root(s) - correct_root);

  gsl_root_fsolver_free(s);
}

void
test_fdf (const gsl_root_fdfsolver_type * T, const char * description,
        gsl_function_fdf *fdf, double root, double correct_root)
{
  int status;
  size_t iterations = 0;
  double prev = 0 ;

  gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T);
  gsl_root_fdfsolver_set (s, fdf, root) ;

  do
    {
      iterations++ ;
      prev = gsl_root_fdfsolver_root(s);
      gsl_root_fdfsolver_iterate (s);
      status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev,
                                   EPSABS, EPSREL);
    }
  while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);

  gsl_test (status, "%s, %s (%g obs vs %g expected) ",
            gsl_root_fdfsolver_name(s), description,
            gsl_root_fdfsolver_root(s), correct_root);

  if (iterations == MAX_ITERATIONS)
    {
      gsl_test (GSL_FAILURE, "exceeded maximum number of iterations");
    }

  /* check the validity of the returned result */

  if (!WITHIN_TOL (gsl_root_fdfsolver_root(s), correct_root,
                   EPSREL, EPSABS))
    {
      gsl_test (GSL_FAILURE, "incorrect precision (%g obs vs %g expected)",
                gsl_root_fdfsolver_root(s), correct_root);

    }
  gsl_root_fdfsolver_free(s);
}

void
test_fdf_e (const gsl_root_fdfsolver_type * T,
            const char * description, gsl_function_fdf *fdf,
            double root, double correct_root)
{
  int status;
  size_t iterations = 0;
  double prev = 0 ;

  gsl_root_fdfsolver * s = gsl_root_fdfsolver_alloc(T);
  status = gsl_root_fdfsolver_set (s, fdf, root) ;

  gsl_test (status, "%s (set), %s", T->name, description);

  do
    {
      iterations++ ;
      prev = gsl_root_fdfsolver_root(s);
      gsl_root_fdfsolver_iterate (s);
      status = gsl_root_test_delta(gsl_root_fdfsolver_root(s), prev,
                                   EPSABS, EPSREL);
    }
  while (status == GSL_CONTINUE && iterations < MAX_ITERATIONS);

  gsl_test (!status, "%s, %s", gsl_root_fdfsolver_name(s),
            description, gsl_root_fdfsolver_root(s) - correct_root);
  gsl_root_fdfsolver_free(s);
}

void
my_error_handler (const char *reason, const char *file, int line, int err)
{
  if (0)
    printf ("(caught [%s:%d: %s (%d)])\n", file, line, reason, err);
}

max_stars_repo_name: ruslankuzmin/julia


text:
// A cascade of IntegralHOGDetectors.
//
// Copyright 2012 Mark Desnoyer (mdesnoyer@gmail.com)

#ifndef __HOG_DETECTOR_INTEGRAL_HOG_CASCADE_TRAINER_H__
#define __HOG_DETECTOR_INTEGRAL_HOG_CASCADE_TRAINER_H__

#include <ros/ros.h>
#include <boost/shared_ptr.hpp>
#include <boost/scoped_ptr.hpp>
#include <functional>
#include <vector>
#include <string>
#include <queue>
#include <limits>
#include <cmath>
#include <opencv2/core/core.hpp>
#include <boost/unordered_map.hpp>
#include <boost/unordered_set.hpp>
#include <gsl/gsl_multimin.h>
#include "hog_detector/integral_hog_cascade.h"
#include "hog_detector/integral_hog_detector.h"
#include "hog_detector/interpolation.h"

namespace hog_detector {

// Class to calculate the various time components of the cascade
class IntegralCascadeTimeCalculator {
public:
  // Creates the calculator based on text files with timing data (in seconds)
  //
  // integralHistTime - single line with the time
  // fillBlockCacheTime - lines of <blockW>,<blockH>,<time>
  // svmEvalTime - lines of <descriptorSize>,<subWinW>,<subWinH>,<time/window>
  // trueHogTime - lines of <nWindows>,<time>
  IntegralCascadeTimeCalculator(const std::string& integralHistTime,
                                const std::string& fillBlockCacheTime,
                                const std::string& svmEvalTime,
                                const std::string& trueHogTime);

  // Time to compute the integral histogram
  double GetIntegralHistTime() const { return integralHistTime_; }

  // Time to fill the block cache for a given block size
  double GetFillCacheTime(const cv::Size& blockSize) const;

  // Time to do the svm evaluation of a single window
  double GetSVMEvalTime(int descriptorSize) const;

  // Time to compute the full HOG detector on N windows
  double GetTrueHogTime(int nWindows) const;

  int GetMaxWindows() const { return maxWindows_; }

  double GetMaxHogTime() const { return GetTrueHogTime(GetMaxWindows()); }

private:
  double integralHistTime_;

  boost::unordered_map<cv::Size, double> cacheFillTime_;

  boost::scoped_ptr<Interpolator> svmTime_;
  boost::scoped_ptr<Interpolator> trueHogTime_;

  int maxWindows_;
};

class IntegralHogCascadeTrainer {
public:
  // Constructor
  //
  // timeCalculator - takes ownership of it
  IntegralHogCascadeTrainer(IntegralCascadeTimeCalculator* timeCalculator,
                            const cv::Size& winSize);

  // Add a candidate IntegralHogDetector.
  void AddCandidateDetector(IntegralHogDetector::Ptr detector) {
    if (detector->winSize() == winSize_) {
      detectorIdx_[detector.get()] = detectors_.size();
      detectors_.push_back(detector);
    } else {
      ROS_WARN_STREAM("The detector has the wrong window size: ("
                      << winSize_.width << ","
                      << winSize_.height  << ") vs. ("
                      << detector->winSize().width << ","
                      << detector->winSize().height << ")");
    }
  }

  // Trains the cascade detector
  //
  // Inputs:
  // filenames - Filenames of images to train with. Uses the center
  //             winSize window in the image.
  // labels - Labels for all the images. 1 is positive, -1 is negative
  // budget - Processing budget in seconds (average per/frame)
  // missCost - The cost of a miss
  // falsePosCost - The cost of a false positive
  // fracNegSampled - fraction of the number of negative examples we're
  //                  actually training with. So the number of negative
  //                  samples seen by this class is trueNegs/fracNegSampled
  void Train(const std::vector<std::string>& filenames,
             const std::vector<float>& labels,
             float budget,
             float missCost,
             float falsePosCost);
  // Train using the supermodular-submodular approach
  void TrainSupSub(const std::vector<std::string>& filenames,
                   const std::vector<float>& label,
                   float missCost,
                   float falsePosCost,
                   float timeCostPerError, // time is per seconds
                   float fracNegSampled);

  // Training can also be done by adding one image at a time and then
  // calling TrainWithLoadedData
  void AddImageForTraining(const cv::Mat& image,
                           const std::vector<cv::Rect>& rois,
                           float label);
  void TrainWithLoadedData(float budget,
                           float missCost,
                           float falsePosCost);
  void TrainSupSubWithLoadedData(float missCost,
                                 float falsePosCost,
                                 float timeCostPerError,
                                 float fracNegSampled);

  const IntegralHogCascade* cascade() const { return cascade_.get(); }

private:
  // Structure to hold the score and selected threshold for a given detector
  struct DetectorStats {
    DetectorStats(float _score, float _thresh, unsigned int _detectorIdx)
      : score(_score), thresh(_thresh), detectorIdx(_detectorIdx) {}

    // Sort primarily based on the score
    bool operator<(const DetectorStats& other) const {
      if (fequal(score, other.score)) {
        if (detectorIdx == other.detectorIdx) {
          return thresh < other.thresh;
        }
        return detectorIdx < other.detectorIdx;
      }
      return score < other.score;
    }
    bool operator==(const DetectorStats& other) const {
      return detectorIdx == other.detectorIdx && fequal(score, other.score)
        && fequal(thresh, other.thresh);
    }
    bool operator>(const DetectorStats& other) const {return other < *this;}

    bool fequal(float a, float b,
                float ep=std::numeric_limits<float>::epsilon()*100) const {
      return std::fabs(a-b) < ep;
    }

    float score;
    float thresh;
    unsigned int detectorIdx;
  };

  // Number of images being trained on
  int nImages_;

  cv::Size winSize_;

  boost::scoped_ptr<IntegralCascadeTimeCalculator> timeCalculator_;

  // All the possible detectors
  std::vector<IntegralHogDetector::Ptr> detectors_;

  // Scores of the training images for all the detectors. Each row is
  // a different detector and each column is a different training
  // image.
  cv::Mat_<float> scores_;

  // Labels for all the loaded images
  cv::Mat_<float> labels_;

  // Inverse index from the detector to its index in scores_ and labels_
  boost::unordered_map<IntegralHogDetector*, int> detectorIdx_;

  // The trained cascade
  boost::scoped_ptr<IntegralHogCascade> cascade_;

  // The true sum of the number of images given that the negatives are
  // sampled more sparsely.
  double trueImageSum_;

  void LoadImagesForTraining(const std::vector<std::string>& filenames,
                             const std::vector<float>& labels);

  // Trains a cascade using a greedy metric.
  //
  // Inputs:
  // budget - Time budget in seconds,
  // mistCost - Cost of a miss
  // falsePosCost - Cost of a false positive
  // metric - The metric to use
  //
  // Outputs:
  // cascade - The trained cascade
  // returns - Total cost of the cascade
  template <typename T>
  float TrainUsingGreedyMetric(float budget, float missCost,
                               float falsePosCost,
                               T metric,
                               IntegralHogCascade* cascade);

  // Two different metrics for selecting the next best detector.
  //
  // Inputs:
  // missCost - Cost of a miss
  // falsePosCost - Cost of a false positive
  // inputWindows - - Boolean matrix identifying which training windows are
  //                still in play. Will be updated
  // allocatedTime - The time already allocated. Will be updated after this
  //                 selection
  // cacheSizesUsed - Set of cache sizes used. Will be updated after this
  //                  selection.
  //
  // Returns:
  // pair of threshold and index of the selected detector
  std::pair<float, int> SelectBestDetectorIgnoringTime(
    float missCost,
    float falsePosCost,
    cv::Mat& inputWindows,
    float& allocatedTime,
    boost::unordered_set<cv::Size>& cacheSizesUsed);

  std::pair<float, int> SelectBestDetectorPerTime(
    float missCost,
    float falsePosCost,
    cv::Mat& inputWindows,
    float& allocatedTime,
    boost::unordered_set<cv::Size>& cacheSizesUsed);

  // Updates all the tracking data for when a detector is chosen
  void ChooseDetector(int detectorI, float thresh,
                      cv::Mat& inputWindows,
                      float& allocatedTime,
                      boost::unordered_set<cv::Size>& cacheSizesUsed);

  // Selects the best detector that is under budget
  void SelectBestDetectorUnderBudget(
    float budget, float missCost,
    float falsePosCost,
    const cv::Mat& inputWindows,
    float allocatedTime,
    const boost::unordered_set<cv::Size>& cacheSizesUsed,
    float* thresh,
    int* chosenDetector,
    float* cost);

  // Selects the best detector according to an approximate submodular measure that is tight at cascade.
  //
  // Inputs:
  //   missCost - Cost of a miss
  //   falsePosCost - Cost of a false positive
  //   timeCostPerError - Conversion factor for 1s of processing time compared
  //                    to an error.
  //   cascade - The best cascade so far. The approximation is tight at
  //               this point.
  //   winPastCascade - Boolean vector specifying the windows that pass the cascade
  //
  // Outputs:
  //  thresh - Threshold of the detector to choose
  //  chosenDetector - Index of the chosen detector
  //  returns - The delta cost for the added detector
  typedef std::priority_queue<DetectorStats, std::vector<DetectorStats>,
                              std::greater<DetectorStats> > DetectorQueue;
  float SelectBestDetectorApprox(
    float missCost,
    float falsePosCost,
    float timeCostPerError,
    float fracNegSampled,
    const IntegralHogCascade& cascade,
    const cv::Mat& winPastApproxCascade,
    DetectorQueue* detectorQueue,
    float* thresh,
    int* chosenDetector);

  // Compute the bayes risk for a given threshold on a set of images
  // whose scores are curScores.
  double ComputeBayesRisk(const cv::Mat& inputWindows,
                          const cv::Mat_<float>& curScores,
                          float thresh,
                          double missCost,
                          double falsePosCost);

  double ComputeProcessingTimeForDetector(
    float fracInputWindows,
    int detectorI,
    const boost::unordered_set<cv::Size>& cacheSizesUsed);

  // For the SupSub proceedure find the cascade that minimizes the
  // approximate cost. The cost is approximate because the D_a and C_M
  // terms are approximated by an upper bound.
  //
  // Inputs:
  // missCost - Cost of a miss
  // falsePosCost - Cost of a false positive, already adjusted with
  //                fracNegSampled
  // timeCostPerError - Conversion factor for 1s of processing time compared
  //                    to an error.
  // bestCascade - The best cascade so far. The approximation is tight at
  //               this point.
  // fracNegSampled - fraction of the number of negative examples we're
  //                  actually training with. So the number of negative
  //                  samples seen by this class is trueNegs/fracNegSampled
  //
  // Outputs:
  // bestCascade - An update
  // Returns - The true cost of the best cascade (not the approximate cost)
  float MinimizeApproxCost(float missCost, float falsePosCost,
                           float timeCostPerError,
                           float fracNegSampled,
                           boost::scoped_ptr<IntegralHogCascade>* ioCascade);

  // Computes the total cost of the cascade
  double ComputeTotalCost(float missCost,
                          float falsePosCost,
                          float timeCostPerError,
                          float fracNegSampled,
                          const IntegralHogCascade& cascade);

  // Initializes the data for computing the approximate submodular cost.
  //
  // Inputs:
  // missCost - Cost of a miss
  // falsePosCost - Cost of a false positive
  // timeCostPerError - Conversion factor for 1s of processing time compared
  //                    to an error.
  // fracNegSampled - fraction of the number of negative examples we're
  //                  actually training with. So the number of negative
  //                  samples seen by this class is trueNegs/fracNegSampled
  // cascade - The cascade to approximate to
  //
  // Outputs:
  // removeCosts - Map for each stage for the change in cost if it is removed
  // pastCascade - Vector specifying which entries got past the cascade from
  //               rows in scores_
  void InitializeApproxSubCost(
    float missCost,
    float falsePosCost,
    float timeCostPerError,
    float fracNegSampled,
    const IntegralHogCascade& cascade,
    boost::unordered_map<IntegralHogCascade::Stage, float>* removeCosts,
    cv::Mat* pastCascade);

  float CalculateDeltaForRemovingStage(
    float falsePosCost,
    float timeCostPerError,
    float fracNegSampled,
    const IntegralHogCascade& cascade,
    const boost::unordered_map<IntegralHogCascade::Stage, float>& approxCosts,
    const IntegralHogCascade::Stage& stage);

  // Find the best approximate cost for a given detector.
  //
  // missCost, falsePosCost, timeCostPerError, fraceNegSampled -
  //   relative weights for the costs.
  // scores - score of the detector evaluated on each candidate region
  // usefulPassApproxCascade - usefule regions that pass the approx cascade
  // falsePosWins - False positive windows from the current cascade
  // windowsPassed - Windows that pass the current cascade
  // usefulWindows - Windows that are positives for training
  // notUsefuleWindows - Windows that are negatives for training
  //
  // Returns: Pair of <cost, bestThreshold>
  std::pair<float, float> GetBestApproxCost(
    unsigned int detectorIdx,
    float missCost,
    float falsePosCost,
    float timeCostPerError,
    float fracNegSampled,
    double baseHogTime,
    float nPassedApproxCascade,
    const cv::Mat& usefulPassApproxCascade,
    const cv::Mat& falsePosWins,
    const cv::Mat& windowsPassed,
    const cv::Mat& usefulWindows,
    const cv::Mat& notUsefulWindows);

};

// Class that finds the threshold to minimize the cost of a filter
class DeltaMinimizer {
public:
  DeltaMinimizer();
  ~DeltaMinimizer();

  // Minimizes the cost of adding a current detector
  //
  // Inputs:
  // missCost, falsePosCost, timeCostPerError, fraceNegSampled -
  //   relative weights for the costs.
  // scores - score of the detector evaluated on each candidate region
  // usefulPassApproxCascade - usefule regions that pass the approx cascade
  // falsePosWins - False positive windows from the current cascade
  // windowsPassed - Windows that pass the current cascade
  // usefulWindows - Windows that are positives for training
  // notUsefuleWindows - Windows that are negatives for training
  // timeCalculator - To calculate the time
  // cpuDelta - cpu cost to run this detector
  float MinimizeCost(float missCost,
                     float falsePosCost,
                     float timeCostPerError,
                     float fracNegSampled,
                     int nImages,
                     float trueImageSum,
                     double baseHogTime,
                     const cv::Mat_<float>& scores,
                     const cv::Mat& usefulPassApproxCascade,
                     const cv::Mat& falsePosWins,
                     const cv::Mat& windowsPassed,
                     const cv::Mat& usefulWindows,
                     const cv::Mat& notUsefulWindows,
                     const IntegralCascadeTimeCalculator* timeCalculator,
                     float cpuDelta);

  float bestCost() const { return bestCost_; }
  float bestThresh() const { return bestThresh_; }

  // Evaluates the cost function and its derivative with respect to
  // the number of entries filtered. Public because gsl is old and
  // needs acess to this.
  void EvalWithDeriv(const gsl_vector* nFiltered, void* params,
                     double* val, gsl_vector* derivs) const;
  double EvalCost(const gsl_vector* nFiltered, void* params) const;
  void EvalDeriv(const gsl_vector* nFiltered, void* params,
                 gsl_vector* derivs) const;

private:
  float bestCost_;
  float bestThresh_;

  // Used for the minimization
  gsl_multimin_fdfminimizer* minimizer_;
  gsl_multimin_function_fdf minFunc_;

  // Data used for the calculator
  double cpuDelta_;
  // Map from the number of new windows thresholded to the threshold.
  boost::scoped_ptr<Interpolator> thresholdInterpolator_;
  // Interpolator for the cost of extra misses with a given threshold
  boost::scoped_ptr<Interpolator> missCostInterpolator_;
  // Interpolator for the cost of the false positives relative to the
  // number of new windows thresholded.
  boost::scoped_ptr<Interpolator> falsePosInterpolator_;
  // Interpolator for the HOG time delta
  boost::scoped_ptr<Interpolator> hogTimeInterpolator_;

  void RunMinimization(double minWindowsThreshed, double maxWindowsThreshed,
                       double startPoint, double firstStep,
                       float* bestCost, float* bestThresh);

};

} // namespace

#endif //__HOG_DETECTOR_INTEGRAL_HOG_CASCADE_TRAINER_H__

ext: h


text:
/***************************************************************************************
 *  Multivariate Normal density function and random number generator
 *  Multivariate Student t density function and random number generator
 *  Wishart random number generator
 *  Using GSL -> www.gnu.org/software/gsl
 *
 *  Copyright (C) 2006  Ralph dos Santos Silva
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 *
 *  AUTHOR
 *     Ralph dos Santos Silva,  [EMAIL PROTECTED]
 *     March, 2006
***************************************************************************************/
#include <math.h>
#include <gsl/gsl_sf_gamma.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include "../hmm/adkGSL.h"
#include "mvnorm_gsl.h"
/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rmvnorm(const gsl_rng *r, const int n, const gsl_vector *mean, const gsl_matrix *var, gsl_vector *result){
/* multivariate normal distribution random number generator */
/*
*	n	dimension of the random vetor
*	mean	vector of means of size n
*	var	variance matrix of dimension n x n
*	result	output variable with a sigle random vector normal distribution generation
*/
int k;
gsl_matrix *work = gsl_matrix_alloc(n,n);

gsl_matrix_memcpy(work,var);
gsl_linalg_cholesky_decomp(work);

for(k=0; k<n; k++)
	gsl_vector_set( result, k, gsl_ran_ugaussian(r) );

gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, work, result);
gsl_vector_add(result,mean);

gsl_matrix_free(work);

return 0;
}

/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rmvnorm_prealloc(const gsl_rng *r, const int n, const gsl_vector *mean, const gsl_matrix *var, gsl_matrix *work, gsl_vector *result){
/* multivariate normal distribution random number generator */
/*
*	n	dimension of the random vetor
*	mean	vector of means of size n
*	var	variance matrix of dimension n x n
*	result	output variable with a sigle random vector normal distribution generation
*/
int k;

gsl_matrix_memcpy(work,var);
gsl_linalg_cholesky_decomp(work);

for(k=0; k<n; k++)
	gsl_vector_set( result, k, gsl_ran_ugaussian(r) );

gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, work, result);
gsl_vector_add(result,mean);

return 0;
}

/*****************************************************************************************************************/
/*****************************************************************************************************************/
double dmvnorm(const int n, const gsl_vector *x, const gsl_vector *mean, const gsl_matrix *var){
/* multivariate normal density function    */
/*
*	n	dimension of the random vetor
*	mean	vector of means of size n
*	var	variance matrix of dimension n x n
*/
int s;
double ax,ay;
gsl_vector *ym, *xm;
gsl_matrix *work = gsl_matrix_alloc(n,n),
           *winv = gsl_matrix_alloc(n,n);
gsl_permutation *p = gsl_permutation_alloc(n);

gsl_matrix_memcpy( work, var );
gsl_linalg_LU_decomp( work, p, &s );
gsl_linalg_LU_invert( work, p, winv );
ax = gsl_linalg_LU_det( work, s );
gsl_matrix_free( work );
gsl_permutation_free( p );

xm = gsl_vector_alloc(n);
gsl_vector_memcpy( xm, x);
gsl_vector_sub( xm, mean );
ym = gsl_vector_alloc(n);
gsl_blas_dsymv(CblasUpper,1.0,winv,xm,0.0,ym);
gsl_matrix_free( winv );
gsl_blas_ddot( xm, ym, &ay);
gsl_vector_free(xm);
gsl_vector_free(ym);
ay = exp(-0.5*ay)/sqrt( pow((2*M_PI),n)*ax );

return ay;
}

/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rmvt(const gsl_rng *r, const int n, const gsl_vector *location, const gsl_matrix *scale, const int dof, gsl_vector *result){
/* multivariate Student t distribution random number generator */
/*
*	n	 dimension of the random vetor
*	location vector of locations of size n
*	scale	 scale matrix of dimension n x n
*	dof	 degrees of freedom
*	result	 output variable with a single random vector normal distribution generation
*/
int k;
gsl_matrix *work = gsl_matrix_alloc(n,n);
double ax = 0.5*dof;

ax = gsl_ran_gamma(r,ax,(1/ax));     /* gamma distribution */

gsl_matrix_memcpy(work,scale);
gsl_matrix_scale(work,(1/ax));       /* scaling the matrix */
gsl_linalg_cholesky_decomp(work);

for(k=0; k<n; k++)
	gsl_vector_set( result, k, gsl_ran_ugaussian(r) );

gsl_blas_dtrmv(CblasLower, CblasNoTrans, CblasNonUnit, work, result);
gsl_vector_add(result, location);

gsl_matrix_free(work);

return 0;
}

/*****************************************************************************************************************/
/*****************************************************************************************************************/
double dmvt(const int n, const gsl_vector *x, const gsl_vector *location, const gsl_matrix *scale, const int dof){
/* multivariate Student t density function */
/*
*	n	 dimension of the random vetor
*	location vector of locations of size n
*	scale	 scale matrix of dimension n x n
*	dof	 degrees of freedom
*/
int s;
double ax,ay,az=0.5*(dof + n);
gsl_vector *ym, *xm;
gsl_matrix *work = gsl_matrix_alloc(n,n),
           *winv = gsl_matrix_alloc(n,n);
gsl_permutation *p = gsl_permutation_alloc(n);

gsl_matrix_memcpy( work, scale );
gsl_linalg_LU_decomp( work, p, &s );
gsl_linalg_LU_invert( work, p, winv );
ax = gsl_linalg_LU_det( work, s );
gsl_matrix_free( work );
gsl_permutation_free( p );

xm = gsl_vector_alloc(n);
gsl_vector_memcpy( xm, x);
gsl_vector_sub( xm, location );
ym = gsl_vector_alloc(n);
gsl_blas_dsymv(CblasUpper,1.0,winv,xm,0.0,ym);
gsl_matrix_free( winv );
gsl_blas_ddot( xm, ym, &ay);
gsl_vector_free(xm);
gsl_vector_free(ym);

ay = pow((1+ay/dof),-az)*gsl_sf_gamma(az)/(gsl_sf_gamma(0.5*dof)*sqrt( pow((dof*M_PI),n)*ax ));

return ay;
}
/*****************************************************************************************************************/
/*****************************************************************************************************************/
int rwishart(const gsl_rng *r, const int n, const int dof, const gsl_matrix *scale, gsl_matrix *work, gsl_matrix *output){
/* Wishart distribution random number generator */
/*
*	n	 gives the dimension of the random matrix
*	dof	 degrees of freedom
*	scale	 scale matrix of dimension n x n
*	result	 output variable with a single random matrix Wishart distribution generation
*/
int k,l;

for(k=0; k<scale->size1; k++){
	gsl_matrix_set( work, k, k, sqrt( gsl_ran_chisq( r, (dof-k) ) ) );
	for(l=0; l<k; l++){
		gsl_matrix_set( work, k, l, gsl_ran_ugaussian(r) );
	}
}
gsl_matrix_memcpy(output,scale);
gsl_linalg_cholesky_decomp(output);
gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasNonUnit,1.0,output,work);
gsl_blas_dsyrk(CblasUpper,CblasNoTrans,1.0,work,0.0,output);
return 0;
}

// ADK additions from here on

/*****************************************************************************************************************/
/*****************************************************************************************************************/
double dmvnorm_prealloc(const int n, const gsl_vector *x, const gsl_vector *mean, const gsl_matrix *var, gsl_matrix *work,
	 gsl_matrix *winv,  gsl_vector *ym,  gsl_vector *xm, gsl_permutation *p){
/* multivariate normal density function    */
/*
*	n	dimension of the random vetor
*	mean	vector of means of size n
*	var	variance matrix of dimension n x n
*
*   added this prealloc version to cut down on memory calls
*/

int s;
double ax,ay;

gsl_matrix_memcpy( work, var );
gsl_linalg_LU_decomp( work, p, &s );
gsl_linalg_LU_invert( work, p, winv );
ax = gsl_linalg_LU_det( work, s );
gsl_vector_memcpy( xm, x);
gsl_vector_sub( xm, mean );
gsl_blas_dsymv(CblasUpper,1.0,winv,xm,0.0,ym);
gsl_blas_ddot( xm, ym, &ay);
ay = exp(-0.5*ay)/sqrt( pow((2*M_PI),n)*ax );

return ay;
}

/*****************************************************************************************************************/
/*****************************************************************************************************************/
double dInverseWishart(const gsl_matrix *x, const gsl_matrix *sigma, const double dof, gsl_matrix *work,
	 gsl_matrix *winv,  gsl_permutation *p){

	double denom, num;
	int s, d, i;
	double detSig,detX,trace;

	//get det of both matrices, take inverse of X
	gsl_matrix_memcpy( work, sigma );
	gsl_linalg_LU_decomp( work, p, &s );
	detSig = gsl_linalg_LU_det( work, s );
	gsl_matrix_memcpy( work, x );
	gsl_linalg_LU_decomp( work, p, &s );
	detX = gsl_linalg_LU_det( work, s );
	gsl_linalg_LU_invert( work, p, winv );
	//get product of gamma's for denominator; calc denom
	d = x->size1;
	denom = 1.0;
	denom *= pow(2,(d*dof/2.0)) * mv_gamma_func(dof/2.0, d);

	//multiply x by inverse of sigma and take trace
	gsl_blas_dgemm(CblasNoTrans, CblasNoTrans,1.0, sigma, winv,0.0, work);
	trace = gsl_matrix_trace(work);
	num =  pow(detSig, dof / 2.0) * pow(detX , (-dof-d-1)/ 2.0) * exp(-0.5 * trace);
	return(num/denom);
}

/*****************************************************************************************************************/
//rwishartGelman-- uses alg described in Gelman et al appendix
void rwishartGelman(const gsl_rng *r, const int n, const int dof, const gsl_matrix *scale,
	gsl_matrix *work,gsl_matrix *work2, gsl_vector *mean,gsl_vector *xm, gsl_matrix *output){
/* Wishart distribution random number generator */
/*
*	n	 gives the dimension of the random matrix
*	dof	 degrees of freedom
*	scale	 scale matrix of dimension n x n
*	result	 output variable with a single random matrix Wishart distribution generation
*/
	int i;

	gsl_vector_set_all(mean,0.0);
	gsl_matrix_set_all(output,0.0);
	for(i = 0; i < dof; i++){
		gsl_vector_set_all(xm,0);
		rmvnorm_prealloc(r, n, mean, scale, work, xm);
		gsl_vector_outer_product(xm,xm,work2);
		gsl_matrix_add(output,work2);
	}
}

//rwishartOdellFeiveson-- uses alg of Odell and Feiveson
void rwishartOdellFeiveson(const gsl_rng *r, const int n, const int dof, const gsl_matrix *scale,
	gsl_matrix *work,gsl_matrix *work2, gsl_vector *mean,gsl_vector *xm, gsl_matrix *output){
/* Wishart distribution random number generator */
/*
*	n	 gives the dimension of the random matrix
*	dof	 degrees of freedom
*	scale	 scale matrix of dimension n x n
*	result	 output variable with a single random matrix Wishart distribution generation
*/
	int i,j, k;
	double sum, tmp;

	gsl_vector_set_all(mean,0.0);
	gsl_matrix_set_all(output,0.0);

	for(j=0;j<n;j++){
		gsl_vector_set(mean,j,gsl_ran_chisq(r,dof-j));
		for(i=0;i<j;i++){
			gsl_matrix_set(work,i,j,gsl_ran_ugaussian(r));
		}
	}

	//set up B matrix in work2
	gsl_matrix_set(work2,0,0,gsl_vector_get(mean,0));
	for(j=0;j<n;j++){
		for(i=0;i<j;i++){
			if(i==0){
				tmp = sqrt(gsl_vector_get(mean,j)) * gsl_matrix_get(work,0,j);
				gsl_matrix_set(work2,0,j,tmp);
			}
			else{
				tmp = sqrt(gsl_vector_get(mean,i)) * gsl_matrix_get(work,i,j);
				for(k=0;k<i;k++){
					tmp += gsl_matrix_get(work,k,i) * gsl_matrix_get(work,k,j);
				}
				gsl_matrix_set(work2,i,j,tmp);
			}
		}
		tmp = gsl_vector_get(mean,j);
		for(k=0;k<j;k++){
			tmp += gsl_matrix_get(work,k,j) * gsl_matrix_get(work,k,j);
		}
		gsl_matrix_set(work2,j,j,tmp);
	}
	gsl_matrix_memcpy(output,scale);
	gsl_linalg_cholesky_decomp(output);
	gsl_blas_dtrmm(CblasLeft,CblasLower,CblasNoTrans,CblasNonUnit,1.0,output,work2);
	gsl_blas_dtrmm(CblasRight,CblasUpper,CblasNoTrans,CblasNonUnit,1.0,work2,output);

}

//mv_gamma_func -- returns the multivariate gamma function
double mv_gamma_func(double a, double d){
	double val = 1.0;
	int i;

	for(i=1;i<=d;i++){
		 val *= gsl_sf_gamma(a- (0.5 * (i-1)));
		}
	val *=  pow(M_PI,(d*(d-1)/4.0));
	return(val);
}

ext: c

max_stars_repo_name: andrewkern/pgLib

max_issues_repo_name: andrewkern/pgLib

max_forks_repo_name: andrewkern/pgLib


text:
/***************************************************************************
    File                 : Fit.h
    Project              : QtiPlot
    --------------------------------------------------------------------
    Copyright            : (C) 2006 by Ion Vasilief
    Email (use @ for *)  : ion_vasilief*yahoo.fr
    Description          : Fit base class

 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *  This program is free software; you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation; either version 2 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the Free Software           *
 *   Foundation, Inc., 51 Franklin Street, Fifth Floor,                    *
 *   Boston, MA  02110-1301  USA                                           *
 *                                                                         *
 ***************************************************************************/
#ifndef FIT_H
#define FIT_H

#include <QObject>

#include "../ApplicationWindow.h"
#include "Filter.h"

#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_multimin.h>

class Table;
class Matrix;

//! Fit base class
class Fit : public Filter
{
	Q_OBJECT

	public:

		typedef double (*fit_function_simplex)(const gsl_vector *, void *);
		typedef int (*fit_function)(const gsl_vector *, void *, gsl_vector *);
		typedef int (*fit_function_df)(const gsl_vector *, void *, gsl_matrix *);
		typedef int (*fit_function_fdf)(const gsl_vector *, void *, gsl_vector *, gsl_matrix *);

		enum Algorithm{ScaledLevenbergMarquardt, UnscaledLevenbergMarquardt, NelderMeadSimplex};
		enum WeightingMethod{NoWeighting, Instrumental, Statistical, Dataset};
        enum FitType{BuiltIn = 0, Plugin = 1, User = 2};

		Fit(ApplicationWindow *parent, Graph *g = 0, const QString& name = QString());
		Fit(ApplicationWindow *parent, Table *t, const QString& name = QString());
		~Fit();

		//! Actually does the fit. Should be reimplemented in derived classes.
		virtual void fit();
        virtual bool run(){fit(); return true;};

		//! Sets the data set to be used for weighting
		bool setWeightingData(WeightingMethod w, const QString& colName = QString::null);

		void setDataCurve(int curve, double start, double end);
		bool setDataFromTable(Table *t, const QString& xColName, const QString& yColName, int from = 1, int to = -1);

		QString resultFormula(){return d_result_formula;};
		QString formula(){return d_formula;};
		virtual void setFormula(const QString&){};

		int numParameters(){return d_p;};
		QStringList parameterNames(){return d_param_names;};
		virtual void setParametersList(const QStringList&){};
        void setParameterExplanations(const QStringList& lst){d_param_explain = lst;};

        double initialGuess(int parIndex){return gsl_vector_get(d_param_init, parIndex);};
		void setInitialGuess(int parIndex, double val){gsl_vector_set(d_param_init, parIndex, val);};
		void setInitialGuesses(double *x_init);

		virtual void guessInitialValues(){};

		void setParameterRange(int parIndex, double left, double right);
		void setAlgorithm(Algorithm s){d_solver = s;};

		//! Specifies weather the result of the fit is a function curve
		void generateFunction(bool yes, int points = 100);

		//! Output string added to the plot as a new legend
		virtual QString legendInfo();

		//! Returns a vector with the fit results
		double* results(){return d_results;};

		//! Returns a vector with the fit residuals
		double* residuals();

		//! Plot residuals and display data values in a column
		QwtPlotCurve* showResiduals();

		void showPredictionLimits(double confidenceLevel);
		void showConfidenceLimits(double confidenceLevel);
		//! Lower Confidence Limit
		double lcl(int parIndex, double confidenceLevel);
		//! Upper Confidence Limit
		double ucl(int parIndex, double confidenceLevel);

		//! Returns a vector with the standard deviations of the results
		double* errors();

		//! Returns the sum of squares of the residuals from the best-fit line
		double chiSquare() {return chi_2;};

		//! Returns R^2
		double rSquare();

		//! Returns adjusted R^2
		double adjustedRSquare(){return d_adjusted_r_square;};

		//! Returns the Residual Sum of Squares
		double rss(){return d_rss;};

		//! Returns the Root Mean Squared Error
		double rmse(){return sqrt(d_rss/(d_n - d_p));};

		//! Specifies wheather the errors must be scaled with sqrt(chi_2/dof)
		void scaleErrors(bool yes = true){d_scale_errors = yes;};

		Table* parametersTable(const QString& tableName);
		void writeParametersToTable(Table *t, bool append = false);

		Matrix* covarianceMatrix(const QString& matrixName);

        bool save(const QString& fileName);
        bool load(const QString& fileName);

        FitType type(){return d_fit_type;};
        void setType(FitType t){d_fit_type = t;};

        QString fileName(){return d_file_name;};
		void setFileName(const QString& fn){d_file_name = fn;};

        //! Frees the memory allocated for the X and Y data sets
        void freeMemory();

        //! Calculates the data for the output fit curve
        virtual double eval(double *, double){return 0.0;};

	private:
		void init();

		//! Pointer to the GSL multifit minimizer (for simplex algorithm)
		gsl_multimin_fminimizer * fitSimplex(gsl_multimin_function f, int &iterations, int &status);

		//! Pointer to the GSL multifit solver
		gsl_multifit_fdfsolver * fitGSL(gsl_multifit_function_fdf f, int &iterations, int &status);

		//! Customs and stores the fit results according to the derived class specifications. Used by exponential fits.
		virtual void customizeFitResults(){};

	protected:
		//! Allocates the memory for the fit workspace
		void initWorkspace(int par);
		//! Frees the memory allocated for the fit workspace
		void freeWorkspace();
		//! Adds the result curve as a FunctionCurve to the plot, if d_gen_function = true
		virtual FunctionCurve * insertFitFunctionCurve(const QString& name, double *x, double *y, int penWidth = 1);

		//! Adds the result curve to the plot
		virtual void generateFitCurve();

        //! Calculates the data for the output fit curve and store itin the X an Y vectors
		virtual void calculateFitCurveData(double *X, double *Y) {Q_UNUSED(X) Q_UNUSED(Y)};

		//! Output string added to the result log
		virtual QString logFitInfo(int iterations, int status);

		fit_function d_f;
		fit_function_df d_df;
		fit_function_fdf d_fdf;
		fit_function_simplex d_fsimplex;

		//! Number of fit parameters
		int d_p;

		//! Initial guesses for the fit parameters
		gsl_vector *d_param_init;

		/*! \brief Tells whether the fitter uses non-linear/simplex fitting
		 * with an initial parameters set, that must be freed in the destructor.
		 */
		bool is_non_linear;

		//! weighting data set used for the fit
		double *d_w;

		//! Names of the fit parameters
		QStringList d_param_names;

		//! Stores a list of short explanations for the significance of the fit parameters
		QStringList d_param_explain;

		//! Specifies weather the result curve is a FunctionCurve or a normal curve with the same x values as the fit data
		bool d_gen_function;

		//! Algorithm type
		Algorithm d_solver;

		//! The fit formula given on input
		QString d_formula;

		//! The result fit formula, where the fit parameters are replaced with the calculated values.
		QString d_result_formula;

		//! Covariance matrix
		gsl_matrix *covar;

		//! The kind of weighting to be performed on the data
		WeightingMethod d_weighting;

		//! The name of the weighting dataset
		QString weighting_dataset;

		//! Stores the result parameters
		double *d_results;

		//! Stores standard deviations of the result parameters
		double *d_errors;

		//! Stores fit residuals
		double *d_residuals;

		//! The sum of squares of the residuals from the best-fit line
		double chi_2;

		//! Residual sum of squares
		double d_rss;

		//! Adjusted R^2
		double d_adjusted_r_square;

		//! Specifies wheather the errors must be scaled with sqrt(chi_2/dof)
		bool d_scale_errors;

		//! Table window used for the output of fit parameters
		Table *d_param_table;

		//! Matrix window used for the output of covariance matrix
		Matrix *d_cov_matrix;

		FitType d_fit_type;

		//! Path of the XML file where the user stores the fit model
        QString d_file_name;

		//! Stores the left limits of the research interval for the result parameters
		double *d_param_range_left;

		//! Stores the right limits of the research interval for the result parameters
		double *d_param_range_right;
};

#endif

max_stars_repo_name: hoehnp/SpaceDesignTool

max_issues_repo_name: hoehnp/SpaceDesignTool

max_forks_repo_name: hoehnp/SpaceDesignTool


text:
#pragma once

#include <type_traits>
#include <gsl/gsl>

#include "sly/macros.h"
#include "sly/types.h"

template <class T>
struct is_unique_ptr : std::false_type
{};

template <class T, class D>
struct is_unique_ptr<std::unique_ptr<T, D>> : std::true_type
{};

template <class T>
struct unique_ptr_type
{};

template <class T, class D>
struct unique_ptr_type<std::unique_ptr<T, D>>
{
    using type = std::unique_ptr<T, D>;
    using underlying_type = T;
    using delete_type = D;
};

template <class T>
struct is_shared_ptr : std::false_type
{};

template <class T>
struct is_shared_ptr<std::shared_ptr<T>> : std::true_type
{};

template <typename T>
struct is_indexable  : std::false_type
{ };

template <typename T, typename A>
struct is_indexable<std::vector<T, A>> : std::true_type
{ };

template <typename T, std::size_t N>
struct is_indexable<std::array<T, N>> : std::true_type
{ };

namespace sly {
    typedef s32 StatusCode;
    typedef const char* ErrorMessage;

    static const StatusCode SUCCESS = (StatusCode)0;
    static const StatusCode UNKNOWN = (StatusCode)-1;

    template <typename T, typename = typename std::is_reference<T>::type>
    struct retval { };

    template <typename T>
    struct retval<T, std::false_type> {
    public:

        using type = T;
        using type_pointer = typename std::remove_reference<T>::type*;
        using type_reference = typename std::remove_reference<T>::type&;

        //using unique_ptr = typename is_unique_ptr<T>::value;

        retval(T value) :
            _value(std::move(value)),
            _statusCode(SUCCESS) {}

        retval(StatusCode code) : _statusCode(code) {}

        //retval(const retval&) = delete;
        retval& operator =(const retval&) = delete;
        ~retval() = default;

        template <class V, typename Y = T, typename std::enable_if_t<is_unique_ptr<Y>::value || is_shared_ptr<Y>::value>* = nullptr>
        V to() const {
            return reinterpret_cast<V>(_value.get());
        }

        template <class V, typename Y = T, typename std::enable_if_t<!is_unique_ptr<Y>::value && !is_shared_ptr<Y>::value>* = nullptr>
        V to() const {
            return reinterpret_cast<V>(_value);
        }

        template<typename V>
        operator V() const {
            return to<V>();
        }

        type result() { return _value; }
        type_reference ref() { return _value; }

        template <typename Y = T, typename std::enable_if_t<is_unique_ptr<Y>::value || is_shared_ptr<Y>::value>* = nullptr>
        operator type_reference() const {
            return *_value.get();
        }

        template <typename Y = T, typename std::enable_if_t<is_unique_ptr<Y>::value || is_shared_ptr<Y>::value>* = nullptr>
        operator typename unique_ptr_type<Y>::underlying_type&() {
            return *_value.get();
        }

        template <typename Y = T, typename std::enable_if_t<!is_unique_ptr<Y>::value && !is_shared_ptr<Y>::value>* = nullptr>
        operator type() { return result(); }

        template <typename Y = T, typename std::enable_if_t<std::is_pointer<Y>::value || is_unique_ptr<Y>::value || is_shared_ptr<Y>::value>* = nullptr>
        type_reference operator->() { return _value; }

        template <typename Y = T, typename std::enable_if_t<!std::is_pointer<Y>::value && !is_unique_ptr<Y>::value && !is_shared_ptr<Y>::value>* = nullptr>
        type_pointer operator->() { return &_value; }

        //template <typename N, typename Y = T, typename std::enable_if_t<std::is_pointer<Y>::value || is_unique_ptr<Y>::value || is_shared_ptr<Y>::value || is_indexable<Y>::value>* = nullptr>
        template<typename N>
        auto operator[](const N index) { return _value[index]; }

        StatusCode statusCode() const { return _statusCode; }

        bool_t succeeded() { return statusCode() == SUCCESS; }
        bool_t failed() { return statusCode() != SUCCESS; }

        template<typename X>
        X as() {
            return (X)result;
        }

    private:
        StatusCode _statusCode;
        T _value;
    };

    template <typename T>
    struct retval<T, std::true_type> {
    public:

        using type = T;
        using type_pointer = typename std::remove_reference<T>::type*;

        retval(T value) :
            _value(std::addressof(value)),
            _statusCode(SUCCESS) {}

        retval(StatusCode code) : _statusCode(code), _value(nullptr) {}

        retval(const retval&) = delete;
        retval& operator =(const retval&) = delete;
        ~retval() = default;

        T result() { return *_value; }
        operator T() { return result(); }
        type_pointer operator->() { return _value; }

        StatusCode statusCode() const { return _statusCode; }

        bool_t succeeded() { return statusCode() == SUCCESS; }
        bool_t failed() { return statusCode() != SUCCESS; }

        template<typename X>
        X as() {
            return (X)result();
        }

    private:
        StatusCode _statusCode;
        type_pointer _value;
    };

    template <>
    struct retval<void, std::false_type> {
    public:
        retval() : _statusCode(SUCCESS) {}
        retval(StatusCode code) : _statusCode(code) {}

        retval(const retval&) = delete;
        retval& operator =(const retval&) = delete;
        ~retval() = default;

        void value() const { }

        StatusCode statusCode() const { return _statusCode; }

        bool_t succeeded() { return statusCode() == SUCCESS; }
        bool_t failed() { return statusCode() != SUCCESS; }

        template<typename X>
        X as() {
            return (X)result;
        }

    private:
        StatusCode _statusCode;
    };

    template <typename T>
    retval<T> value(const T& value) {
        auto result = retval<T> ( std::move(value) ) ;
        return result;
    }

    template <typename T>
    retval<T&> reference(T& value) {
        return retval<T&> ( value );
    }

    retval<void> success();
    retval<void> failed();

    template <typename T>
    retval<T> failed(StatusCode statusCode, ErrorMessage message) {

        Expects(statusCode != SUCCESS);
        return retval<T> { statusCode };
    }

    template <typename T>
    retval<T> failed() {
        return retval<T> { UNKNOWN };
    }
}

max_stars_repo_name: Gibbeon/sly

max_issues_repo_name: Gibbeon/sly

max_forks_repo_name: Gibbeon/sly


text:
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_spmatrix.h>
#include <assert.h>

#include "graph_struct.h"

// pour améliorer la complexité:
// -> Les matrices de graph sont symétriques, on peut donc n'effectuer les opérations lors de la création du graph que sur une moitié de celui-ci
// -> Faire l'attribution des couleurs en même temps que la matrice d'adjacence

struct graph_t* graph__empty (int vertices){

  struct graph_t* g = malloc(sizeof(*g));
  g -> num_vertices = vertices;
  g -> t = gsl_spmatrix_alloc(vertices, vertices);
  g -> o = gsl_spmatrix_alloc(2, vertices);

  return g;
}

struct graph_t* graph_copy(struct graph_t * g){
  struct graph_t *copy = malloc(sizeof(struct graph_t));
  copy -> num_vertices = g -> num_vertices;
  copy -> o = gsl_spmatrix_alloc(2, g->num_vertices);
  gsl_spmatrix_memcpy(copy -> o, g -> o);
  copy -> t = gsl_spmatrix_alloc(g -> num_vertices, g -> num_vertices);
  gsl_spmatrix_memcpy(copy -> t, g -> t);
  return copy;
}

struct graph_t* make_graph(int m, char c){
  struct graph_t* g;

  switch (c){
  case 'c' :
    g = square_graph__empty(m);
    return g;
  case 't':
    g = triangular_graph__empty(m);
    return g;
  case 'h':
    g = hexa_graph__empty(m);
    return g;
  default:
    return NULL;
  }
}

struct graph_t* square_graph__empty(int m){
  int vertices = (m+1) * (m+1);
  struct graph_t* g = graph__empty (vertices);

  int line_length = m+1;

  for (int i = 0; i < vertices; i++){
    int index = i%line_length;

    // Adds right links
    if (index+1 < line_length)
      gsl_spmatrix_set(g -> t, i, i+1, 1);
    // Adds left links
    if (index-1 >= 0)
      gsl_spmatrix_set(g -> t, i, i-1, 1);
    // Adds down links
    if (i+line_length < vertices)
      gsl_spmatrix_set(g -> t, i, i+line_length, 1);
    // Adds up links
    if (i-line_length >= 0)
      gsl_spmatrix_set(g -> t, i, i-line_length, 1);
  }

  // Ownership of initial boxes
  for (int i = 0; i < line_length-1; i++){
    gsl_spmatrix_set(g -> o, 0, i, 1);
    gsl_spmatrix_set(g -> o, 0, vertices-1-i, 1);
    gsl_spmatrix_set(g -> o, 1, (i+1)*line_length, 1);
    gsl_spmatrix_set(g -> o, 1, vertices-1-(i+1)*line_length, 1);
  }

  return g;
}

struct graph_t* triangular_graph__empty(int m){
  int num_vertices = 6 * m * m;
  struct graph_t* g = graph__empty(num_vertices);
  int first_vertice_row=0;
  int row = 0;
  while(row < m){
    //Edges in a row
    for(int vertice = first_vertice_row; vertice < first_vertice_row + m + row; vertice++){
      gsl_spmatrix_set(g -> t, vertice, vertice + m + row, 1);
      gsl_spmatrix_set(g -> t, vertice + m + row, vertice, 1);
      gsl_spmatrix_set(g -> t, vertice, vertice + m + row + 1, 1);
      gsl_spmatrix_set(g -> t, vertice + m + row + 1, vertice, 1);

      gsl_spmatrix_set(g -> t, num_vertices - 1 - (vertice + m + row), num_vertices - 1 - vertice, 1);
      gsl_spmatrix_set(g -> t, num_vertices - 1 - vertice, num_vertices - 1 - (vertice + m + row ), 1);
      gsl_spmatrix_set(g -> t, num_vertices - 1 - (vertice + m + row + 1), num_vertices - 1 - vertice, 1);
      gsl_spmatrix_set(g -> t, num_vertices - 1 - vertice, num_vertices-1 - (vertice + m + row + 1), 1);
    }
    //Edges between 2 rows
    for(int vertice = first_vertice_row + m + row; vertice <= first_vertice_row + 2*(m + row); vertice++){
      gsl_spmatrix_set(g -> t, vertice, vertice + m + row + 1, 1);
      gsl_spmatrix_set(g -> t, vertice + m + row + 1, vertice, 1);

      gsl_spmatrix_set(g -> t, num_vertices - 1 - vertice, num_vertices - 1 - (vertice + m + row + 1), 1);
      gsl_spmatrix_set(g -> t, num_vertices - 1 - (vertice + m + row + 1), num_vertices - 1 - vertice, 1);
    }
    row++;
    first_vertice_row += 2*m + row*(row+1)/2;
  }
  //black starting tiles
  for(int vertice = 0; vertice < m; vertice ++){
    gsl_spmatrix_set(g -> o, BLACK, vertice, 1);
    gsl_spmatrix_set(g -> o, BLACK, vertice + m + 1, 1);

    gsl_spmatrix_set(g -> o, BLACK, num_vertices - 1 - vertice, 1);
    gsl_spmatrix_set(g -> o, BLACK, num_vertices - 1 - (vertice + m + 1), 1);
  }
  for(int row = 1; row < m/2 ; row++){
    gsl_spmatrix_set(g -> o, BLACK, m-1 + 2*row*m + row*(row+1), 1);
    gsl_spmatrix_set(g -> o, BLACK, m-1 + 2*row*m + row*(row+1) + m + row + 1, 1);

    gsl_spmatrix_set(g -> o, BLACK, num_vertices - 1 - (m-1 + 2*row*m + row*(row+1)), 1);
    gsl_spmatrix_set(g -> o, BLACK, num_vertices - 1 - (m-1 + 2*row*m + row*(row+1) + m + row + 1), 1);
  }
  //white starting tiles
  for(int row = 0; row < m; row++){
    gsl_spmatrix_set(g -> o, WHITE, m + 2*row*m + row*(row+1), 1);
    gsl_spmatrix_set(g -> o, WHITE, m + 2*row*m + row*(row+1) + m + row + 1, 1);

    gsl_spmatrix_set(g -> o, WHITE, num_vertices - 1 - (m + 2*row*m + row*(row+1)), 1);
    gsl_spmatrix_set(g -> o, WHITE, num_vertices - 1 - (m + 2*row*m + row*(row+1) + m + row + 1), 1);
  }
  for(int row = m/2+1; row < m; row++){
    gsl_spmatrix_set(g -> o, WHITE, m-1 + 2*row*m + row*(row+1), 1);
    gsl_spmatrix_set(g -> o, WHITE, m-1 + 2*row*m + row*(row+1) - (m + row), 1);

    gsl_spmatrix_set(g -> o, WHITE, num_vertices - 1 - (m-1 + 2*row*m + row*(row+1)), 1);
    gsl_spmatrix_set(g -> o, WHITE, num_vertices - 1 - (m-1 + 2*row*m + row*(row+1) - (m + row)),1);
  }
  return g;
}

struct graph_t* hexa_graph__empty(int m){
  int vertices = (m+1) * (m+1);
  struct graph_t* g = graph__empty (vertices);

  int line_length = m+1;

  for (int i = 0; i < vertices; i++){
    int index = i%line_length;

    // Adds right links
    if (index+1 < line_length)
      gsl_spmatrix_set(g -> t, i, i+1, 1);
    // Adds left links
    if (index-1 >= 0)
      gsl_spmatrix_set(g -> t, i, i-1, 1);
    // Adds down links
    if (i/line_length < line_length-1 && index > 0)
      gsl_spmatrix_set(g -> t, i, i+line_length - 1, 1);
    if (i+line_length < vertices && i/line_length < line_length-1)
      gsl_spmatrix_set(g -> t, i, i+line_length, 1);
    // Adds up links
    if (i/m > 0 && index < line_length-1)
      gsl_spmatrix_set(g -> t, i, i-line_length + 1, 1);
    if (i-line_length >= 0 && i/line_length > 0)
      gsl_spmatrix_set(g -> t, i, i-line_length, 1);

    // Ownership of initial boxes
    for (int i = 0; i < line_length-1; i++){
      gsl_spmatrix_set(g -> o, 0, i, 1);
      gsl_spmatrix_set(g -> o, 0, vertices-1-i, 1);
      gsl_spmatrix_set(g -> o, 1, (i+1)*line_length, 1); // Pour m = 3, l_l = 4 : 4, 8, 12
      gsl_spmatrix_set(g -> o, 1, vertices-1-(i+1)*line_length, 1); //
    }
  }
  return g;
}

void add__move_to_graph (struct graph_t* g, int vertice, int player){
  assert (g != NULL);

  gsl_spmatrix_set(g -> o, player, vertice, 1);
}

int get__box_owner (struct graph_t* g, int vertice){
  for (int p = 0; p < 2; p++){
    if (gsl_spmatrix_get(g -> o, p, vertice) == 1)
      return p;
  }
  return EMPTY_BOX;
}

void graph__free (struct graph_t* g){
  assert (g != NULL);

  gsl_spmatrix_free(g -> o);
  gsl_spmatrix_free(g -> t);

  free(g);
}

void print_graph (struct graph_t* g, char c){
  switch (c){
  case 'c' :
    square_graph__print (g);
    break;
  case 't':
    triangular_graph__print (g);
    break;
  case 'h':
    hexa_graph__print (g);
    break;
  default:
    printf("Don't find this type of graph\n");
    break;
  }
}

void square_graph__print (struct graph_t* g){
  assert (g != NULL);

  int n = sqrt(g -> num_vertices);

  for (int i = 0; i < n; i++){
    printf("\t");
    for (int j = 0; j < n; j++){
      int box_owner = get__box_owner (g, i*n + j);
      if (box_owner == 0)
	printf("\033[34;01m■ \033[00m");
      else if (box_owner == 1)
	printf("\033[31;01m■ \033[00m");
      else
	printf("■ ");
    }
    printf("\n");
  }
}

void triangular_graph__print (struct graph_t* g){(void) g;}

void hexa_graph__print (struct graph_t* g){
    assert (g != NULL);

  int n = sqrt(g -> num_vertices);
  char space[1000] = "";

  for (int i = 0; i < n; i++){
    printf("\t%s", space);

    for (int space = 0; space < i; space++){
      printf(" ");
    }
    for (int j = 0; j < n; j++){
      int box_owner = get__box_owner (g, i*n + j);
      if (box_owner == 0)
	printf("\033[34;01m⬢ \033[00m");
      else if (box_owner == 1)
	printf("\033[31;01m⬢ \033[00m");
      else
	printf("⬢ ");
    }
    printf("\n");
  }
}

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: ielomariala/Hex-Game


text:
/**
 *
 * @file example_dgelqf.c
 *
 *  PLASMA testing routines
 *  PLASMA is a software package provided by Univ. of Tennessee,
 *  Univ. of California Berkeley and Univ. of Colorado Denver
 *
 * @brief Example using LQ factorization
 *
 * @version 2.6.0
 * @author Bilel Hadri
 * @date 2010-11-15
 * @generated d Tue Jan  7 11:45:20 2014
 *
 **/
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>

#ifndef max
#define max(a, b) ((a) > (b) ? (a) : (b))
#endif
#ifndef min
#define min(a, b) ((a) < (b) ? (a) : (b))
#endif

int check_orthogonality(int, int, int, double*);
int check_factorization(int, int, double*, double*, int, double*);

int IONE=1;
int ISEED[4] = {0,0,0,1};   /* initial seed for dlarnv() */

int main ()
{

    int cores = 2;
    int M     = 10;
    int N     = 15;
    int LDA   = 10;
    int K = min(M, N);
    int info;
    int info_ortho, info_factorization;
    int i,j;
    int LDAxN = LDA*N;

    double *A1 = (double *)malloc(LDA*N*sizeof(double));
    double *A2 = (double *)malloc(LDA*N*sizeof(double));
    double *Q  = (double *)malloc(LDA*N*sizeof(double));
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ((!A1)||(!A2)||(!Q)){
        printf("Out of Memory \n ");
        return EXIT_SUCCESS;
    }

    /* Plasma Initialization */
    PLASMA_Init(cores);
    printf("-- PLASMA is initialized to run on %d cores. \n",cores);

    /* Allocate T */
    PLASMA_Alloc_Workspace_dgelqf(M, N, &T);

    /* Initialize A1 and A2 */
    LAPACKE_dlarnv_work(IONE, ISEED, LDAxN, A1);
    for (i = 0; i < M; i++)
        for (j = 0; j < N; j++)
            A2[LDA*j+i] = A1[LDA*j+i] ;

    /* Factorization QR of the matrix A2 */
    info = PLASMA_dgelqf(M, N, A2, LDA, T);

    /* Building the economy-size Q */
    memset((void*)Q, 0, LDA*N*sizeof(double));
    for (i = 0; i < K; i++)
        Q[LDA*i+i] = 1.0;

    PLASMA_dorglq(M, N, K, A2, LDA, T, Q, LDA);

    /* Check the orthogonality, factorization and the solution */
    info_ortho = check_orthogonality(M, N, LDA, Q);
    info_factorization = check_factorization(M, N, A1, A2, LDA, Q);

    if ((info_ortho != 0)|(info_factorization != 0)|(info != 0))
       printf("-- Error in DGELQF example ! \n");
    else
       printf("-- Run of DGELQF example successful ! \n");

    free(A1); free(A2); free(Q); free(T);

    PLASMA_Finalize();

    return EXIT_SUCCESS;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

int check_orthogonality(int M, int N, int LDQ, double *Q)
{
    double alpha, beta;
    double normQ;
    int info_ortho;
    int i;
    int minMN = min(M, N);
    double eps;
    double *work = (double *)malloc(minMN*sizeof(double));

    eps = LAPACKE_dlamch_work('e');
    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix USE DLASET?*/
    double *Id = (double *) malloc(minMN*minMN*sizeof(double));
    memset((void*)Id, 0, minMN*minMN*sizeof(double));
    for (i = 0; i < minMN; i++)
        Id[i*minMN+i] = (double)1.0;

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_dsyrk(CblasColMajor, CblasUpper, CblasTrans, N, M, alpha, Q, LDQ, beta, Id, N);
    else
        cblas_dsyrk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);

    normQ = LAPACKE_dlansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'u', minMN, Id, minMN, work);

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }

    free(work); free(Id);

    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the factorization QR
 */

int check_factorization(int M, int N, double *A1, double *A2, int LDA, double *Q)
{
    double Anorm, Rnorm;
    double alpha, beta;
    int info_factorization;
    int i,j;
    double eps;

    eps = LAPACKE_dlamch_work('e');

    double *Ql       = (double *)malloc(M*N*sizeof(double));
    double *Residual = (double *)malloc(M*N*sizeof(double));
    double *work              = (double *)malloc(max(M,N)*sizeof(double));

    alpha=1.0;
    beta=0.0;

    if (M >= N) {
        /* Extract the R */
        double *R = (double *)malloc(N*N*sizeof(double));
        memset((void*)R, 0, N*N*sizeof(double));
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);

        /* Perform Ql=Q*R */
        memset((void*)Ql, 0, M*N*sizeof(double));
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, (alpha), Q, LDA, R, N, (beta), Ql, M);
        free(R);
    }
    else {
        /* Extract the L */
        double *L = (double *)malloc(M*M*sizeof(double));
        memset((void*)L, 0, M*M*sizeof(double));
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);

    /* Perform Ql=LQ */
        memset((void*)Ql, 0, M*N*sizeof(double));
        cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, (alpha), L, M, Q, LDA, (beta), Ql, M);
        free(L);
    }

    /* Compute the Residual */
    for (i = 0; i < M; i++)
        for (j = 0 ; j < N; j++)
            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];

    Rnorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, Residual, M, work);
    Anorm = LAPACKE_dlange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A2, LDA, work);

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }

    if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else {
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(work); free(Ql); free(Residual);

    return info_factorization;
}


text:
//
//  mcrat.h
//  Global header file for all of MCRaT
//
//  Created by Tyler Parsotan on 7/23/21.
//

//define on and off for switches
#define ON  1
#define OFF 0

//define switches for MCRaT to continue a simulation 'c' or initalize a simulation
#define INITALIZE 'i'
#define CONTINUE  'c'

//define the codes for the different types of hydro simulations that we can use
#define FLASH           0
#define PLUTO_CHOMBO    1
#define PLUTO           2 //have separate compiler directive for file type
#define RIKEN           3 //need to replace this when done testing, just keeping it in the background for now and minimize its effect

//define PLUTO file types here
#define FILE_DBL         0
#define FILE_FLT         1
#define FILE_DBL_H5      2
#define FILE_FLT_H5      3
#define FILE_VTK         4

//define types of simulations that can be run
#define SCIENCE                         0
#define CYLINDRICAL_OUTFLOW             1
#define SPHERICAL_OUTFLOW               2
#define STRUCTURED_SPHERICAL_OUTFLOW    3

//define the geometries that we can handle
#define CARTESIAN   0 //2D, 3D
#define SPHERICAL   1 //2D, 3D
#define CYLINDRICAL 2 //2D
#define POLAR       3 //only in 3D, technically Cylindrical coordinates but PLUTO calls this POLAR

//define the hydro dimensions that we can handle
#define TWO                 0
#define TWO_POINT_FIVE      1
#define THREE               2

//define the types of things that we can assume for the thermal synchrotron emission and how we calculate the B field
#define INTERNAL_E  0
#define TOTAL_E     1
#define SIMULATION  2 //need to add this to statement below to take care of defaults

//define photon types
#define INJECTED_PHOTON 'i'
#define COMPTONIZED_PHOTON 'k'
#define CS_POOL_PHOTON 'p'
#define UNABSORBED_CS_PHOTON 'c'
#define REBINNED_PHOTON 'r'

extern const double C_LIGHT;
extern const double A_RAD;
extern const double PL_CONST;
extern const double K_B;
extern const double M_P;
extern const double THOM_X_SECT;
extern const double M_EL;
extern const double FINE_STRUCT;
extern const double CHARGE_EL;
extern const double R_EL;

#define STR_BUFFER 2000

struct photon
{
    char type; //was the photon injected as blackbody or wien, 'i', or was it emitted as cyclo-synchrotron or was it a cyclo-synchrotron photon that was compton scattered
    double p0; //E/c, 4 momentum is in lab frame
    double p1; // p_x
    double p2; //p_y
    double p3; //p_z
    double comv_p0; //E/c, 4 momentum is in comoving frame
    double comv_p1; // p_x
    double comv_p2; //p_y
    double comv_p3; //p_z
    double r0; //x in MCRaT coordinates these get saved
    double r1; //y
    double r2; //z
    double s0; //stokes I always 1
    double s1; //stokes Q/I +1 is in positive y_tilde coordinate (z_hat X photon 4 momentum)
    double s2; //stokes U/I
    double s3; //Stokes V/I
    double num_scatt;
    double weight; //each photon should have equal weight, sp this shouldnt matter, weight in mc.par file but across injections can have varying weights
    int nearest_block_index; //index that  allows for extraction of information of the hydro grid block that the photon si located within
} ; //structure to hold photon information

struct hydro_dataframe
{
    /*
     Coordinate System  |   Coordinate unit vector order (r0,r1,r2)/(v0,v1,v2)
     3D Cartesian       |       x, y, z
     3D Spherical       |       r, theta, phi
     3D Polar           |       r, phi, z
     2D Cartesian       |       x, z
     2D Cylindrical     |       r, z (phi) //in PLUTO its possible to save 3D vectors  in 2.5 dims where the final unit vector is phi hat
     2D Spherical       |       r, theta, (phi)
     */
    int num_elements; //number of elements in each array
    double *r0; //coodinates in hydro coodinate system that user provides,
    double *r1;
    double *r2;
    double *r0_size;//size of fluid elements
    double *r1_size;
    double *r2_size;
    double *r;//spherical coordinates of fluid elements for photon injection, mdifying fluid value purposes
    double *theta;
    double *v0; //velocity in hydro coordinate system
    double *v1;
    double *v2;
    double *dens;
    double *dens_lab;
    double *pres;
    double *temp;
    double *gamma;
    double *B0; //magentic field in hydro coordinate system
    double *B1;
    double *B2;

    //also hold global simulation information
    double r0_domain[2]; //holds min and max values of r0 coordinates for hydro sim
    double r1_domain[2];
    double r2_domain[2];
    double fps; //frames per second of the simulation
    int scatt_frame_number;
    int inj_frame_number;
    int last_frame;
    int increment_inj_frame; //the change in between each injection frame which may change if different number of fps is used in each portion of hydro code
    int increment_scatt_frame; //same as above expect for frames that photons are acattered in
}; // structure to hold all information for a given hydro simulation

//include all libraries needed
#include "hdf5.h"
#include "mpi.h"

#include <omp.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_histogram2d.h>
#include <gsl/gsl_integration.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <stdbool.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <glob.h>
#include <unistd.h>
#include <dirent.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_permutation.h>
#include <gsl/gsl_sort.h>
#include <gsl/gsl_sort_vector.h>

//include all other header files in mcrat
#include "mcrat_input.h"

#include "mclib.h"
#include "mclib_riken.h"
#include "mclib_pluto.h"
#include "mclib_flash.h"
#include "mc_cyclosynch.h"
#include "geometry.h"
#include "mcrat_scattering.h"
#include "mcrat_io.h"

//set default hydro v scale to be the speed of light
#define HYDRO_V_SCALE C_LIGHT

//take care of hydro scales in pressure, magnetic field, etc
#define HYDRO_P_SCALE   HYDRO_D_SCALE*HYDRO_V_SCALE*HYDRO_V_SCALE
#if B_FIELD_CALC==SIMULATION
    #define HYDRO_B_SCALE sqrt(4*M_PI*HYDRO_P_SCALE)
#endif

//take care of default PLUTO file types here
#ifdef PLUTO
    #ifndef PLUTO_FILETYPE
        #define PLUTO_FILETYPE  FILE_DBL
    #endif
#endif

//take care of synchrotron defaults here
#ifdef CYCLOSYNCHROTRON_SWITCH

    //if the percentage of max photon that will be used to create the energy bins isnt defined, define it to be 10%, also applies to emiting synch photons
    #ifndef CYCLOSYNCHROTRON_REBIN_E_PERC
        #define CYCLOSYNCHROTRON_REBIN_E_PERC 0.1
    #endif

    //if the polar angle bins that the rebinned synch photons isnt defined use 0.5 degree increments by default
    #ifndef CYCLOSYNCHROTRON_REBIN_ANG
        #define CYCLOSYNCHROTRON_REBIN_ANG 0.5
    #endif

    //if the azimuthal angle bins that the rebinned synch photons isnt defined use 0.5 degree increments by default
    #ifndef CYCLOSYNCHROTRON_REBIN_ANG_PHI
        #define CYCLOSYNCHROTRON_REBIN_ANG_PHI 10 //0.5
    #endif

    //if the user hasnt defined anything for how to calculate the B field, assume that they want it calculated from the total energy
    #ifndef B_FIELD_CALC
        #define B_FIELD_CALC TOTAL_E
    #endif
    //it is defined therefore see if EPSILON_B has been set and B_FIELD_CALC != SIMULATION
    #if B_FIELD_CALC == TOTAL_E || B_FIELD_CALC == INTERNAL_E
        //see if epsilon_b has been set
        #ifndef EPSILON_B
            //if not set it to be 0.5 by default
            #define EPSILON_B 0.5
        #endif
    #endif
#else
    //if its not defined set it to be off by default
    #define CYCLOSYNCHROTRON_SWITCH OFF

#endif

//take care of stokes switch and comv_switch and save_type defaults too
#ifndef STOKES_SWITCH
    #define STOKES_SWITCH   OFF
#endif

#ifndef COMV_SWITCH
    #define COMV_SWITCH   OFF
#endif

#ifndef SAVE_TYPE
    #define SAVE_TYPE   OFF
#endif

//throw errors during compilation if other switches are not defined
#ifndef SIM_SWITCH
#error Need to define hydro simulation type in mcrat_input.h file using SIM_SWITCH
#endif

#ifndef DIMENSIONS
#error Need to define hydro simulation dimensions in mcrat_input.h file using DIMENSIONS (should be set to two for this version of MCRaT)
#endif

#ifndef GEOMETRY
#error Need to define hydro simulation geometry in mcrat_input.h file using GEOMETRY
#endif

#ifndef HYDRO_L_SCALE
#error Need to define hydro simulation length scaling in mcrat_input.h file using HYDRO_L_SCALE
#endif

#ifndef HYDRO_D_SCALE
#error Need to define hydro simulation density scaling in mcrat_input.h file using HYDRO_D_SCALE
#endif

#ifndef MCPAR
#error Need to define name of MCRaT parameter file in mcrat_input.h file using MCPAR (it is typically called mc.par, see e.g. the MCRaT manual)
#endif

max_stars_repo_name: lazzati-astro/MCRaT

max_issues_repo_name: lazzati-astro/MCRaT

max_forks_repo_name: lazzati-astro/MCRaT


text:
#pragma once

#include "halley/text/halleystring.h"
#include "halley/utils/utils.h"
#include <gsl/span>

namespace Halley
{
	class Path
	{
	public:
		Path();
		Path(const char* name);
		Path(const std::string& name);
		Path(const String& name);

		Path(const Path& other) = default;
		Path(Path&& other) noexcept = default;
		Path& operator=(const Path& other) = default;
		Path& operator=(Path&& other) noexcept = default;

		Path& operator=(const std::string& other);
		Path& operator=(const String& other);

		Path getRoot() const;
		Path getFront(size_t n) const;
		Path getFilename() const;
		Path getStem() const;
		String getExtension() const;
		String getString() const;
		String getNativeString() const;
		String toString() const;

		size_t getNumberPaths() const;

		Path dropFront(int numberFolders) const;

		Path parentPath() const;
		Path replaceExtension(String newExtension) const;

		Path operator/(const char* other) const;
		Path operator/(const Path& other) const;
		Path operator/(const String& other) const;
		Path operator/(const std::string& other) const;

		bool operator==(const char* other) const;
		bool operator==(const String& other) const;
		bool operator==(const Path& other) const;

		bool operator!=(const Path& other) const;
		bool operator<(const Path& other) const;

		std::string string() const;

		static void writeFile(const Path& path, gsl::span<const gsl::byte> data);
		static void writeFile(const Path& path, const Bytes& data);
		static void writeFile(const Path& path, const String& data);
		static Bytes readFile(const Path& path);
		static void removeFile(const Path& path);

		Path makeRelativeTo(const Path& path) const;
		Path changeRelativeRoot(const Path& currentParent, const Path& newParent) const;

		bool isDirectory() const;
		bool isFile() const;
		bool isAbsolute() const;
		bool isEmpty() const;

	private:
		Vector<String> pathParts;
		void normalise();
		void setPath(const String& value);

		explicit Path(Vector<String> parts);
	};

	using TimestampedPath = std::pair<Path, int64_t>;
}

max_stars_repo_name: amrezzd/halley

max_issues_repo_name: amrezzd/halley

max_forks_repo_name: amrezzd/halley


text:
/*
 * Copyright 2020 Makani Technologies LLC
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *      http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#include "lib/util/random_matrix.h"

#include <assert.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_vector.h>
#include <stdint.h>

#include "common/c_math/linalg.h"

void MatArrRandomGaussian(int32_t nr, int32_t nc, const gsl_rng *rng,
                          double *d) {
  for (int32_t i = 0; i < nr; ++i) {
    for (int32_t j = 0; j < nc; ++j) {
      d[i*nc + j] = gsl_ran_gaussian(rng, 1.0);
    }
  }
}

void MatArrRandomOrthogonal(int32_t n, const gsl_rng *rng, double *d) {
  MatArrRandomGaussian(n, n, rng, d);
  gsl_matrix *Q = gsl_matrix_alloc((size_t)n, (size_t)n);
  gsl_matrix *R = gsl_matrix_alloc((size_t)n, (size_t)n);
  gsl_vector *tau = gsl_vector_alloc((size_t)n);

  gsl_matrix_view A = gsl_matrix_view_array(d, (size_t)n, (size_t)n);
  gsl_linalg_QR_decomp(&A.matrix, tau);
  gsl_linalg_QR_unpack(&A.matrix, tau, Q, R);
  gsl_matrix_memcpy(&A.matrix, Q);

  gsl_vector_free(tau);
  gsl_matrix_free(R);
  gsl_matrix_free(Q);
}

void MatArrRandomWithRank(int32_t nr, int32_t nc, int32_t rank,
                          const gsl_rng *rng, double *d) {
  assert(rank <= nr && rank <= nc);

  MatArrZero(nr, nc, d);
  double *u = (double *)malloc((size_t)nr * sizeof(*d));
  double *v = (double *)malloc((size_t)nc * sizeof(*d));
  for (int32_t k = 0; k < rank; ++k) {
    MatArrRandomGaussian(1, nr, rng, u);
    MatArrRandomGaussian(1, nc, rng, v);

    for (int32_t i = 0; i < nr; ++i) {
      for (int32_t j = 0; j < nc; ++j) {
        d[i*nc + j] += u[i] * v[j];
      }
    }
  }
  free(u);
  free(v);
}


text:
#pragma once

#include <gsl.h>

#include <boost/graph/graph_traits.hpp>
#include <boost/range/adaptor/filtered.hpp>

#include "VesselForestData.h"

namespace crimson
{
struct ActiveVesselsFilter {
    ActiveVesselsFilter() = default;
    ActiveVesselsFilter(const VesselForestData* owner)
        : _owner(owner)
    {
    }

    bool operator()(const VesselForestData::VesselPathUIDType& uid) const { return _owner->getVesselUsedInBlending(uid); }

private:
    const VesselForestData* _owner = nullptr;
};

struct ActiveBooleanOperationsFilter {
    ActiveBooleanOperationsFilter() = default;
    ActiveBooleanOperationsFilter(const VesselForestData* owner)
        : _owner(owner)
    {
    }

    bool operator()(const VesselForestData::BooleanOperationInfo& bop) const
    {
        return _owner->getVesselUsedInBlending(bop.vessels.first) && _owner->getVesselUsedInBlending(bop.vessels.second);
    }

private:
    const VesselForestData* _owner;
};

struct OutEdgeIteratorFilter {
    OutEdgeIteratorFilter() = default;
    OutEdgeIteratorFilter(const VesselForestData::VesselPathUIDType* uid)
        : _uid(uid)
    {
    }

    bool operator()(const VesselForestData::BooleanOperationInfo& bop) const
    {
        return bop.vessels.first == *_uid || bop.vessels.second == *_uid;
    }

private:
    const VesselForestData::VesselPathUIDType* _uid = nullptr;
};

struct ExtractOutEdgeFromBooleanOperationInfo {
    ExtractOutEdgeFromBooleanOperationInfo() = default;
    ExtractOutEdgeFromBooleanOperationInfo(const VesselForestData::VesselPathUIDType* sourceVertexUID)
        : _sourceVertexUID(sourceVertexUID)
    {
    }

    VesselForestData::VesselPathUIDPair operator()(const VesselForestData::BooleanOperationInfo& bop) const
    {
        if (*_sourceVertexUID == bop.vessels.first) {
            return bop.vessels;
        }
        Expects(*_sourceVertexUID == bop.vessels.second);
        return std::make_pair(bop.vessels.second, bop.vessels.first);
    }

private:
    const VesselForestData::VesselPathUIDType* _sourceVertexUID;
};

} // namespace crimson

namespace boost
{
template <>
struct graph_traits<crimson::VesselForestData> {
    /////////////////////////////////
    // Graph concept
    /////////////////////////////////
    using vertex_descriptor = crimson::VesselForestData::VesselPathUIDType;
    using edge_descriptor = crimson::VesselForestData::VesselPathUIDPair;

    using directed_category = boost::undirected_tag;
    using edge_parallel_category = boost::disallow_parallel_edge_tag;

    struct traversal_category : boost::vertex_list_graph_tag, boost::incidence_graph_tag {
    };

    /////////////////////////////////
    // Incidence graph concept
    /////////////////////////////////
    using out_edge_iterator = boost::transform_iterator<
        crimson::ExtractOutEdgeFromBooleanOperationInfo,
        boost::filter_iterator<
            crimson::OutEdgeIteratorFilter,
            boost::filter_iterator<crimson::ActiveBooleanOperationsFilter,
                                   crimson::VesselForestData::BooleanOperationContainerType::const_iterator>>>;
    using degree_size_type = std::size_t;

    // out_edges(v, g)
    // source(e, g)
    // target(e, g)
    // out_degree(v, g)

    /////////////////////////////////
    // Vertex list graph concept
    /////////////////////////////////
    using vertex_iterator = boost::filter_iterator<crimson::ActiveVesselsFilter,
                                                   crimson::VesselForestData::VesselPathUIDContainerType::const_iterator>;
    using vertices_size_type = crimson::VesselForestData::VesselPathUIDContainerType::size_type;

    // vertices(g)
    // num_vertices(g)

    static vertex_descriptor null_vertex() { return {}; }
};
} // namespace boost

namespace crimson
{

inline boost::graph_traits<crimson::VesselForestData>::vertex_descriptor
source(const boost::graph_traits<crimson::VesselForestData>::edge_descriptor& e, const crimson::VesselForestData& g)
{
    return e.first;
}

inline boost::graph_traits<crimson::VesselForestData>::vertex_descriptor
target(const boost::graph_traits<crimson::VesselForestData>::edge_descriptor& e, const crimson::VesselForestData& g)
{
    return e.second;
}

inline std::pair<boost::graph_traits<crimson::VesselForestData>::out_edge_iterator,
                 boost::graph_traits<crimson::VesselForestData>::out_edge_iterator>
out_edges(const boost::graph_traits<crimson::VesselForestData>::vertex_descriptor& v, const crimson::VesselForestData& g)
{
    auto filtIter1 = boost::make_filter_iterator(crimson::ActiveBooleanOperationsFilter{&g}, g.getBooleanOperations().begin(),
                                                 g.getBooleanOperations().end());
    auto endFiltIter1 = boost::make_filter_iterator(crimson::ActiveBooleanOperationsFilter{&g}, g.getBooleanOperations().end(),
                                                    g.getBooleanOperations().end());

    auto filtIter2 = boost::make_filter_iterator(crimson::OutEdgeIteratorFilter{&v}, filtIter1, endFiltIter1);
    auto endFiltIter2 = boost::make_filter_iterator(crimson::OutEdgeIteratorFilter{&v}, endFiltIter1, endFiltIter1);

    return {boost::make_transform_iterator(filtIter2, crimson::ExtractOutEdgeFromBooleanOperationInfo{&v}),
            boost::make_transform_iterator(endFiltIter2, crimson::ExtractOutEdgeFromBooleanOperationInfo{&v})};
}

inline boost::graph_traits<crimson::VesselForestData>::degree_size_type
out_degree(const boost::graph_traits<crimson::VesselForestData>::vertex_descriptor& v, const crimson::VesselForestData& g)
{
    auto iterPair = out_edges(v, g);
    return static_cast<boost::graph_traits<crimson::VesselForestData>::degree_size_type>(
        std::distance(iterPair.first, iterPair.second));
}

inline std::pair<boost::graph_traits<crimson::VesselForestData>::vertex_iterator,
                 boost::graph_traits<crimson::VesselForestData>::vertex_iterator>
vertices(const crimson::VesselForestData& g)
{
    auto iter = boost::make_filter_iterator(crimson::ActiveVesselsFilter{&g}, g.getVessels().begin(), g.getVessels().end());
    auto endIter = boost::make_filter_iterator(crimson::ActiveVesselsFilter{&g}, g.getVessels().end(), g.getVessels().end());

    return {iter, endIter};
}

inline boost::graph_traits<crimson::VesselForestData>::vertices_size_type num_vertices(const crimson::VesselForestData& g)
{
    auto iterPair = vertices(g);
    return static_cast<boost::graph_traits<crimson::VesselForestData>::vertices_size_type>(
        std::distance(iterPair.first, iterPair.second));
}

// static void _test_vertex_list_graph_concept()
//{
//    BOOST_CONCEPT_ASSERT((boost::IncidenceGraphConcept<crimson::VesselForestData>));
//    BOOST_CONCEPT_ASSERT((boost::VertexListGraphConcept<crimson::VesselForestData>));
//}

} // namespace crimson

max_stars_repo_name: carthurs/CRIMSONGUI

max_issues_repo_name: carthurs/CRIMSONGUI

max_forks_repo_name: carthurs/CRIMSONGUI


text:
#pragma once

#include <cstdint>
#include <functional>
#include <memory>
#include <string>
#include <tuple>
#include <utility>
#include <vector>

#include <gsl/gsl>
#include <nonstd/optional.hpp>

#include "chainerx/array_body.h"
#include "chainerx/array_fwd.h"
#include "chainerx/device.h"
#include "chainerx/dtype.h"
#include "chainerx/graph.h"
#include "chainerx/macro.h"
#include "chainerx/shape.h"

namespace chainerx {

class BackwardContext;
class Device;

using BackwardFunction = std::function<void(BackwardContext&)>;

namespace internal {

class ArrayNode;
class OpNode;

struct ArrayProps {
    explicit ArrayProps(const Array& array);
    explicit ArrayProps(const ArrayNode& array_node);
    explicit ArrayProps(const ArrayBody& array_body);

    Shape shape;
    Dtype dtype;
    Device& device;
};

class OpNodeBackwardEntry {
public:
    OpNodeBackwardEntry(OpNode& op_node, std::vector<size_t> input_array_node_indices, BackwardFunction backward_func);

    OpNode& op_node() const { return op_node_; }

    size_t input_array_node_count() const { return input_array_node_indices_.size(); }

    const std::vector<size_t>& input_array_node_indices() const { return input_array_node_indices_; }

    const BackwardFunction& backward_func() const { return backward_func_; }

private:
    friend class OpNode;

    OpNode& op_node_;

    // The index mapping from local (this backward function) to global (op node).
    // Can be unset if the input array does not require grad.
    std::vector<size_t> input_array_node_indices_;

    BackwardFunction backward_func_;
};

// Creates an output array node at the specified index and adds edges between the output array node and the op node.
// Undefined behavior if the output array node already exists.
// This function is used by BackwardContext::GetRetainedOutput().
std::shared_ptr<ArrayNode> FabricateOutputArrayNode(std::shared_ptr<OpNode> op_node, size_t output_array_node_index);

class OpNode {
public:
    // Creates a new op node that has output array nodes corresponding to the given outputs.
    static std::shared_ptr<OpNode> CreateWithOutputArrayNodes(
            std::string name, BackpropId backprop_id, size_t input_count, const std::vector<ConstArrayRef>& outputs);

    OpNode(const OpNode&) = delete;
    OpNode(OpNode&&) = delete;
    OpNode& operator=(const OpNode&) = delete;
    OpNode& operator=(OpNode&&) = delete;

    OpNodeBackwardEntry& RegisterBackwardFunction(
            std::vector<std::tuple<size_t, std::shared_ptr<ArrayNode>>> input_array_nodes, BackwardFunction backward_func);

    // Adds links to input array nodes of other graphs.
    // The size of the vector must be equal to the number of inputs.
    void AddEdgesToInputArrayNodesOfOuterGraph(
            const BackpropId& outer_backprop_id, std::vector<std::shared_ptr<ArrayNode>> outer_graphs_input_array_nodes);

    // Adds links to output array nodes of other graphs.
    // The size of the vector must be equal to the number of outputs.
    void AddEdgesToOutputArrayNodesOfOuterGraph(
            const BackpropId& outer_backprop_id, std::vector<std::shared_ptr<ArrayNode>> outer_graphs_output_array_nodes);

    void Unchain() {
        backward_entries_.clear();
        std::fill(input_array_nodes_.begin(), input_array_nodes_.end(), std::shared_ptr<ArrayNode>{});
        AssertConsistency();
    }

    bool HasInputArrayNode(size_t input_index) const { return input_array_nodes_[input_index] != nullptr; }

    std::string name() const { return name_; }

    std::vector<std::shared_ptr<ArrayNode>>& input_array_nodes();

    const std::vector<std::shared_ptr<ArrayNode>>& input_array_nodes() const;

    gsl::span<OpNodeBackwardEntry> backward_entries() { return backward_entries_; }

    gsl::span<const OpNodeBackwardEntry> backward_entries() const { return backward_entries_; }

    size_t input_array_node_count() const { return input_array_nodes_.size(); }

    size_t output_array_node_count() const { return output_array_props_.size(); }

    int64_t rank() const { return rank_; }

    BackpropId backprop_id() const { return backprop_id_; }

    const ArrayProps& GetOutputArrayProps(size_t i) const {
        CHAINERX_ASSERT(i < output_array_props_.size());
        return output_array_props_[i];
    }

    // Returns the list of output array nodes on "this" graph.
    const std::vector<nonstd::optional<std::weak_ptr<ArrayNode>>>& output_array_nodes() const { return output_array_nodes_; }

    // Returns the list of output array nodes on "this" graph.
    std::vector<nonstd::optional<std::weak_ptr<ArrayNode>>>& output_array_nodes() { return output_array_nodes_; }

    // Returns the input array nodes of all graphs.
    const std::vector<std::tuple<BackpropId, std::vector<std::shared_ptr<ArrayNode>>>>& outer_graphs_input_array_nodes() const {
        return outer_graphs_input_array_nodes_;
    }

    // Returns the output array nodes of all graphs.
    const std::vector<std::tuple<BackpropId, std::vector<std::shared_ptr<ArrayNode>>>>& outer_graphs_output_array_nodes() const {
        return outer_graphs_output_array_nodes_;
    }

private:
    OpNode(std::string name, BackpropId backprop_id, size_t input_array_node_count);

    void AssertConsistency() const;

    std::string name_;

    // Backprop ID.
    // Backprop ID is also held in the first entry of output_array_nodes_, but the reference to it may be invalidated, whereas this member
    // is stable during the lifetime of this OpNode instance.
    BackpropId backprop_id_;

    int64_t rank_{0};

    // List of input array nodes.
    std::vector<std::shared_ptr<ArrayNode>> input_array_nodes_;

    // List of output array nodes of this graph.
    std::vector<nonstd::optional<std::weak_ptr<ArrayNode>>> output_array_nodes_;

    // List of input/output array nodes of outer graphs.
    // Outer graphs refer to graphs with lower ordinals.
    // Each entry is a pair of backprop ID and list of input/output array nodes.
    std::vector<std::tuple<BackpropId, std::vector<std::shared_ptr<ArrayNode>>>> outer_graphs_input_array_nodes_;
    std::vector<std::tuple<BackpropId, std::vector<std::shared_ptr<ArrayNode>>>> outer_graphs_output_array_nodes_;

    // Array props of output array nodes. This is used for creating dummy gradients.
    std::vector<ArrayProps> output_array_props_;

    std::vector<OpNodeBackwardEntry> backward_entries_;
};

}  // namespace internal
}  // namespace chainerx

max_issues_repo_name: nolfwin/chainer

max_forks_repo_name: nolfwin/chainer