text:
/* multifit/multiwlinear.c
 *
 * Copyright (C) 2015 Patrick Alken
 *
 * 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 <gsl/gsl_errno.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

#include "linear_common.c"

/* General weighted case */

static int
multifit_wlinear_svd (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      double tol,
                      int balance,
                      size_t * rank,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double *chisq, gsl_multifit_linear_workspace * work)
{
  if (X->size1 != y->size)
    {
      GSL_ERROR
        ("number of observations in y does not match rows of matrix X",
         GSL_EBADLEN);
    }
  else if (X->size2 != c->size)
    {
      GSL_ERROR ("number of parameters c does not match columns of matrix X",
                 GSL_EBADLEN);
    }
  else if (w->size != y->size)
    {
      GSL_ERROR ("number of weights does not match number of observations",
                 GSL_EBADLEN);
    }
  else if (cov->size1 != cov->size2)
    {
      GSL_ERROR ("covariance matrix is not square", GSL_ENOTSQR);
    }
  else if (c->size != cov->size1)
    {
      GSL_ERROR
        ("number of parameters does not match size of covariance matrix",
         GSL_EBADLEN);
    }
  else if (X->size1 != work->n || X->size2 != work->p)
    {
      GSL_ERROR
        ("size of workspace does not match size of observation matrix",
         GSL_EBADLEN);
    }
  else
    {
      const size_t n = X->size1;
      const size_t p = X->size2;

      size_t i, j, p_eff;

      gsl_matrix *A = work->A;
      gsl_matrix *Q = work->Q;
      gsl_matrix *QSI = work->QSI;
      gsl_vector *S = work->S;
      gsl_vector *t = work->t;
      gsl_vector *xt = work->xt;
      gsl_vector *D = work->D;

      /* Scale X,  A = sqrt(w) X */

      gsl_matrix_memcpy (A, X);

      for (i = 0; i < n; i++)
        {
          double wi = gsl_vector_get (w, i);

          if (wi < 0)
            wi = 0;

          {
            gsl_vector_view row = gsl_matrix_row (A, i);
            gsl_vector_scale (&row.vector, sqrt (wi));
          }
        }

      /* Balance the columns of the matrix A if requested */

      if (balance)
        {
          gsl_linalg_balance_columns (A, D);
        }
      else
        {
          gsl_vector_set_all (D, 1.0);
        }

      /* Decompose A into U S Q^T */

      gsl_linalg_SV_decomp_mod (A, QSI, Q, S, xt);

      /* Solve sqrt(w) y = A c for c, by first computing t = sqrt(w) y */

      for (i = 0; i < n; i++)
        {
          double wi = gsl_vector_get (w, i);
          double yi = gsl_vector_get (y, i);
          if (wi < 0)
            wi = 0;
          gsl_vector_set (t, i, sqrt (wi) * yi);
        }

      gsl_blas_dgemv (CblasTrans, 1.0, A, t, 0.0, xt);

      /* Scale the matrix Q,  Q' = Q S^-1 */

      gsl_matrix_memcpy (QSI, Q);

      {
        double s0 = gsl_vector_get (S, 0);
        p_eff = 0;

        for (j = 0; j < p; j++)
          {
            gsl_vector_view column = gsl_matrix_column (QSI, j);
            double sj = gsl_vector_get (S, j);
            double alpha;

            if (sj <= tol * s0)
              {
                alpha = 0.0;
              }
            else
              {
                alpha = 1.0 / sj;
                p_eff++;
              }

            gsl_vector_scale (&column.vector, alpha);
          }

        *rank = p_eff;
      }

      gsl_vector_set_zero (c);

      /* solution */
      gsl_blas_dgemv (CblasNoTrans, 1.0, QSI, xt, 0.0, c);

      /* unscale the balancing factors */
      gsl_vector_div (c, D);

      /* compute chisq, from residual r = y - X c */
      {
        double r2 = 0;

        for (i = 0; i < n; i++)
          {
            double yi = gsl_vector_get (y, i);
            double wi = gsl_vector_get (w, i);
            gsl_vector_const_view row = gsl_matrix_const_row (X, i);
            double y_est, ri;
            gsl_blas_ddot (&row.vector, c, &y_est);
            ri = yi - y_est;
            r2 += wi * ri * ri;
          }

        *chisq = r2;

        /* Form covariance matrix cov = (X^T W X)^-1 = (Q S^-1) (Q S^-1)^T */

        for (i = 0; i < p; i++)
          {
            gsl_vector_view row_i = gsl_matrix_row (QSI, i);
            double d_i = gsl_vector_get (D, i);

            for (j = i; j < p; j++)
              {
                gsl_vector_view row_j = gsl_matrix_row (QSI, j);
                double d_j = gsl_vector_get (D, j);
                double s;

                gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);

                gsl_matrix_set (cov, i, j, s / (d_i * d_j));
                gsl_matrix_set (cov, j, i, s / (d_i * d_j));
              }
          }
      }

      return GSL_SUCCESS;
    }
}

int
gsl_multifit_wlinear (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double *chisq, gsl_multifit_linear_workspace * work)
{
  size_t rank;
  int status = gsl_multifit_wlinear_tsvd(X, w, y, GSL_DBL_EPSILON, c, cov, chisq, &rank, work);

  return status;
}

int
gsl_multifit_wlinear_tsvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           const double tol,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double * chisq,
                           size_t * rank,
                           gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (y->size != n)
    {
      GSL_ERROR("number of observations in y does not match matrix",
                GSL_EBADLEN);
    }
  else if (w->size != n)
    {
      GSL_ERROR("number of weights in w does not match matrix",
                GSL_EBADLEN);
    }
  else if (p != c->size)
    {
      GSL_ERROR ("number of parameters c does not match matrix",
                 GSL_EBADLEN);
    }
  else if (tol <= 0)
    {
      GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
    }
  else
    {
      int status;
      double rnorm, snorm;
      gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
      gsl_vector_view b = gsl_vector_subvector(work->t, 0, n);

      /* compute A = sqrt(W) X, b = sqrt(W) y */
      status = gsl_multifit_linear_applyW(X, w, y, &A.matrix, &b.vector);
      if (status)
        return status;

      /* compute SVD of A */
      status = gsl_multifit_linear_bsvd(&A.matrix, work);
      if (status)
        return status;

      status = multifit_linear_solve(X, &b.vector, tol, 0.0, rank,
                                     c, &rnorm, &snorm, work);
      if (status)
        return status;

      *chisq = rnorm * rnorm;

      /* variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */
      {
        const size_t p = X->size2;
        size_t i, j;
        gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
        gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);

        for (i = 0; i < p; i++)
          {
            gsl_vector_view row_i = gsl_matrix_row (&QSI.matrix, i);
            double d_i = gsl_vector_get (&D.vector, i);

            for (j = i; j < p; j++)
              {
                gsl_vector_view row_j = gsl_matrix_row (&QSI.matrix, j);
                double d_j = gsl_vector_get (&D.vector, j);
                double s;

                gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);

                gsl_matrix_set (cov, i, j, s / (d_i * d_j));
                gsl_matrix_set (cov, j, i, s / (d_i * d_j));
              }
          }
      }
    }

  return GSL_SUCCESS;
}

int
gsl_multifit_wlinear_svd (const gsl_matrix * X,
                          const gsl_vector * w,
                          const gsl_vector * y,
                          double tol,
                          size_t * rank,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double *chisq, gsl_multifit_linear_workspace * work)
{
  int status  = multifit_wlinear_svd (X, w, y, tol, 1, rank, c,
                                      cov, chisq, work);
  return status;

}

int
gsl_multifit_wlinear_usvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           double tol,
                           size_t * rank,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double *chisq, gsl_multifit_linear_workspace * work)
{
  int status  = multifit_wlinear_svd (X, w, y, tol, 0, rank, c,
                                      cov, chisq, work);
  return status;

}


text:
#ifndef EXAMPLE_GENERATOR_H
#define EXAMPLE_GENERATOR_H

#include <opencv/cv.h>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

#include "helper/bounding_box.h"
#include "loader/loader_imagenet_det.h"
#include "loader/video.h"
#include "helper/Common.h"
#include "helper/CommonCV.h"
#include "helper/Constants.h"

#include <stdlib.h>     /* srand, rand */
#include <time.h>       /* time */
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h> /* GAUSSIAN*/

struct BBParams {
  double lambda_shift;
  double lambda_scale;
  double min_scale;
  double max_scale;
};

// Generates additional training examples by taking random crops of the target object,
// causing apparent translation and scale changes.
class ExampleGenerator
{
public:
  ExampleGenerator(const double lambda_shift, const double lambda_scale,
                   const double min_scale, const double max_scale);

  // Set up to train on the previous and current image, and the previous and current bounding boxes.
  void Reset(const BoundingBox& bbox_prev, const BoundingBox& bbox_curr,
             const cv::Mat& image_prev, const cv::Mat& image_curr);

  // Shift the whole bounding box for the current frame
  // (simulates camera motion)
  void MakeTrainingExampleBBShift(const bool visualize_example,
                                  cv::Mat* image_rand_focus,
                                  cv::Mat* target_pad,
                                  BoundingBox* bbox_gt_scaled) const;
  void MakeTrainingExampleBBShift(cv::Mat* image_rand_focus,
                                  cv::Mat* target_pad,
                                  BoundingBox* bbox_gt_scaled) const;

  // Focus the current image at the location of the previous frame's
  // bounding box (true motion).
  void MakeTrueExample(cv::Mat* image_focus, cv::Mat* target_pad,
                       BoundingBox* bbox_gt_scaled) const;
  void MakeTrueExampleTight(cv::Mat* image_focus, cv::Mat* target_tight,
                      BoundingBox* bbox_gt_scaled) const;

  // Make batch_size training examples according to the input parameters.
  void MakeTrainingExamples(const int num_examples, std::vector<cv::Mat>* images,
                            std::vector<cv::Mat>* targets,
                            std::vector<BoundingBox>* bboxes_gt_scaled);

  // Helper function to generate one moved box, non-gaussian version
  // trans_range, scale_range are defaults for uniform sampling of uniform samples
  static BoundingBox GenerateOneRandomCandidate(BoundingBox &bbox, gsl_rng* rng, int W, int H,
                                                const string method = "uniform", const double trans_range = POS_TRANS_RANGE, const double scale_range = POS_SCALE_RANGE,
                                                const double sd_x = SD_X, const double sd_y = SD_Y, const double sd_scale = SD_SCALE);

  // Make candidates given one frame
  // candidates: 200 neg and 50 pos candidates
  // labels: vector of scalar 1 means pos and 0 means neg
  // fixed number output, candidates and labels will have size == num_pos + num_neg
  void MakeCandidatesAndLabels(vector<Mat> *candidates, vector<double> *labels,
                               const int num_pos = POS_CANDIDATES,
                               const int num_neg = NEG_CANDIDATES);

  // Actual woker under MakeCandidatesAndLabels
  void MakeCandidatesAndLabelsBBox(vector<BoundingBox> *candidate_bboxes, vector<double> *labels,
                                   const int num_pos = POS_CANDIDATES,
                                   const int num_neg = NEG_CANDIDATES);

  void MakeCandidatesPos(vector<BoundingBox> *candidates, const int num = POS_CANDIDATES,
                                  const string method = "gaussian", const double trans_range = POS_TRANS_RANGE, const double scale_range = POS_SCALE_RANGE,
                                  const double sd_x = SD_X, const double sd_y = SD_Y, const double sd_scale = SD_SCALE, const double pos_iou_th = POS_IOU_TH );

  void MakeCandidatesNeg(vector<BoundingBox> *candidates, const int num = NEG_CANDIDATES,
                                  const string method = "uniform", const double trans_range = 2 * SD_X, const double scale_range = POS_SCALE_RANGE,
                                  const double sd_x = SD_X, const double sd_y = SD_Y, const double sd_scale = SD_SCALE);

  void set_indices(const int video_index, const int frame_index) {
    video_index_ = video_index; frame_index_ = frame_index;
  }

private:
  void MakeTrainingExampleBBShift(const bool visualize_example,
                                  const BBParams& bbparams,
                                  cv::Mat* image_rand_focus,
                                  cv::Mat* target_pad,
                                  BoundingBox* bbox_gt_scaled) const;

  void VisualizeExample(const cv::Mat& target_pad,
                        const cv::Mat& image_rand_focus,
                        const BoundingBox& bbox_gt_scaled) const;

  void get_default_bb_params(BBParams* default_params) const;

  // To generate synethic examples, shift the bounding box by an exponential with the given lambda parameter.
  double lambda_shift_;

  // To generate synethic examples, shift the bounding box by an exponential with the given lambda parameter.
  double lambda_scale_;

  // Do not scale the synthetic examples by more than max_scale_ shrink them by more than min_scale_.
  double min_scale_;
  double max_scale_;

  // Current training image.
  cv::Mat image_curr_;

  // Location of the target within the current and previous images.
  BoundingBox bbox_curr_gt_;
  BoundingBox bbox_prev_gt_;

  // Cropped and scaled image of the target object from the previous image.
  cv::Mat target_pad_;
  // tight previous target in previous image
  cv::Mat target_tight_;

  // Video and frame index from which the current example was generated.
  // These values are only used when saving images to a file, to assign them
  // a unique identifier.
  int video_index_;
  int frame_index_;

  gsl_rng* rng_;

  static std::mt19937 engine_;
};

#endif // EXAMPLE_GENERATOR_H

ext: h


text:
#ifndef _FUNCTIONS_
#define _FUNCTIONS_

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>
#include <limits.h>

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

#include "structs.c"

#define PATH_LENGTH 32
#define LINE_LENGTH 50000

#define FRECHET 1
#define DTW 2

#define THRESHOLD 0.0000005

// Declare functions of functions.c file
int is_centroid(int *, int , int );

double rand_gaussian(void);

int centroids_transposition(segment **, int *, int *, int , double **);

double **create_distances_array(segment **, int , int );

#endif

max_stars_repo_path: Part 3/Roads/src/functions.h

max_issues_repo_path: Part 3/Roads/src/functions.h

max_forks_repo_path: Part 3/Roads/src/functions.h


text:
#pragma once

#include "types/type_util.h"

#include <expected.hpp>
#include <gsl/gsl>
#include <sstream>
#include <variant>

namespace util {

template<typename... Args>
std::string
concat(const Args&... args)
{
  std::stringstream ss;
  (ss << ... << args);
  return ss.str();
}

namespace detail {
template<typename T>
T
str_to(const std::string& str, std::size_t* pos, int base)
{
  static_assert(AlwaysFalse<T>::value, "try_parse_number not defined for type T!");
}

#define GEN_STR_TO_SPECIALIZATION(type, func)                                                      \
  template<>                                                                                \
  inline type str_to(const std::string& str, std::size_t* pos, int base)                           \
  {                                                                                \
    const auto parsed = func(str, pos, base);                                                      \
    if (parsed > std::numeric_limits<type>::max()) {                                               \
      throw std::out_of_range{ "Value is out of range for type " #type };                          \
    }                                                                                \
    return gsl::narrow_cast<type>(parsed);                                                         \
  }

GEN_STR_TO_SPECIALIZATION(uint8_t, std::stoul)
GEN_STR_TO_SPECIALIZATION(int8_t, std::stol)
GEN_STR_TO_SPECIALIZATION(uint16_t, std::stoul)
GEN_STR_TO_SPECIALIZATION(int16_t, std::stol)
GEN_STR_TO_SPECIALIZATION(uint32_t, std::stoul)
GEN_STR_TO_SPECIALIZATION(int32_t, std::stol)
GEN_STR_TO_SPECIALIZATION(uint64_t, std::stoull)
GEN_STR_TO_SPECIALIZATION(int64_t, std::stoll)

template<>
inline float
str_to(const std::string& str, std::size_t* pos, int)
{
  return std::stof(str, pos);
}

template<>
inline double
str_to(const std::string& str, std::size_t* pos, int)
{
  return std::stod(str, pos);
}

template<>
inline long double
str_to(const std::string& str, std::size_t* pos, int)
{
  return std::stold(str, pos);
}

}

/**
 * Wrapper around the std::stoX functions. Tries to parse the given string as a number
 * of the given type and returns a tl::expected with the resulting number of the reason
 * for failure. The arguments are the same as for the std functions, with the exception that
 * base is only valid if T is an integer type
 */
template<typename T>
tl::expected<T, std::variant<std::invalid_argument, std::out_of_range>>
try_parse_number(const std::string& str, std::size_t* pos = nullptr, int base = 10)
{
  try {
    const auto number = detail::str_to<T>(str, pos, base);
    return { number };
  } catch (const std::invalid_argument& ex) {
    return tl::make_unexpected(std::variant<std::invalid_argument, std::out_of_range>{ ex });
  } catch (const std::out_of_range& ex) {
    return tl::make_unexpected(std::variant<std::invalid_argument, std::out_of_range>{ ex });
  }
}

} // namespace util

max_stars_repo_name: igd-geo/schwarzwald

max_issues_repo_name: igd-geo/schwarzwald

max_forks_repo_name: igd-geo/schwarzwald


text:
#ifndef _SPC_RESP_H

#define _SPC_RESP_H

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

#include <gsl/gsl_spline.h>
#include <gsl/gsl_nan.h>
#include <gsl/gsl_sys.h>

#include "fitsio.h"
#include "aXe_grism.h"
#include "spc_cfg.h"
#include "spc_spc.h"

#define RESPBUFFERSIZE 1024

// number of points used
// for smoothing
#define RESP_SMOOTH_LENGTH 50

// the number of sigmas
// to smooth over
#define RESP_SMOOTH_NSIG 3.0

#define MAX_ITER_TL 1000

/*
 * Struct: interpolator
 */
typedef struct
{
  double xmin;
  double xmax;
  interpolator *resp_values;
  interpolator *resp_errors;
}
  response_function;

typedef struct
{
  double wavelength;
  const calib_function *wl_calibration;
}
trlength_search;

extern spectrum *
get_response_function_from_FITS(char filename[], int hdunum);

extern void
apply_response_function(spectrum *spec, spectrum *resp, const int quant_cont);

extern void
apply_response_functionII(spectrum *spec, response_function *resp_func, const int quant_cont);

extern d_point
get_smooth_pars_for_beam(const aperture_conf *conf, const int smooth_conv, beam actbeam);

extern int
check_conf_for_smoothing(const aperture_conf *conf, const int smooth_conv);

extern void
apply_smoothed_response(const calib_function *wl_calibration, const int for_grism,
                        const int quant_cont, response_function *resp_func,
                        const d_point smooth_pars, spectrum *spec);

extern double
find_wavelength(double x, void *params);

extern double
get_tlength_prismwav(const double wavelength, const calib_function *wl_calibration);

extern double
get_central_tracelength(const double wavelength, const calib_function *wl_calibration,
                        const int for_grism);

extern void
get_smoothed_response(const double wavelength, const d_point smooth_pars,
                           const calib_function *wl_calibration, const int for_grism,
                           const gsl_vector *weights, response_function *resp_func,
                           double *resp_vals);
//extern void
//get_smoothed_response(const double wavelength, const double sigma_wav, const gsl_vector *weights,
//                     response_function *resp_func, double *resp_vals);

extern void
get_weighted_sensitivity(const gsl_vector *pixvalues, const gsl_vector *errvalues, const gsl_vector *weights,
                         const gsl_vector *pmask, double *resp_vals);

extern void
fill_weight(gsl_vector *weights);

extern void
get_troughput_table_name(char *filename, int beamID, char *table_name);

extern double
get_response_value_plus(const spectrum *resp, const double wavelength,
                        int *nguess);

extern void
get_response_values(response_function *resp_func, double wavelength, double* resp_vals);

extern  response_function *
create_response_function(char *filename);

extern  void
free_response_function(response_function *resp_func);

extern calib_function *
get_calfunc_for_beam(const beam actbeam, const int for_grism, char CONF_file[], const aperture_conf * conf);

#endif

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe


text:
/* multifit/multiwlinear.c
 *
 * Copyright (C) 2015 Patrick Alken
 *
 * 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 <gsl/gsl_errno.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>

#include "linear_common.c"

int
gsl_multifit_wlinear (const gsl_matrix * X,
                      const gsl_vector * w,
                      const gsl_vector * y,
                      gsl_vector * c,
                      gsl_matrix * cov,
                      double *chisq, gsl_multifit_linear_workspace * work)
{
  size_t rank;
  int status = gsl_multifit_wlinear_tsvd(X, w, y, GSL_DBL_EPSILON, c, cov, chisq, &rank, work);

  return status;
}

int
gsl_multifit_wlinear_tsvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           const double tol,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double * chisq,
                           size_t * rank,
                           gsl_multifit_linear_workspace * work)
{
  const size_t n = X->size1;
  const size_t p = X->size2;

  if (y->size != n)
    {
      GSL_ERROR("number of observations in y does not match matrix",
                GSL_EBADLEN);
    }
  else if (w->size != n)
    {
      GSL_ERROR("number of weights in w does not match matrix",
                GSL_EBADLEN);
    }
  else if (p != c->size)
    {
      GSL_ERROR ("number of parameters c does not match matrix",
                 GSL_EBADLEN);
    }
  else if (tol <= 0)
    {
      GSL_ERROR ("tolerance must be positive", GSL_EINVAL);
    }
  else
    {
      int status;
      double rnorm, snorm;
      gsl_matrix_view A = gsl_matrix_submatrix(work->A, 0, 0, n, p);
      gsl_vector_view b = gsl_vector_subvector(work->t, 0, n);

      /* compute A = sqrt(W) X, b = sqrt(W) y */
      status = gsl_multifit_linear_applyW(X, w, y, &A.matrix, &b.vector);
      if (status)
        return status;

      /* compute SVD of A */
      status = gsl_multifit_linear_bsvd(&A.matrix, work);
      if (status)
        return status;

      status = multifit_linear_solve(X, &b.vector, tol, 0.0, rank,
                                     c, &rnorm, &snorm, work);
      if (status)
        return status;

      *chisq = rnorm * rnorm;

      /* variance-covariance matrix cov = s2 * (Q S^-1) (Q S^-1)^T */
      {
        const size_t p = X->size2;
        size_t i, j;
        gsl_matrix_view QSI = gsl_matrix_submatrix(work->QSI, 0, 0, p, p);
        gsl_vector_view D = gsl_vector_subvector(work->D, 0, p);

        for (i = 0; i < p; i++)
          {
            gsl_vector_view row_i = gsl_matrix_row (&QSI.matrix, i);
            double d_i = gsl_vector_get (&D.vector, i);

            for (j = i; j < p; j++)
              {
                gsl_vector_view row_j = gsl_matrix_row (&QSI.matrix, j);
                double d_j = gsl_vector_get (&D.vector, j);
                double s;

                gsl_blas_ddot (&row_i.vector, &row_j.vector, &s);

                gsl_matrix_set (cov, i, j, s / (d_i * d_j));
                gsl_matrix_set (cov, j, i, s / (d_i * d_j));
              }
          }
      }
    }

  return GSL_SUCCESS;
}

int
gsl_multifit_wlinear_svd (const gsl_matrix * X,
                          const gsl_vector * w,
                          const gsl_vector * y,
                          double tol,
                          size_t * rank,
                          gsl_vector * c,
                          gsl_matrix * cov,
                          double *chisq, gsl_multifit_linear_workspace * work)
{
  int status = gsl_multifit_wlinear_tsvd(X, w, y, tol, c, cov, chisq, rank, work);
  return status;
}

int
gsl_multifit_wlinear_usvd (const gsl_matrix * X,
                           const gsl_vector * w,
                           const gsl_vector * y,
                           double tol,
                           size_t * rank,
                           gsl_vector * c,
                           gsl_matrix * cov,
                           double *chisq, gsl_multifit_linear_workspace * work)
{
  /* FIXME: this call does actually perform balancing, but this function is deprecated anyway */
  int status = gsl_multifit_wlinear_tsvd(X, w, y, tol, c, cov, chisq, rank, work);
  return status;
}

ext: c

max_stars_repo_name: ielomariala/Hex-Game

max_issues_repo_name: ielomariala/Hex-Game

max_forks_repo_name: karanbirsandhu/nu-sense


text:
/* cdf/gsl_cdf.h
 *
 * Copyright (C) 2002 Jason H. Stover.
 *
 * 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:  J. Stover */

#ifndef __GSL_CDF_H__
#define __GSL_CDF_H__

#include <gsl/gsl_types.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS          /* empty */
# define __END_DECLS            /* empty */
#endif

__BEGIN_DECLS

GSL_EXPORT double gsl_cdf_ugaussian_P (const double x);
GSL_EXPORT double gsl_cdf_ugaussian_Q (const double x);

GSL_EXPORT double gsl_cdf_ugaussian_Pinv (const double P);
GSL_EXPORT double gsl_cdf_ugaussian_Qinv (const double Q);

GSL_EXPORT double gsl_cdf_gaussian_P (const double x, const double sigma);
GSL_EXPORT double gsl_cdf_gaussian_Q (const double x, const double sigma);

GSL_EXPORT double gsl_cdf_gaussian_Pinv (const double P, const double sigma);
GSL_EXPORT double gsl_cdf_gaussian_Qinv (const double Q, const double sigma);

GSL_EXPORT double gsl_cdf_gamma_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_gamma_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_gamma_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_gamma_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_cauchy_P (const double x, const double a);
GSL_EXPORT double gsl_cdf_cauchy_Q (const double x, const double a);

GSL_EXPORT double gsl_cdf_cauchy_Pinv (const double P, const double a);
GSL_EXPORT double gsl_cdf_cauchy_Qinv (const double Q, const double a);

GSL_EXPORT double gsl_cdf_laplace_P (const double x, const double a);
GSL_EXPORT double gsl_cdf_laplace_Q (const double x, const double a);

GSL_EXPORT double gsl_cdf_laplace_Pinv (const double P, const double a);
GSL_EXPORT double gsl_cdf_laplace_Qinv (const double Q, const double a);

GSL_EXPORT double gsl_cdf_rayleigh_P (const double x, const double sigma);
GSL_EXPORT double gsl_cdf_rayleigh_Q (const double x, const double sigma);

GSL_EXPORT double gsl_cdf_rayleigh_Pinv (const double P, const double sigma);
GSL_EXPORT double gsl_cdf_rayleigh_Qinv (const double Q, const double sigma);

GSL_EXPORT double gsl_cdf_chisq_P (const double x, const double nu);
GSL_EXPORT double gsl_cdf_chisq_Q (const double x, const double nu);

GSL_EXPORT double gsl_cdf_chisq_Pinv (const double P, const double nu);
GSL_EXPORT double gsl_cdf_chisq_Qinv (const double Q, const double nu);

GSL_EXPORT double gsl_cdf_exponential_P (const double x, const double mu);
GSL_EXPORT double gsl_cdf_exponential_Q (const double x, const double mu);

GSL_EXPORT double gsl_cdf_exponential_Pinv (const double P, const double mu);
GSL_EXPORT double gsl_cdf_exponential_Qinv (const double Q, const double mu);

GSL_EXPORT double gsl_cdf_exppow_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_exppow_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_tdist_P (const double x, const double nu);
GSL_EXPORT double gsl_cdf_tdist_Q (const double x, const double nu);

GSL_EXPORT double gsl_cdf_tdist_Pinv (const double P, const double nu);
GSL_EXPORT double gsl_cdf_tdist_Qinv (const double Q, const double nu);

GSL_EXPORT double gsl_cdf_fdist_P (const double x, const double nu1, const double nu2);
GSL_EXPORT double gsl_cdf_fdist_Q (const double x, const double nu1, const double nu2);

GSL_EXPORT double gsl_cdf_fdist_Pinv (const double P, const double nu1, const double nu2);
GSL_EXPORT double gsl_cdf_fdist_Qinv (const double Q, const double nu1, const double nu2);

GSL_EXPORT double gsl_cdf_beta_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_beta_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_beta_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_beta_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_flat_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_flat_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_flat_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_flat_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_lognormal_P (const double x, const double zeta, const double sigma);
GSL_EXPORT double gsl_cdf_lognormal_Q (const double x, const double zeta, const double sigma);

GSL_EXPORT double gsl_cdf_lognormal_Pinv (const double P, const double zeta, const double sigma);
GSL_EXPORT double gsl_cdf_lognormal_Qinv (const double Q, const double zeta, const double sigma);

GSL_EXPORT double gsl_cdf_gumbel1_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_gumbel1_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_gumbel1_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_gumbel1_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_gumbel2_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_gumbel2_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_gumbel2_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_gumbel2_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_weibull_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_weibull_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_weibull_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_weibull_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_pareto_P (const double x, const double a, const double b);
GSL_EXPORT double gsl_cdf_pareto_Q (const double x, const double a, const double b);

GSL_EXPORT double gsl_cdf_pareto_Pinv (const double P, const double a, const double b);
GSL_EXPORT double gsl_cdf_pareto_Qinv (const double Q, const double a, const double b);

GSL_EXPORT double gsl_cdf_logistic_P (const double x, const double a);
GSL_EXPORT double gsl_cdf_logistic_Q (const double x, const double a);

GSL_EXPORT double gsl_cdf_logistic_Pinv (const double P, const double a);
GSL_EXPORT double gsl_cdf_logistic_Qinv (const double Q, const double a);

GSL_EXPORT double gsl_cdf_binomial_P (const unsigned int k, const double p, const unsigned int n);
GSL_EXPORT double gsl_cdf_binomial_Q (const unsigned int k, const double p, const unsigned int n);

GSL_EXPORT double gsl_cdf_poisson_P (const unsigned int k, const double mu);
GSL_EXPORT double gsl_cdf_poisson_Q (const unsigned int k, const double mu);

GSL_EXPORT double gsl_cdf_geometric_P (const unsigned int k, const double p);
GSL_EXPORT double gsl_cdf_geometric_Q (const unsigned int k, const double p);

GSL_EXPORT double gsl_cdf_negative_binomial_P (const unsigned int k, const double p, const double n);
GSL_EXPORT double gsl_cdf_negative_binomial_Q (const unsigned int k, const double p, const double n);

GSL_EXPORT double gsl_cdf_pascal_P (const unsigned int k, const double p, const unsigned int n);
GSL_EXPORT double gsl_cdf_pascal_Q (const unsigned int k, const double p, const unsigned int n);

GSL_EXPORT double gsl_cdf_hypergeometric_P (const unsigned int k, const unsigned int n1,
                                            const unsigned int n2, const unsigned int t);
GSL_EXPORT double gsl_cdf_hypergeometric_Q (const unsigned int k, const unsigned int n1,
                                            const unsigned int n2, const unsigned int t);

__END_DECLS

#endif /* __GSL_CDF_H__ */

max_stars_repo_name: dynaryu/vaws

max_issues_repo_name: dynaryu/vaws

max_forks_repo_name: dynaryu/vaws


text:
/*
 * Copyright (c) 2020 Microsoft. All rights reserved.
 */

#pragma once

#include <array>

#include <algorithm>
#include <gsl/span>

#include "common/common.h"
#include "common/logger.h"

// A wrapper around std::array that is specialized to uint8_t and enforces
// initialization
template <size_t TLength> class Blob : public CanBeSerialsed {
  using CArrayType = uint8_t[TLength];
  using CppArrayType = std::array<uint8_t, TLength>;

  CppArrayType _data;

public:
  static constexpr size_t array_size = TLength;
  using array_type = CppArrayType;

  constexpr Blob(const CArrayType &data) {
    std::copy(std::begin(data), std::end(data), _data.begin());
  }

  constexpr Blob(const gsl::span<const uint8_t> &data) {
    const size_t size = std::min(_data.size(), data.size());
    std::copy_n(data.begin(), size, _data.begin());
  }

  constexpr Blob() noexcept : _data() {
    for (auto &b : _data) {
      b = 0;
    }
  }

  CppArrayType &array() { return _data; }

  const CppArrayType &array() const { return _data; }

  MemRef span() { return MakeMemRef(&_data); }

  CMemRef span() const { return MakeCMemRef(&_data); }

  uint8_t *data() noexcept { return _data.data(); }

  const uint8_t *data() const noexcept { return _data.data(); }

  uint8_t *end() {
    const MemRef span = gsl::make_span(data(), TLength + 1);
    const MemRef sub = span.subspan(TLength, 1);
    return sub.data();
  }

  size_t length() const { return TLength; }

  bool operator==(const Blob<TLength> &d) const { return d.array() == array(); }

  bool operator!=(const Blob<TLength> &d) const { return d.array() != array(); }

  // we cannot check the whole object, but we can check the array.
  static_assert(OkAsParam<CppArrayType, TLength>::value);
};

// A wrapper around std::array that is specialized to char and enforces
// initialization
template <size_t TLength> class CharBlob : public CannotBeSerialsed {
  using CArrayType = char[TLength];
  using CppArrayType = std::array<char, TLength>;

private:
  CppArrayType _data;
  size_t _length; // Danger!

public:
  constexpr CharBlob(const gsl::span<const char> &data) {
    const size_t size = data.size();

    if (size > TLength) {
      HW_LOG_FAIL("Blob not long enough");
      while (1)
        ;
    }

    _data.fill(0);
    std::copy_n(data.begin(), size, _data.begin());
    _length = size;
  }

  constexpr CharBlob(const char *data, size_t length) {
    if (length > TLength) {
      HW_LOG_FAIL("Blob not long enough");
      while (1)
        ;
    }

    const CSignedMemRef aspan = MakeCSignedMemRef(data, length);
    _data.fill(0);
    std::copy_n(aspan.begin(), length, _data.begin());
    _length = length;
  }

  constexpr CharBlob() noexcept : _data() {
    _data.fill(0);
    _length = 0;
  }

  CppArrayType &array() { return _data; }

  const CppArrayType &array() const { return _data; }

  SignedMemRef span() { return MakeSignedMemRef(&_data); }

  CSignedMemRef span() const { return MakeCSignedMemRef(&_data); }

  char *data() { return _data.data(); }

  const char *data() const { return _data.data(); }

  size_t maxlength() const { return TLength; }

  size_t length() const { return _length; }

  bool operator==(const CharBlob<TLength> &d) {
    if (d.length() != length()) {
      return false;
    }

    // not using the strcmp as these buffers may not have null-terminating
    // character
    const size_t len = length();
    const CSignedMemRef _dspan = MakeCSignedMemRef(data(), len);
    const CSignedMemRef dspan = MakeCSignedMemRef(d.data(), len);

    return (_dspan == dspan);
  }

  bool operator!=(const CharBlob<TLength> &d) {
    if (d.length() != length()) {
      return true;
    }

    // not using the strcmp as these buffers may not have null-terminating
    // character
    const size_t len = length();
    const CSignedMemRef _dspan = MakeCSignedMemRef(data(), len);
    const CSignedMemRef dspan = MakeCSignedMemRef(d.data(), len);

    return (_dspan != dspan);
  }
};


text:
/*! \file getdata.h
    \brief rのメッシュと、そのメッシュにおける電子密度を与えるクラスの宣言

    Copyright © 2015-2019 @dc1394 All Rights Reserved.
    This software is released under the BSD 2-Clause License.
*/

#ifndef _GETRMESHANDRHO_H_
#define _GETRMESHANDRHO_H_

#pragma once

#include "../utility/property.h"
#include <cstdint>          // for std::int32_t, std::uint32_t
#include <memory>           // for std::unique_ptr
#include <string>           // for std::string
#include <gsl/gsl_spline.h> // for gsl_interp_accel, gsl_interp_accel_free, gsl_spline, gsl_spline_free

namespace getdata {
    using namespace utility;

    //! A class.
    /*!
        rのメッシュと、そのメッシュにおける動径波動関数を与えるクラス
    */
    class GetData final {
        // #region 列挙型

    public:
        //!  A enumerated type
        /*!
            密度か動径波動関数かを表す列挙型
        */
        enum class Rho_Wf_type {
            // 密度
            RHO,
            // 動径波動関数
            WF
        };

        // #endregion 列挙型

        // #region コンストラクタ・デストラクタ

        //! A constructor.
        /*!
            唯一のコンストラクタ
            \param filename rのメッシュと、そのメッシュにおける電子密度が記録されたデータファイル名
        */
        GetData(std::string const & filename);

        //! A destructor.
        /*!
            デフォルトデストラクタ
        */
        ~GetData() = default;

        // #endregion コンストラクタ・デストラクタ

        // #region メンバ関数

        //!  A public member function (const).
        /*!
            関数の値を返す
            \param r rの値
            \return 関数の値
        */
        double operator()(double r) const;

        //!  A public member function (const).
        /*!
            関数の微分の値を返す
            \param r rの値
            \return 関数の微分の値
        */
        double dphidr(double r) const;

        // #endregion メンバ関数

        // #region プロパティ

        //! A property.
        /*!
            元素名
        */
        Property<std::string const&> Atomname;

        //!  A property.
        /*!
            方位量子数へのプロパティ
        */
        Property<std::uint32_t> const L;

        //!  A property.
        /*!
            主量子数へのプロパティ
        */
        Property<std::int32_t> const N;

        //!  A property.
        /*!
            軌道へのプロパティ
        */
        Property<std::string> const Orbital;

        //! A property.
        /*!
            波動関数の最大値のプロパティ
        */
        Property<double> const Phimax;

        //! A property.
        /*!
            波動関数が最大値を取るときのrへのプロパティ
        */
        Property<double> const R2rhomaxr;

        //!  A private member variable.
        /*!
            解く方程式のタイプへのプロパティ
        */
        Property<GetData::Rho_Wf_type> const Rho_wf_type;

        //! A property.
        /*!
            rのメッシュの最大値のプロパティ
        */
        Property<double> const R_meshmax;

        //! A property.
        /*!
            rのメッシュの最小値のプロパティ
        */
        Property<double> const R_meshmin;

        // #endregion プロパティ

        // #region メンバ変数

    private:
        //! A private member variable.
        /*!
            gsl_interp_accelへのスマートポインタ
        */
        std::unique_ptr<gsl_interp_accel, decltype(&gsl_interp_accel_free)> const acc_;

        //!  A private member variable.
        /*!
            元素名
        */
        std::string atomname_;

        //!  A private member variable.
        /*!
            方位量子数
        */
        std::uint32_t l_;

        //!  A private member variable.
        /*!
            主量子数
        */
        std::int32_t n_;

        //!  A private member variable.
        /*!
            軌道
        */
        std::string orbital_;

        //!  A private member variable.
        /*!
            波動関数φ
        */
        std::vector<double> phi_;

        //!  A private member variable.
        /*!
            波動関数の最大値
        */
        double phimax_;

        //!  A private member variable.
        /*!
            波動関数が最大値を取るときのr
        */
        double r2rhomaxr_;

        //!  A private member variable.
        /*!
            解く方程式のタイプ
        */
        Rho_Wf_type rho_wf_type_;

        //!  A private member variable.
        /*!
            rのメッシュの最大値
        */
        double r_meshmax_;

        //!  A private member variable.
        /*!
            rのメッシュの最小値
        */
        double r_meshmin_;

        //!  A private member variable.
        /*!
           rのメッシュ
        */
        std::vector<double> r_mesh_;

        //! A private member variable.
        /*!
            gsl_interp_typeへのスマートポインタ
        */
        std::unique_ptr<gsl_spline, decltype(&gsl_spline_free)> spline_;

        // #endregion メンバ変数

        // #region 禁止されたコンストラクタ・メンバ関数
    public:

        //! A public constructor (deleted).
        /*!
            デフォルトコンストラクタ（禁止）
        */
        GetData() = delete;

        //! A public copy constructor (deleted).
        /*!
            コピーコンストラクタ（禁止）
            \param dummy コピー元のオブジェクト（未使用）
        */
        GetData(GetData const & dummy) = delete;

        //! A public member function (deleted).
        /*!
            operator=()の宣言（禁止）
            \param dummy コピー元のオブジェクト（未使用）
            \return コピー元のオブジェクト
        */
        GetData& operator=(GetData const & dummy) = delete;

        // #endregion 禁止されたコンストラクタ・メンバ関数
    };

    // #region メンバ関数

    inline double GetData::operator()(double r) const
    {
        return gsl_spline_eval(spline_.get(), r, acc_.get());
    }

    inline double GetData::dphidr(double r) const
    {
        return gsl_spline_eval_deriv(spline_.get(), r, acc_.get());
    }

    // #endregion メンバ関数
}

#endif  // _GETRMESHANDRHO_H_


text:
/**
 * File: trace_conf.h
 */

#ifndef _TRACE_CONF_H
#define _TRACE_CONF_H

#include <string.h>
#include <math.h>
#include <gsl/gsl_vector.h>
#include "aXe_grism.h"
#include "aXe_errors.h"
#include "spc_cfg.h"
#include "disp_conf.h"

/**
	A structure to contain the coefficient of a 2D
	polynomial that can be used to compute the given
	coefficient of the polynomial dispersion relation
	as at a particular i.j location in the image
*/
typedef struct
{
     gsl_vector *pol;		/* A vector containing the 2D polynomial coefficients */
     d_point offset;		/* X and Y offsets to apply to input coordinates */
     d_point cpoint;		/* The detector location at which this structure was computed */
     int ID;			/* The ID of the beam for which this coefficient is defined */
     char file[MAXCHAR];
}
tracestruct;

extern int
get_beam_trace_norder (char *filename, int beamID);

extern gsl_vector *
get_beam_trace_order (char *filename, int beamID, int order);

extern float
get_trace_coeff_at_pos (char *filename, int beamID, int order, d_point p);

extern gsl_vector *
get_trace_coeffs_at_pos (char *filename, int beamID, d_point p);

extern gsl_vector *
get_beam_trace_xoff (char *filename, int beamID);

extern gsl_vector *
get_beam_trace_yoff (char *filename, int beamID);

extern float
get_trace_xoff_at_pos (char *filename, int beamID, d_point p);

extern float
get_trace_yoff_at_pos (char *filename, int beamID, d_point p);

extern float
eval_trace_off_at_pos (gsl_vector *coeffs, d_point p, int beamID);

extern tracestruct *
get_tracestruct_at_pos (char *filename, int beamID, d_point p);

extern void
tracestruct_fprintf (FILE * file, tracestruct * trace);

extern void
free_tracestruct (tracestruct * trace);
#endif


text:
/* AUTORIGHTS
Copyright (C) 2007 Princeton University

This file is part of Ferret Toolkit.

Ferret Toolkit 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, 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 <assert.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <fenv.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_cdf.h>
#include <gsl/gsl_vector.h>
#include "LSH.h"
#include "perturb.h"

#define MISS_GAMMA	1

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif

#define MAX_M	64
#define MAX_T	600

static inline double p_col_helper (double x)
{
	double result;
	result = 2.0*gsl_cdf_ugaussian_P(x) - 1.0;
	result += sqrt(2.0/M_PI) * (exp(-x*x/2.0)-1.0)/x;
	return result;
}

static inline double p_col_p_helper (double x, double k)
{
	return gsl_cdf_ugaussian_P((1.0 + k) * x) - gsl_cdf_ugaussian_P(k * x);
}

int LSH_recall_init (LSH_recall_t *recall, int d_step, float d_min, float d_max,
		int M, int L, int T, float W)
{
	ptb_vec_t *ptb_score;
	ptb_vec_t *ptb_set;
	float **table;
	double dist, r, p;
	int i, j, k;
	int l;

	ptb_score = gen_score(M);
	assert(ptb_score != NULL);
	ptb_set = type_calloc(ptb_vec_t, T);
	assert(ptb_set != NULL);
	gen_perturb_set(ptb_score, ptb_set, M, T);

	recall->d_step = d_step;
	recall->d_min = d_min;
	recall->d_max = d_max;
	recall->T = T;
	recall->table = table = type_matrix_alloc(float, T + 1, d_step);

	for (i = 0; i < d_step; i++)
	{
		dist = d_min + (d_max - d_min) * i / d_step;
		dist = W / sqrt(dist);
		p = p_col_helper(dist);
		table[0][i] = log(1.0 - pow(p, M));
		for (j = 1; j <= T; j++)
		{
			r = 1.0;
			l = 0;
			for (k = 0; k < 2 * M; k++)
			if (ptb_set[j-1].set & (1 << k))
			{
				r *= p_col_p_helper(dist, ptb_score[k].key1);
				l++;
			}
			assert(l > 0);
			l = M - l;
			r *= pow(p, l);
			table[j][i] = table[j-1][i] + log(1.0-r);
		}
		for (j = 0; j <= T; j++)
		{
			table[j][i] = 1.0 - exp(L * table[j][i]);
		}
	}

	free(ptb_score);
	free(ptb_set);
	return 0;
}

int LSH_recall_dump (LSH_recall_t *recall, CASS_FILE *fout)
{
	int ret;
	ret = cass_write_size(&recall->d_step, 1, fout);
	ret += cass_write_size(&recall->T, 1, fout);
	ret += cass_write_float(&recall->d_min, 1, fout);
	ret += cass_write_float(&recall->d_max, 1, fout);

	if (ret != 4) return CASS_ERR_IO;

	if (type_matrix_dump_stream(float, fout, recall->T + 1, recall->d_step, recall->table) != 0) return CASS_ERR_IO;

	return 0;
}

int LSH_recall_load (LSH_recall_t *recall, CASS_FILE *fin)
{
	int ret;
	cass_size_t row, col;
	ret = cass_read_size(&recall->d_step, 1, fin);
	ret += cass_read_size(&recall->T, 1, fin);
	ret += cass_read_float(&recall->d_min, 1, fin);
	ret += cass_read_float(&recall->d_max, 1, fin);

	if (ret != 4) return CASS_ERR_IO;

	type_matrix_load_stream(float, fin, &row, &col, &recall->table);

	assert(row == recall->T +1);
	assert(col == recall->d_step);

	return 0;
}

max_stars_repo_name: OleksiiOleksenko/fex

max_issues_repo_name: OleksiiOleksenko/fex

max_forks_repo_name: OleksiiOleksenko/fex